分享matlab冲击响应谱程序

ChaChing

Ref : http://www.mathworks.com/matlabc ... ctrum/content/srs.m

1. disp(' ')
   
2. disp(' srs.m ver 2.0 July 3, 2006')
   
3. disp(' by Tom Irvine Email: tomirvine@aol.com')
   
4. disp(' ')
   
5. disp(' This program calculates the shock response spectrum')
   
6. disp(' of an acceleration time history, which is pre-loaded into Matlab.')
   
7. disp(' The time history must have two columns: time(sec) & acceleration')
   
8. disp(' ')
   
9. %
   
10. clear t;
    
11. clear y;
    
12. clear yy;
    
13. clear n;
    
14. clear fn;
    
15. clear a1;
    
16. clear a2
    
17. clear b1;
    
18. clear b2;
    
19. clear jnum;
    
20. clear THM;
    
21. clear resp;
    
22. clear x_pos;
    
23. clear x_neg;
    
24. %
    
25. iunit=input(' Enter acceleration unit: 1= G 2= m/sec^2 ');
    
26. %
    
27. disp(' ')
    
28. disp(' Select file input method ');
    
29. disp(' 1=external ASCII file ');
    
30. disp(' 2=file preloaded into Matlab ');
    
31. file_choice = input('');
    
32. %
    
33. if(file_choice==1)
    
34. [filename, pathname] = uigetfile('*.*');
    
35. filename = fullfile(pathname, filename);
    
36. %
    
37. fid = fopen(filename,'r');
    
38. THM = fscanf(fid,'%g %g',[2 inf]);
    
39. THM=THM';
    
40. else
    
41. THM = input(' Enter the matrix name: ');
    
42. end
    
43. %
    
44. t=double(THM(:,1));
    
45. y=double(THM(:,2));
    
46. %
    
47. tmx=max(t);
    
48. tmi=min(t);
    
49. n = length(y);
    
50. %
    
51. out1 = sprintf('\n %d samples \n',n);
    
52. disp(out1)
    
53. %
    
54. dt=(tmx-tmi)/(n-1);
    
55. sr=1./dt;
    
56. %
    
57. out1 = sprintf(' SR = %g samples/sec dt = %g sec \n',sr,dt);
    
58. disp(out1)
    
59. %
    
60. fn(1)=input(' Enter the starting frequency (Hz) ');
    
61. if fn(1)>sr/30.
    
62. fn(1)=sr/30.;
    
63. end
    
64. %
    
65. idamp=input(' Enter damping format: 1= damping ratio 2= Q ');
    
66. %
    
67. disp(' ')
    
68. if(idamp==1)
    
69. damp=input(' Enter damping ratio (typically 0.05) ');
    
70. else
    
71. Q=input(' Enter the amplification factor (typically Q=10) ');
    
72. damp=1./(2.*Q);
    
73. end
    
74. %
    
75. disp(' ')
    
76. disp(' Select algorithm: ')
    
77. disp(' 1=Kelly-Richman 2=Smallwood ');
    
78. ialgorithm=input(' ');
    
79. %
    
80. tmax=(tmx-tmi) + 1./fn(1);
    
81. limit = round( tmax/dt );
    
82. n=limit;
    
83. yy=zeros(1,limit);
    
84. for i=1:length(y)
    
85. yy(i)=y(i);
    
86. end
    
87. %
    
88. disp(' ')
    
89. disp(' Calculating response..... ')
    
90. %
    
91. % SRS engine
    
92. %
    
93. for j=1:1000
    
94. %
    
95. omega=2.*pi*fn(j);
    
96. omegad=omega*sqrt(1.-(damp^2));
    
97. cosd=cos(omegad*dt);
    
98. sind=sin(omegad*dt);
    
99. domegadt=damp*omega*dt;
    
100. %
     
101. if(ialgorithm==1)
     
102. a1(j)=2.*exp(-domegadt)*cosd;
     
103. a2(j)=-exp(-2.*domegadt);
     
104. b1(j)=2.*domegadt;
     
105. b2(j)=omega*dt*exp(-domegadt);
     
106. b2(j)=b2(j)*( (omega/omegad)*(1.-2.*(damp^2))*sind -2.*damp*cosd );
     
107. b3(j)=0;
     
108. %
     
109. else
     
110. E=exp(-damp*omega*dt);
     
111. K=omegad*dt;
     
112. C=E*cos(K);
     
113. S=E*sin(K);
     
114. Sp=S/K;
     
115. %
     
116. a1(j)=2*C;
     
117. a2(j)=-E^2;
     
118. b1(j)=1.-Sp;
     
119. b2(j)=2.*(Sp-C);
     
120. b3(j)=E^2-Sp;
     
121. end
     
122. forward=[ b1(j), b2(j), b3(j) ];
     
123. back =[ 1, -a1(j), -a2(j) ];
     
124. %
     
125. resp=filter(forward,back,yy);
     
126. %
     
127. x_pos(j)= max(resp);
     
128. x_neg(j)= min(resp);
     
129. %
     
130. jnum=j;
     
131. if fn(j) > sr/8.
     
132. break
     
133. end
     
134. fn(j+1)=fn(1)*(2. ^ (j*(1./12.)));
     
135. end
     
136. %
     
137. % Output options
     
138. %
     
139. disp(' ')
     
140. disp(' Select output option ');
     
141. choice=input(' 1=plot only 2=plot & output text file ' );
     
142. disp(' ')
     
143. %
     
144. if choice == 2
     
145. %%
     
146. [writefname, writepname] = uiputfile('*','Save SRS data as');
     
147. writepfname = fullfile(writepname, writefname);
     
148. writedata = [fn' x_pos' (abs(x_neg))' ];
     
149. fid = fopen(writepfname,'w');
     
150. fprintf(fid,' %g %g %g\n',writedata');
     
151. fclose(fid);
     
152. %%
     
153. % disp(' Enter output filename ');
     
154. % SRS_filename = input(' ','s');
     
155. %
     
156. % fid = fopen(SRS_filename,'w');
     
157. % for j=1:jnum
     
158. % fprintf(fid,'%7.2f %10.3f %10.3f \n',fn(j),x_pos(j),abs(x_neg(j)));
     
159. % end
     
160. % fclose(fid);
     
161. end
     
162. %
     
163. % Plot SRS
     
164. %
     
165. disp(' ')
     
166. disp(' Plotting output..... ')
     
167. %
     
168. % Find limits for plot
     
169. %
     
170. srs_max = max(x_pos);
     
171. if max( abs(x_neg) ) > srs_max
     
172. srs_max = max( abs(x_neg ));
     
173. end
     
174. srs_min = min(x_pos);
     
175. if min( abs(x_neg) ) < srs_min
     
176. srs_min = min( abs(x_neg ));
     
177. end
     
178. %
     
179. figure(1);
     
180. plot(fn,x_pos,fn,abs(x_neg),'-.');
     
181. %
     
182. if iunit==1
     
183. ylabel('Peak Accel (G)');
     
184. else
     
185. ylabel('Peak Accel (m/sec^2)');
     
186. end
     
187. xlabel('Natural Frequency (Hz)');
     
188. Q=1./(2.*damp);
     
189. out5 = sprintf(' Acceleration Shock Response Spectrum Q=%g ',Q);
     
190. title(out5);
     
191. grid;
     
192. set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','log');
     
193. legend ('positive','negative',2);
     
194. %
     
195. ymax= 10^(round(log10(srs_max)+0.8));
     
196. ymin= 10^(round(log10(srs_min)-0.6));
     
197. %
     
198. fmax=max(fn);
     
199. fmin=fmax/10.;
     
200. %
     
201. fmax= 10^(round(log10(fmax)+0.5));
     
202. %
     
203. if fn(1) >= 0.1
     
204. fmin=0.1;
     
205. end
     
206. if fn(1) >= 1
     
207. fmin=1;
     
208. end
     
209. if fn(1) >= 10
     
210. fmin=10;
     
211. end
     
212. if fn(1) >= 100
     
213. fmin=100;
     
214. end
     
215. axis([fmin,fmax,ymin,ymax]);
     
216. %
     
217. disp(' ')
     
218. disp(' Plot pseudo velocity? ');
     
219. vchoice=input(' 1=yes 2=no ' );
     
220. if(vchoice==1)
     
221. figure(2);
     
222. %
     
223. % Convert to pseudo velocity
     
224. %
     
225. for j=1:jnum
     
226. if iunit==1
     
227. x_pos(j)=386.*x_pos(j)/(2.*pi*fn(j));
     
228. x_neg(j)=386.*x_neg(j)/(2.*pi*fn(j));
     
229. else
     
230. x_pos(j)=x_pos(j)/(2.*pi*fn(j));
     
231. x_neg(j)=x_neg(j)/(2.*pi*fn(j));
     
232. end
     
233. end
     
234. %
     
235. srs_max = max(x_pos);
     
236. if max( abs(x_neg) ) > srs_max
     
237. srs_max = max( abs(x_neg ));
     
238. end
     
239. srs_min = min(x_pos);
     
240. if min( abs(x_neg) ) < srs_min
     
241. srs_min = min( abs(x_neg ));
     
242. end
     
243. %
     
244. plot(fn,x_pos,fn,abs(x_neg),'-.');
     
245. %
     
246. if iunit==1
     
247. ylabel('Velocity (in/sec)');
     
248. else
     
249. ylabel('Velocity (m/sec)');
     
250. end
     
251. xlabel('Natural Frequency (Hz)');
     
252. Q=1./(2.*damp);
     
253. out5 = sprintf(' Pseudo Velocity Shock Response Spectrum Q=%g ',Q);
     
254. title(out5);
     
255. grid;
     
256. set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log','YScale','log');
     
257. legend ('positive','negative',2);
     
258. %
     
259. ymax= 10^(round(log10(srs_max)+0.8));
     
260. ymin= 10^(round(log10(srs_min)-0.6));
     
261. %
     
262. axis([fmin,fmax,ymin,ymax]);
     
263. end

复制代码

牛小贱

先收藏了！！！

jxx19907

谢谢楼主！

王利明

楼主，想问下这个代码原理，是用的Laplace小波来做点乘寻找冲击吗？

ChaChing

王利明 发表于 2016-2-3 20:51
楼主，想问下这个代码原理，是用的Laplace小波来做点乘寻找冲击吗？

水平有限,不清楚Laplace小波,但个人以为没那麼复杂