求助各位大虾冲击响应谱的matlab程序

QQ-WYH

冲击信号是一个时域信号，是做试验用示波器得到的数据，怎么样使用MATLAB计算出其冲击响应谱，画出冲击响应谱的波形。小弟谢过大虾们了！邮箱：764112884@qq.com

犟牛

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

复制代码