Fourier變換-Gabor變換-Wigner分布-小波變換實例分析

1、1、分別用短時Fourier，Gabor變換分析下列信號，要求提供程序，圖形結(jié)果并對它們的結(jié)果進行對比分析。采樣頻率FS=1920HZ，采樣長度N=512.Matlab程序如下：fs=1920;%采樣頻率N=512; %采樣長度 t=0:1/fs:(N-1)/fs; %時間序列x1=(1+0.2*sin(2*pi*15*t).*(cos(2*pi*30*t)+0.5*sin(2*pi*15*t)+sin(2*pi*120*t);%信號figure(1)plot(t,x1);%畫想（t）的圖像 y1=fft(x1,N); %對信號進行快速Fourier變換mag1=abs(y1);%求變換后的幅

2、值 k=0:N-1;f1=k*fs/N; figure(2)grid on stem(f1,mag1);%繪制N點DFI的幅頻特性圖xlabel(f1);ylabel(幅值);axis(0,256,0,2*max(abs(y1);%x,y的范圍grid on figure(3)h=window(321,hamming);sig=x1;tfrstft(sig,1:512,512,h);%短時Fourier變換xlabel(時間(秒);ylabel(頻率(Hz); figure(4)q=16;h=window(211,gauss);h=h/norm(h);tfrgabor(x1,128,q,h);

3、%Gabor變換xlabel(時間(秒); ylabel(頻率(Hz);1.1信號的圖形圖1-1 信號時域波形圖1.2信號N點的幅頻特性圖圖1-2 信號的幅頻特性圖對信號進行分析，信號共有個頻率分別是HZ,15HZ,30HZ,45HZ,120HZ，用火柴棍狀表示出來。1.3短時Fourier變換圖圖1-3 短時傅里葉變換圖1.4 Gabor變換圖圖1-4 Gabor變換圖通過上面兩圖可以看出就顯示兩個頻率，分別是30HZ和120HZ，15HZ比較模糊，而0HZ 和45HZ的信號淹沒了，經(jīng)過分析原因可能一是信號的強度不一樣，顯示的清晰度也不一樣；二是采樣頻率過大顯示的比較擁擠。有以上兩圖知：短時

4、Fourier變換和Gabor變換均能顯示在特定時刻該信號的頻率，與Fourier變換相比具有定位的功能。但短時Fourier窗函數(shù)寬度的選擇對時間和頻率分辨率的影響比較大，不能使時間分辨率和頻率分辨率都能提高。2、分別用wigner-ville分布，偽wigner-ville分布，平滑偽wigner-ville分布和Cohen分布分析下列信號：其中=0.25， =5s，。要求提供圖形結(jié)果，并對它們的結(jié)果進行對比分析。Matlab程序如下：a=0.25;t0=5;fs=10;w0=1.57;n=128;t=0:1/fs:(n-1)/fs;x=exp(-a*(t+t0).2+w0*t*j)+ex

5、p(-a*(t-t0).2+w0*t*j);figure(1)plot(abs(x);x=x.;set(gca,xlim,0,n)set(gca,xtick,0:n/4:n) figure(2);tfrwv(hilbert(x);title(Wigner-ville分布)axis(xy);xlabel(時間(秒); ylabel(頻率(Hz); figure(3);tfrpwv(hilbert(x);title(偽Wigner-ville分布)axis(xy);xlabel(時間(秒); ylabel(頻率(Hz); figure(4);tfrspwv(hilbert(x);title(平滑W

6、igner-ville分布)axis(xy);xlabel(時間(秒); ylabel(頻率(Hz);figure(5);tfrcw(hilbert(x);title(cohen時頻分布)axis(xy);xlabel(時間(秒); ylabel(頻率(Hz);2.1信號的時域波形圖圖2-1 信號的時域波形圖2.2 Wiger-ville分布圖2-2 信號的Wiger-Ville分布圖2.3 偽Wiger-ville分布圖2-3 信號的偽Wiger-Ville分布圖2.4 平滑偽Wiger-ville分布圖2-4信號的平滑偽Wiger-Ville分布圖2.5 Cohen時頻分布圖2-5 信號的

7、Cohen時頻分布WVD分布有明顯的缺點，就是有交叉項的存在，而以后的偽WVD分布，平滑WVD分布，平滑偽WVD分布，Cohen類分布都是對WVD分布的改進，通過加窗來抑制交叉項，使圖像變得更平滑。3 對工程信號進行插值和抽取分析所用信號是旋轉(zhuǎn)機械中軸承故障提取信號，所有的信號有個數(shù)據(jù)，為了方便分析和運行，截取了其中的2048個數(shù)據(jù),進行插值與抽取。程序如下：s=xlsread(G:數(shù)字信號處理應(yīng)用現(xiàn)代信號處理課件及程序程序軸承故障.xls);S=fft(s);N=length(s);n=0:N-1;fs=1000;f=n*fs/N;figure(1)subplot(4,1,1);plot(a

8、bs(S);title(原信號頻譜幅值);s4=interp(s,4); %對信號進行4倍插值S4=fft(s4); subplot(4,1,2);plot(abs(S4);title(四倍插值后的信號頻譜幅值);s8=interp(s,8); %對信號進行8倍插值S8=fft(s8); subplot(4,1,3);plot(abs(S8);title(8倍插值后的信號頻譜幅值);s16=interp(s,16);%對信號進行16倍插值S16=fft(s4);subplot(4,1,4);plot(abs(S16);title(十六倍插值后的信號頻譜幅值);figure(2) subplot

9、(4,1,1);plot(abs(S);title(原信號頻譜幅值);b4=decimate(s,4);%對信號進行四倍抽取B4=fft(b4); subplot(4,1,2);plot(abs(B4);title(四倍抽取后的信號頻譜幅值);b8=decimate(s,8);%對信號進行八倍抽取B8=fft(b8); subplot(4,1,3);plot(abs(B8);title(八倍抽取后的信號頻譜幅值);b16=decimate(s,16);%對信號進行16倍抽取B16=fft(b16);subplot(4,1,4);plot(abs(B16);title(十六倍抽取后的信號頻譜幅值

10、)3.1 插值后的圖像圖3-1 信號插值后的圖像3.2抽取后的圖像圖3-2 信號抽取后的圖像4 對工程信號進行各種時頻分析后的程序和圖形結(jié)果分析所用信號是旋轉(zhuǎn)機械中軸承故障提取信號，所有的信號有個數(shù)據(jù)，為了方便分析和運行，截取了其中的2048個數(shù)據(jù)，然后分別對信號做各種時頻分析，如短時傅里葉變換、gabor變換、cohen變換和小波變換等，給出程序和圖形結(jié)果。1、短時傅里葉變換：s=xlsread(G:數(shù)字信號處理應(yīng)用現(xiàn)代信號處理課件及程序程序軸承故障.xls);x1=s;fs=1000;N=length(x1);t=1:N;figure(1);plot(t,x1,LineWidth,2);x

11、label(時間t/s);ylabel(幅值A(chǔ));y1=fft(x1,N); %對信號進行快速Fourier變換mag1=abs(y1);%求得Fourier變換后的振幅k=0:N-1;f1=k*fs/N; figure(2)grid on %網(wǎng)格開啟stem(f1,mag1);%繪制N點DFI的幅頻特性圖xlabel(f1);ylabel(幅度);axis(0,512,0,1.2*max(abs(y1);grid on figure(2)h=window(1111,hamming);sig=x1;tfrstft(sig,1:1024,1024,h);xlabel(時間(秒);ylabel(頻

12、率(Hz);1-1 信號的波形圖圖1-1 工程信號的波形圖1-2 信號的幅頻特性圖圖1-2 信號的幅頻特性圖1-3 信號的短時傅里葉變換圖圖1-3 信號的短時傅里葉變換圖2、gabor變換：s=xlsread( G:數(shù)字信號處理應(yīng)用現(xiàn)代信號處理課件及程序程序軸承故障.xls);x1=s;fs=1000;N=length(x1);t=1:N;figure(1);plot(t,x1,LineWidth,2);xlabel(時間t/s);ylabel(振幅A);figure(3)q=16;h=window(1911,gauss);h=h/norm(h);tfrgabor(x1,128,q,h);xl

13、abel(時間(秒);ylabel(頻率(Hz);2-1 對信號進行g(shù)abor變換圖2-1 信號的gabor變換圖對信號進行g(shù)abor變換的過程中，在過抽樣的情況下，對信號進行時頻分析。gabor變換是短時傅里葉變換加窗后的一種特殊情況。3、Cohen類時頻分布s=xlsread( G:數(shù)字信號處理應(yīng)用現(xiàn)代信號處理課件及程序程序軸承故障.xls);x1=s;fs=1000;N=length(x1);t=1:N;figure(1);plot(t,x1,LineWidth,2);xlabel(時間t/s);ylabel(振幅A);figure(2);tfrwv(hilbert(x);title(W

14、igner-ville分布)axis(xy);xlabel(時間(秒); ylabel(頻率(Hz); figure(3);tfrpwv(hilbert(x);title(偽Wigner-ville分布)axis(xy);xlabel(時間(秒); ylabel(頻率(Hz); figure(4);tfrspwv(hilbert(x);title(平滑偽Wigner-ville分布)axis(xy);xlabel(時間(秒); ylabel(頻率(Hz); figure(5);tfrcw(hilbert(x);title(cohen時頻分布)axis(xy);xlabel(時間(秒); yla

15、bel(頻率(Hz);3-1 Wiger-ville分布圖3-1 信號的Wiger-Ville分布圖3-2 偽Wiger-ville分布圖3-2 信號的偽Wiger-Ville分布圖3-3 平滑偽Wiger-ville分布圖3-3 信號的平滑偽Wiger-Ville分布圖3-4 Cohen時頻分布圖3-4 Cohen時頻分布4、小波變換s=xlsread( G:數(shù)字信號處理應(yīng)用現(xiàn)代信號處理課件及程序程序軸承故障.xls);x1=s;fs=1000;N=length(x1);t=1:N;figure(1);plot(t,x1,LineWidth,2);xlabel(時間t/s);ylabel(振

16、幅A);c,l = wavedec(x,6,db3); %重構(gòu)第1-6層逼近系數(shù)a6 = wrcoef(a,c,l,db3,6);a5 = wrcoef(a,c,l,db3,5);a4 = wrcoef(a,c,l,db3,4);a3 = wrcoef(a,c,l,db3,3);a2 = wrcoef(a,c,l,db3,2);a1 = wrcoef(a,c,l,db3,1);%顯示逼近系數(shù)figure(2)subplot(6,1,1);plot(a6,LineWidth,2);ylabel(a6);subplot(6,1,2);plot(a5,LineWidth,2);ylabel(a5);

17、subplot(6,1,3);plot(a4,LineWidth,2);ylabel(a4);subplot(6,1,4);plot(a3,LineWidth,2);ylabel(a3);subplot(6,1,5);plot(a2,LineWidth,2);ylabel(a2);subplot(6,1,6);plot(a1,LineWidth,2);ylabel(a1); %重構(gòu)第1-6層細(xì)節(jié)系數(shù)d6= wrcoef(d,c,l,db3,6);d5 = wrcoef(d,c,l,db3,5);d4 = wrcoef(d,c,l,db3,4);d3 = wrcoef(d,c,l,db3,3);d2 = wrcoef(d,c,l,db3,2);d1 = wrcoef(d,c,l,db3,1); %顯示細(xì)節(jié)系數(shù)figure(3)subplot(6,1,1);plot(d6,LineWidth,2);ylabel(d6);axis(0 N -5 5);subplot(6,1,2);plot(d5,LineWidth,2);ylabel(d5);axis(0 N -5 5);subplot(6,1,3);plot(d4,LineWidth,2);