matlab 信號(hào)處理

1、一號(hào)樓樓板在環(huán)境激勵(lì)下的振動(dòng)信號(hào)處理振動(dòng)是一種普遍存在的自然現(xiàn)象。為了了解結(jié)構(gòu)對(duì)振動(dòng)的反應(yīng)，我們對(duì)一號(hào)樓樓板進(jìn)行了測(cè)驗(yàn)，本文將對(duì)所測(cè)得的數(shù)據(jù)進(jìn)行處理。具體包括振動(dòng)信號(hào)的預(yù)處理、時(shí)域處理頻域處理、數(shù)字信號(hào)的生成、試驗(yàn)?zāi)B(tài)參數(shù)的頻域識(shí)別和時(shí)域識(shí)別以及整體識(shí)別。本文將對(duì)所測(cè)得的加速度數(shù)據(jù)進(jìn)行處理。一、 振動(dòng)信號(hào)預(yù)處理振動(dòng)信號(hào)預(yù)處理是將振動(dòng)測(cè)試中采集到的數(shù)據(jù)盡可能真實(shí)地還原成實(shí)際振動(dòng)狀況的最基本的數(shù)據(jù)加工方法。下面將對(duì)其進(jìn)行消除多項(xiàng)式趨勢(shì)項(xiàng)處理，采樣數(shù)據(jù)的平滑處理。1、 消除多項(xiàng)式趨勢(shì)項(xiàng)處理在振動(dòng)測(cè)試中采集到的振動(dòng)信號(hào)數(shù)據(jù)，由于放大器隨溫度變化產(chǎn)生的零點(diǎn)漂移、傳感器頻率范圍外低頻性能的不穩(wěn)定以及傳感器

2、周圍的環(huán)境干擾，往往會(huì)偏離基線，甚至偏離基線的大小還會(huì)隨時(shí)間變化。偏離基線隨時(shí)間變化的整個(gè)過程被稱為信號(hào)的趨勢(shì)項(xiàng)。趨勢(shì)項(xiàng)直接影響信號(hào)的正確性，應(yīng)將其去除。本文將借助MATLAB進(jìn)行最小二乘法消除多項(xiàng)式趨勢(shì)項(xiàng)。a、最小二乘法消除多項(xiàng)式趨勢(shì)項(xiàng)的程序 %clear % 清除內(nèi)存中所有變量和函數(shù)clc % 清除工作窗口中所顯示的內(nèi)容close all hidden % 關(guān)閉所有隱藏的窗口%提示用鍵盤輸入輸入數(shù)據(jù)文件名fni=input('消除多項(xiàng)式趨勢(shì)項(xiàng)-輸入數(shù)據(jù)文件名:','s'); %以只讀方式打開數(shù)據(jù)文件fid=fopen(fni,'r'); sf

3、 = fscanf(fid,'%f',1); %讀入采樣頻率值m = fscanf(fid,'%d',1); %讀入擬合多項(xiàng)式階數(shù)fno = fscanf(fid,'%s',1);%讀入輸出數(shù)據(jù)文件名x = fscanf(fid,'%f',inf);%讀入時(shí)程數(shù)據(jù)存成列向量%關(guān)閉數(shù)據(jù)文件status=fclose(fid); %取信號(hào)數(shù)據(jù)長(zhǎng)度 n=length(x); %建立離散時(shí)間列向量 t=(0:1/sf:(n-1)/sf)' %計(jì)算趨勢(shì)項(xiàng)的多項(xiàng)式待定系數(shù)向量aa=polyfit(t,x,m); %用x減去多項(xiàng)式系數(shù)a

4、生成的趨勢(shì)項(xiàng)y=x-polyval(a,t); %將分成2行1列的圖形窗口的第1列設(shè)為當(dāng)前繪圖區(qū)域subplot(2,1,1); %繪制x對(duì)于t的時(shí)程曲線圖形plot(t,x); %在圖幅上添加坐標(biāo)網(wǎng)格grid on; %將分成2行1列的圖形窗口的第2列設(shè)為當(dāng)前繪圖區(qū)域subplot(2,1,2); %繪制y對(duì)于t的時(shí)程曲線圖形 plot(t,y); %在圖幅上添加坐標(biāo)網(wǎng)格grid on; %以寫的方式打開文件或建立一個(gè)新文件 fid=fopen(fno,'w'); %進(jìn)行n次循環(huán)將計(jì)算結(jié)果寫到輸出數(shù)據(jù)文件中 for k=1:n %每行輸出兩個(gè)實(shí)型數(shù)據(jù)，t為時(shí)間，y為消除趨勢(shì)項(xiàng)

5、后的結(jié)果 fprintf(fid,'%f %fn',t(k),y(k); %循環(huán)體結(jié)束語句 end %關(guān)閉數(shù)據(jù)文件status=fclose(fid); b、處理結(jié)果2、 采樣數(shù)據(jù)的平滑處理通過數(shù)據(jù)采集器采樣的道德振動(dòng)信號(hào)數(shù)據(jù)往往疊加有噪聲信號(hào)。由于隨機(jī)干擾信號(hào)的平帶較寬，有時(shí)高頻成分所占比例很大，使得采集到離散數(shù)據(jù)繪成的振動(dòng)曲線上呈現(xiàn)許多毛刺，很不光滑。為了消除干擾信號(hào)的影響，本文進(jìn)行了五點(diǎn)滑動(dòng)平均法平滑處理。a、五點(diǎn)滑動(dòng)平均法平滑處理的程序%clearclcclose all hidden%fni=input('五點(diǎn)滑動(dòng)平均法平滑處理-輸入數(shù)據(jù)文件名:',&

6、#39;s'); fid=fopen(fni,'r'); sf = fscanf(fid,'%f',1); %采樣頻率 m = fscanf(fid,'%d',1); %平滑次數(shù) fno = fscanf(fid,'%s',1);%輸出數(shù)據(jù)文件名x = fscanf(fid,'%f',inf);%輸入數(shù)據(jù)存成列向量status=fclose(fid); %取信號(hào)數(shù)據(jù)長(zhǎng)度 n=length(x); %建立離散時(shí)間列向量 t=(0:1/sf:(n-1)/sf)' %將x賦值給a a=x; %循環(huán)m次進(jìn)行

7、平滑處理計(jì)算 for k=1:m b(1)=(3*a(1)+2*a(2)+a(3)-a(4)/5; b(2)=(4*a(1)+3*a(2)+2*a(3)+a(4)/10; for j=3:n-2 b(j)=(a(j-2)+a(j-1)+a(j)+a(j+1)+a(j+2)/5; end b(n-1)=(a(n-3)+2*a(n-2)+3*a(n-1)+4*a(n)/10; b(n)=(-a(n-3)+a(n-2)+2*a(n-1)+3*a(n)/5; a=b;end%將a賦值給yy=a;%將分成2行1列的圖形窗口的第1列設(shè)為當(dāng)前繪圖區(qū)域subplot(2,1,1); %繪制平滑前的時(shí)程曲線圖形

8、plot(t,x); %添加橫向坐標(biāo)軸的標(biāo)注xlabel('時(shí)間 (s)'); %添加縱向坐標(biāo)軸的標(biāo)注 ylabel('加速度(g)'); %在圖幅上添加坐標(biāo)網(wǎng)格grid on; %將分成2行1列的圖形窗口的第2列設(shè)為當(dāng)前繪圖區(qū)域subplot(2,1,2);%繪制平滑后的時(shí)程曲線圖形 plot(t,y); %添加橫向坐標(biāo)軸的標(biāo)注xlabel('時(shí)間 (s)'); %添加縱向坐標(biāo)軸的標(biāo)注 ylabel('加速度(g)'); %在圖幅上添加坐標(biāo)網(wǎng)格grid on; %打開文件輸出平滑后的數(shù)據(jù)fid=fopen(fno,'w&

9、#39;); for k=1:n %每行寫兩個(gè)實(shí)型數(shù)據(jù)，t為時(shí)間，y為平滑后的數(shù)據(jù) fprintf(fid,'%f %fn',t(k),y(k); end status=fclose(fid); b、處理結(jié)果二、 振動(dòng)信號(hào)時(shí)域處理振動(dòng)信號(hào)時(shí)域處理是對(duì)振動(dòng)波形的分析。從記錄的時(shí)程信號(hào)中提取各種有用的信息或?qū)⒂涗浀臅r(shí)程信號(hào)轉(zhuǎn)換成所需要的形式。通過不同時(shí)域處理方法，可以確定實(shí)測(cè)波形的最大幅值和時(shí)間歷程，求出相位滯后和波形的時(shí)間滯后，有選擇地濾除或保留實(shí)測(cè)波形的某些頻率成分，消除波形的畸變狀況，再現(xiàn)真實(shí)波形面貌。數(shù)字濾波是通過數(shù)學(xué)手段從所采集的離散信號(hào)中選取人們所感興趣的一部分信號(hào)的處

10、理方法。它的主要作用有濾除測(cè)試信號(hào)中的噪聲或虛假成分、提高信噪比、平滑分析數(shù)據(jù)、抑制干擾信號(hào)、分離頻率分量等。數(shù)字進(jìn)入濾波器后，部分頻率可以通過，部分則受阻擋。1、數(shù)字濾波（1）、頻域低通和帶通濾波a、頻域低通和帶通濾波的程序%clearclcclose all hidden%fni=input('頻域帶通濾波-輸入數(shù)據(jù)文件名:','s'); fid = fopen(fni,'r'); sf = fscanf(fid,'%f',1); %采樣頻率fmin = fscanf(fid,'%f',1); %最小截止頻率f

11、max = fscanf(fid,'%f',1); %最大截止頻率sx = fscanf(fid,'%s',1); %讀入橫向坐標(biāo)軸的標(biāo)注 sy = fscanf(fid,'%s',1); %讀入縱向坐標(biāo)軸的標(biāo)注 fno = fscanf(fid,'%s',1); %輸出數(shù)據(jù)文件名x = fscanf(fid,'%f',1,inf);%輸入數(shù)據(jù)存成行向量status=fclose(fid);%取信號(hào)數(shù)據(jù)長(zhǎng)度 n=length(x); %建立離散時(shí)間列向量 t=(0:1/sf:(n-1)/sf)' %取大于并

12、最接近n的2的冪次方為FFT長(zhǎng)度nfft=2nextpow2(n); %四舍五入取整求最小截止頻率對(duì)應(yīng)數(shù)組元素的下標(biāo)ni=round(fmin*nfft/sf+1); %四舍五入取整求最大截止頻率對(duì)應(yīng)數(shù)組元素的下標(biāo)na=round(fmax*nfft/sf+1); %進(jìn)行FFT變換,結(jié)果存于yy=fft(x,nfft); %建立一個(gè)長(zhǎng)度為nfft元素全為0的向量 a=zeros(1,nfft); %將y的正頻率帶通內(nèi)的元素賦值給a a(ni:na)=y(ni:na); %將y的負(fù)頻率帶通內(nèi)的元素賦值給aa(nfft-na+1:nfft-ni+1)=y(nfft-na+1:nfft-ni+1);

13、 %進(jìn)行FFT逆變換,結(jié)果存于yy=ifft(a,nfft); %取逆變換的實(shí)部n個(gè)元素為濾波結(jié)果列向量y=(real(y(1:n)' %繪制濾波前的時(shí)程曲線圖形subplot(2,1,1); plot(t,x); %添加橫向坐標(biāo)軸的標(biāo)注xlabel(sx); %添加縱向坐標(biāo)軸的標(biāo)注ylabel(sy); grid on; %繪制濾波后的時(shí)程曲線圖形subplot(2,1,2); plot(t,y); %添加橫向坐標(biāo)軸的標(biāo)注xlabel(sx); %添加縱向坐標(biāo)軸的標(biāo)注ylabel(sy); grid on;%打開文件輸出濾波后的數(shù)據(jù) fid=fopen(fno,'w'

14、);for k=1:n fprintf(fid,'%f %fn',t(k),y(k); end status=fclose(fid); b、處理結(jié)果（2）、頻域高通和帶阻濾波a、頻域高通和帶阻濾波都的程序%頻域高通和帶阻濾波%clearclcclose all hidden%fni=input('頻域帶阻濾波-輸入數(shù)據(jù)文件名:','s'); fid=fopen(fni,'r'); sf = fscanf(fid,'%f',1); %采樣頻率fmin = fscanf(fid,'%f',1); %最小

15、截止頻率fmax = fscanf(fid,'%f',1); %最大截止頻率sx = fscanf(fid,'%s',1); %讀入橫向坐標(biāo)軸的標(biāo)注 sy = fscanf(fid,'%s',1); %讀入縱向坐標(biāo)軸的標(biāo)注 fno = fscanf(fid,'%s',1); %輸出數(shù)據(jù)文件名x = fscanf(fid,'%f',1,inf);%輸入數(shù)據(jù)存成行向量status=fclose(fid);%取信號(hào)數(shù)據(jù)長(zhǎng)度 n=length(x); %建立離散時(shí)間列向量 t=(0:1/sf:(n-1)/sf)'

16、%取大于并最接近n的2的冪次方為FFT長(zhǎng)度nfft=2nextpow2(n); %四舍五入取整求最小截止頻率對(duì)應(yīng)數(shù)組元素的下標(biāo)ni=round(fmin*nfft/sf+1); %四舍五入取整求最大截止頻率對(duì)應(yīng)數(shù)組元素的下標(biāo)na=round(fmax*nfft/sf+1); %進(jìn)行FFT變換,結(jié)果存于yy=fft(x,nfft); %建立一個(gè)長(zhǎng)度為nfft元素全為0的向量 a=zeros(1,nfft); %將y的正頻率帶阻內(nèi)的元素賦值為0y(ni:na)=a(ni:na); %將y的負(fù)頻率帶阻內(nèi)的元素賦值為0y(nfft-na+1:nfft-ni+1)=a(nfft-na+1:nfft-ni

17、+1); %進(jìn)行IFFT變換,結(jié)果存于aa=ifft(y,nfft); %取逆變換的實(shí)部n個(gè)元素為濾波結(jié)果列向量y=real(a(1:n);%繪制濾波前的時(shí)程曲線圖形subplot(2,1,1);plot(t,x); %添加橫向坐標(biāo)軸的標(biāo)注xlabel(sx); %添加縱向坐標(biāo)軸的標(biāo)注ylabel(sy); grid on;%繪制濾波后的時(shí)程曲線圖形subplot(2,1,2); plot(t,y);%添加橫向坐標(biāo)軸的標(biāo)注xlabel(sx); %添加縱向坐標(biāo)軸的標(biāo)注ylabel(sy); grid on; %打開文件輸出濾波后的數(shù)據(jù) fid=fopen(fno,'w'); f

18、or k=1:n fprintf(fid,'%f %fn',t(k),y(k); end status=fclose(fid); b、處理結(jié)果2、振動(dòng)信號(hào)的積分和微分變換在振動(dòng)信號(hào)測(cè)試過程中，由于儀器設(shè)備或測(cè)試環(huán)境的限制，有的物理量往往需要通過對(duì)采集到的其他物理量進(jìn)行轉(zhuǎn)換處理才能得到。頻域積分a、頻域積分的程序%頻域積分%clearclcclose all hidden%fni=input('頻域積分-輸入數(shù)據(jù)文件名:','s'); fid=fopen(fni,'r'); sf = fscanf(fid,'%f',

19、1); %采樣頻率fmin = fscanf(fid,'%f',1); %最小截止頻率fmax = fscanf(fid,'%f',1); %最大截止頻率c = fscanf(fid,'%f',1); %單位變換系數(shù)it = fscanf(fid,'%f',1); %積分次數(shù)sx = fscanf(fid,'%s',1); %橫向坐標(biāo)軸的標(biāo)注 sy1 = fscanf(fid,'%s',1); %縱向坐標(biāo)軸輸入單位的標(biāo)注 sy2 = fscanf(fid,'%s',1); %縱向坐標(biāo)

20、軸輸出單位的標(biāo)注 fno = fscanf(fid,'%s',1); %輸出數(shù)據(jù)文件名x = fscanf(fid,'%f',1,inf);%輸入數(shù)據(jù)存成行向量 status=fclose(fid);n=length(x); %建立時(shí)間向量t=0:1/sf:(n-1)/sf; %大于并最接近n的2的冪次方為FFT長(zhǎng)度 nfft=2nextpow2(n); %FFT變換y=fft(x,nfft); %計(jì)算頻率間隔（Hz/s）df=sf/nfft;%計(jì)算指定頻帶對(duì)應(yīng)頻率數(shù)組的下標(biāo)ni=round(fmin/df+1); na=round(fmax/df+1);%計(jì)算

21、圓頻率間隔（rad/s） dw=2*pi*df; %建立正的離散圓頻率向量w1=0:dw:2*pi*(0.5*sf-df); %建立負(fù)的離散圓頻率向量w2=2*pi*(0.5*sf-df):-dw:0; %將正負(fù)圓頻率向量組合成一個(gè)向量w=w1,w2; %以積分次數(shù)為指數(shù)，建立圓頻率變量向量w=w.it; %進(jìn)行積分的頻域變換a=zeros(1,nfft);a(2:nfft-1)=y(2:nfft-1)./w(2:nfft-1); if it=2 %進(jìn)行二次積分的相位變換 y=-a; else %進(jìn)行一次積分的相位變換 real(y)=imag(a); imag(y)=-real(a); en

22、da=zeros(1,nfft);%消除指定正頻帶外的頻率成分 a(ni:na)=y(ni:na); %消除指定負(fù)頻帶外的頻率成分 a(nfft-na+1:nfft-ni+1)=y(nfft-na+1:nfft-ni+1); %IFFT變換 y=ifft(a,nfft); %取逆變換的實(shí)部n個(gè)元素并乘以單位變換系數(shù)為積分結(jié)果y=real(y(1:n)*c; %繪制積分前的時(shí)程曲線圖形 subplot(2,1,1); plot(t,x); %添加橫向坐標(biāo)軸的標(biāo)注xlabel(sx); %添加縱向坐標(biāo)軸的標(biāo)注 ylabel(sy1); grid on; %繪制積分后的時(shí)程曲線圖形 subplot(

23、2,1,2); plot(t,y);%添加橫向坐標(biāo)軸的標(biāo)注xlabel(sx); %添加縱向坐標(biāo)軸的標(biāo)注 ylabel(sy2); grid on; %打開文件輸出積分后的數(shù)據(jù) fid=fopen(fno,'w'); for k=1:n fprintf(fid,'%f %fn',t(k),y(k); endstatus=fclose(fid);（b）、處理結(jié)果三、 振動(dòng)信號(hào)頻域處理我們?cè)诖藢?duì)所測(cè)數(shù)據(jù)進(jìn)行譜分析處理。頻域處理也稱譜分析，是建立在傅利葉變換基礎(chǔ)上的時(shí)頻變換處理，所得的結(jié)果是以頻率為變量的函數(shù)，稱譜函數(shù)。下面進(jìn)行信號(hào)譜分析（a）、信號(hào)譜分析的程序%隨機(jī)

24、信號(hào)譜分析%clearclcclose all hiddenformat long%fni=input('隨機(jī)信號(hào)譜分析-輸入數(shù)據(jù)文件名:','s'); fid=fopen(fni,'r'); %分析內(nèi)容(1=自譜,2=互譜,3=頻響函數(shù),4=相干函數(shù))fun = fscanf(fid,'%d',1); sf = fscanf(fid,'%f',1); %采樣頻率nfft = fscanf(fid,'%d',1);%FFT長(zhǎng)度%窗函數(shù)(1=矩形,2=漢寧,3=海明,4=布萊克曼,5=三角)win =

25、 fscanf(fid,'%d',1); fno = fscanf(fid,'%s',1); %輸出數(shù)據(jù)文件名%按行讀入數(shù)據(jù),第1行激勵(lì),第2行響應(yīng) a = fscanf(fid,'%f',2,inf); status=fclose(fid);%取激勵(lì)信號(hào)向量x=a(1,:); %取相對(duì)的響應(yīng)信號(hào)向量y=a(2,:)-a(1,:);%建立頻率向量f=0:sf/nfft:sf/2-sf/nfft; %加窗處理switch win case 1 %矩形窗 w=boxcar(nfft); case 2 %漢寧窗 w=hanning(nfft); cas

26、e 3 %海明窗 w=hamming(nfft); case 4 %布萊克曼窗 w=blackman(nfft); case 5 %三角窗 w=triang(nfft); otherwise w=boxcar(nfft);endswitch fun case 1 %計(jì)算自譜函數(shù) z=psd(y,nfft,sf,w,nfft/2); case 2 %計(jì)算互譜函數(shù) z=csd(x,y,nfft,sf,w,nfft/2); case 3 %計(jì)算頻響函數(shù) z=tfe(x,y,nfft,sf,w,nfft/2); case 4 %計(jì)算相干函數(shù) z=cohere(x,y,nfft,sf,w,nfft/2)

27、; otherwise ;end%繪制幅頻曲線圖 nn=1:nfft/4;subplot(2,1,1);plot(f(nn),abs(z(nn);xlabel('頻率 (Hz)'); ylabel('幅值'); grid on; if fun>1&fun<4 %繪制相頻曲線圖 subplot(2,1,2); plot(f(nn),angle(z(nn); xlabel('頻率 (Hz)'); ylabel('相位'); grid on; end%以寫的方式打開文件或建立一個(gè)新文件fid=fopen(fno,&#

28、39;w'); %輸出自譜和相干函數(shù)數(shù)據(jù)if fun>1&fun<4 for k=1:nfft/2%每行輸出2個(gè)實(shí)型數(shù)據(jù)，f為頻率，z為幅值 fprintf(fid,'%f %fn',f(k),abs(z(k); end%輸出互譜和頻響函數(shù)數(shù)據(jù)else for k=1:nfft/2%每行輸出3個(gè)實(shí)型數(shù)據(jù)，f為頻率，z的實(shí)、虛部 fprintf(fid,'%f %f %fn',f(k),real(z(k),imag(z(k); endendstatus=fclose(fid);（b）、處理結(jié)果無3、 生成白噪聲信號(hào)（a）、程序%生成白噪

29、聲信號(hào)%clearclcclose all hidden%fni=input('生成白噪聲信號(hào)-輸入數(shù)據(jù)文件名:','s'); fid=fopen(fni,'r'); sf = fscanf(fid,'%f',1); %采樣頻率fi = fscanf(fid,'%f',1); %最小截止頻率fa = fscanf(fid,'%f',1); %最大截止頻率tl = fscanf(fid,'%f',1); %時(shí)間長(zhǎng)度(秒)am = fscanf(fid,'%f',1);

30、%最大幅值fno = fscanf(fid,'%s',1); %輸出數(shù)據(jù)文件名status=fclose(fid); %計(jì)算白噪聲信號(hào)數(shù)據(jù)長(zhǎng)度nn=round(sf*tl)+1; %建立信號(hào)離散時(shí)間向量t t=0:1/sf:(n-1)/sf; %大于并最接近n的2的冪次方為FFT長(zhǎng)度 nfft=2nextpow2(n); %四舍五入取整求最小截止頻率對(duì)應(yīng)數(shù)組元素的下標(biāo)ni=round(fi*nfft/sf+1); %四舍五入取整求最大截止頻率對(duì)應(yīng)數(shù)組元素的下標(biāo)na=round(fa*nfft/sf+1); %建立一個(gè)長(zhǎng)度為nfft/2元素全為0的向量 r=zeros(1,nff

31、t/2); %將r的正頻率帶通內(nèi)的元素賦值1，生成幅值譜 r(ni:na)=1; %生成02pi的隨機(jī)數(shù)為隨機(jī)相位a=2*pi*rand(1,nfft/2); %將幅值譜和相位譜轉(zhuǎn)換成實(shí)部和虛部b=r.*exp(i*a)*nfft/2; %將正負(fù)圓頻率向量組合成一個(gè)向量 a=b,b(nfft/2:-1:1); %IFFT變換，并取實(shí)部為生成白噪聲信號(hào)x=real(ifft(a,nfft); %取指定長(zhǎng)度的數(shù)據(jù)并使數(shù)據(jù)的最大值為指定的幅值y=am*x(1:n)/max(abs(x(1:n); %定義自功率譜的FFT長(zhǎng)度nft=1024; %計(jì)算生成白噪聲信號(hào)自功率譜y1=psd(y,nft,sf

32、,boxcar(nft),nft/2); %定義自功率譜的離散頻率向量 f=0:sf/nft:(nft/2-1)*sf/nft; %繪制白噪聲信號(hào)隨時(shí)間變化的曲線圖subplot(2,1,1); plot(t,y);xlabel('時(shí)間(s)'); ylabel('幅值'); grid on;%繪制自功率譜曲線圖subplot(2,1,2); plot(f(1:nft/4),y1(1:nft/4); xlabel('頻率(Hz)'); ylabel('幅值'); grid on;%打開文件輸出生成的白噪聲信號(hào)數(shù)據(jù) fid=fope

33、n(fno,'w'); for k=1:n %每一行輸出兩個(gè)實(shí)型數(shù)據(jù)，t為時(shí)間，y為白噪聲信號(hào)值 fprintf(fid,'%f %fn',t(k),y(k); end status=fclose(fid); （b）、運(yùn)行結(jié)果四、 試驗(yàn)?zāi)B(tài)參數(shù)識(shí)別的頻域識(shí)別試驗(yàn)?zāi)B(tài)參數(shù)識(shí)別是振動(dòng)信號(hào)處理的一個(gè)重要的組成部分，它的主要任務(wù)是從測(cè)試所得的數(shù)據(jù)中，確定振動(dòng)系統(tǒng)的模態(tài)參數(shù)的估計(jì)，其中包括模態(tài)固有頻率、模態(tài)阻尼比、模態(tài)振形以及模態(tài)質(zhì)量和模態(tài)剛度等。下面進(jìn)行有理分式多項(xiàng)式法模態(tài)參數(shù)識(shí)別（a） 程序%有理分式多項(xiàng)式法模態(tài)參數(shù)識(shí)別%clearclcclose all hidde

34、nformat long%fni=input('有理分式多項(xiàng)式法模態(tài)參數(shù)識(shí)別-輸入數(shù)據(jù)文件名:','s'); fid=fopen(fni,'r'); mn = fscanf(fid,'%d',1); %模態(tài)階數(shù)df = fscanf(fid,'%f',1); %頻率間隔fno = fscanf(fid,'%s',1); %輸出數(shù)據(jù)文件名b = fscanf(fid,'%f',2,inf); %實(shí)測(cè)頻響函數(shù)實(shí)、虛部數(shù)據(jù)status=fclose(fid);%定義冪多項(xiàng)式的階次nm =

35、mn*2;%取頻響函數(shù)的長(zhǎng)度n=length(b(1,:);%建立離散頻率向量f=0:df:(n-1)*df;%建立離散圓頻率向量w=2*pi*f; %建立歸一化離散頻率向量wi = w/max(w);%建立實(shí)測(cè)頻響函數(shù)復(fù)數(shù)向量H=b(1,:)+i*b(2,:); %計(jì)算擬合頻響函數(shù)的分子和分母系數(shù)向量A,B = invfreqs(H,wi,nm,nm,100);%冪多項(xiàng)式方程求根(零點(diǎn))P = roots(B);%計(jì)算模態(tài)頻率向量 F1 = abs(P)*max(w)/(2*pi);%計(jì)算阻尼比向量 D1 = -real(P)./(abs(P);%計(jì)算振型系數(shù)向量for k=1:nm if

36、k=1 p(1:nm-1)=P(2:nm); else p(1:k-1)=P(1:k-1); p(k:nm-1)=P(k+1:nm); end y=poly(p); S1(k)=polyval(A,P(k)/polyval(y,P(k); end%計(jì)算擬合的頻響函數(shù)H1 = freqs(A,B,wi);%繪制頻響函數(shù)實(shí)部擬合曲線圖figure(2)nn=1:n;subplot(2,1,1);plot(f,real(H(nn),':',f,real(H1(nn),'r');xlabel('頻率 (Hz)'); ylabel('實(shí)部'

37、); legend('實(shí)測(cè)','擬合');grid on;%繪制頻響函數(shù)虛部擬合曲線圖subplot(2,1,2);plot(f,b(2,:),':',f,imag(H1(nn),'r');xlabel('頻率 (Hz)'); ylabel('虛部'); legend('實(shí)測(cè)','擬合'); grid on;%將模態(tài)頻率從小到大排序F2,I=sort(F1); %剔除方程解中的非模態(tài)項(xiàng)(非共軛根)和共軛項(xiàng)(重復(fù)項(xiàng))m=0;for k=1:nm-1 if F2(k)=F2(k+1) continue; end m=m+1; l=I(k); F(m)=F1(l); %模態(tài)頻率 D(m)=D1(l); %阻尼比 S(m)=S1(l); %振型系數(shù) end%打開文件輸出模態(tài)參數(shù)數(shù)據(jù)fid=fopen