Fourier变换-Gabor变换-Wigner分布-小波变换实例分析

1、分别用短时Fourier，Gabor变换分析下列信号，要求提供程序，图形结果并对它们的结果进行对比分析。采样频率FS=1920HZ，采样长度N=512.()(10.2sin(215))cos(2300.5sin(215))sin(2120)xttttt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);%求变换后的幅值k=0:N-1;f1=k*fs/N;figure(2)gridonstem(f1,mag1);%绘制N点DFI的幅频特性图xlabel('f1');ylabel('幅值’);axis([0,256,0,2*max(abs(y1))]);%x,y的范围gridonfigure(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);%Gabor变换xlabel('时间(秒)');ylabel('频率(Hz)');1.1信号的图形图1-1信号时域波形图1.2信号N点的ＤＦＩ幅频特性图图1-2信号的幅频特性图对信号进行分析，信号共有５个频率分别是０HZ,15HZ,30HZ,45HZ,120HZ，用火柴棍状表示出来。1.3短时Fourier变换图图1-3短时傅里叶变换图1.4Gabor变换图图1-4Gabor变换图通过上面两图可以看出就显示两个频率，分别是30HZ和120HZ，15HZ比较模糊，而0HZ和45HZ的信号淹没了，经过分析原因可能一是信号的强度不一样，显示的清晰度也不一样；二是采样频率过大显示的比较拥挤。有以上两图知：短时Fourier变换和Gabor变换均能显示在特定时刻该信号的频率，与Fourier变换相比具有定位的功能。但短时Fourier窗函数宽度的选择对时间和频率分辨率的影响比较大，不能使时间分辨率和频率分辨率都能提高。2、分别用wigner-ville分布，伪wigner-ville分布，平滑伪wigner-ville分布和Cohen分布分析下列信号：220000()()1ttjtttjtxee其中=0.25，0t=5s，01.57/rads。要求提供图形结果，并对它们的结果进行对比分析。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)+exp(-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('平滑Wigner-ville分布')axis('xy');xlabel('时间(秒)');ylabel('频率(Hz)');figure(5);tfrcw(hilbert(x));title('cohen时频分布')axis('xy');xlabel('时间(秒)');ylabel('频率(Hz)');2.1信号的时域波形图图2-1信号的时域波形图2.2Wiger-ville分布图2-2信号的Wiger-Ville分布图2.3伪Wiger-ville分布图2-3信号的伪Wiger-Ville分布图2.4平滑伪Wiger-ville分布图2-4信号的平滑伪Wiger-Ville分布图2.5Cohen时频分布图2-5信号的Cohen时频分布WVD分布有明显的缺点，就是有交叉项的存在，而以后的伪WVD分布，平滑WVD分布，平滑伪WVD分布，Cohen类分布都是对WVD分布的改进，通过加窗来抑制交叉项，使图像变得更平滑。3对工程信号进行插值和抽取分析所用信号是旋转机械中轴承故障提取信号，所有的信号有120000个数据，为了方便分析和运行，截取了其中的2048个数据,进行插值与抽取。程序如下：s=xlsread('G:\数字信号处理应用\现代信号处理课件及程序\程序\轴承故障.xls');S=fft(s);N=length(s);n=0:N-1;fs=1000;f=n'*fs/N;figure(1)subplot(4,1,1);plot(abs(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(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('十六倍抽取后的信号频谱幅值')3.1插值后的图像图3-1信号插值后的图像3.2抽取后的图像图3-2信号抽取后的图像4对工程信号进行各种时频分析后的程序和图形结果分析所用信号是旋转机械中轴承故障提取信号，所有的信号有120000个数据，为了方便分析和运行，截取了其中的2048个数据，然后分别对信号做各种时频分析，如短时傅里叶变换、gabor变换、cohen变换和小波变换等，给出程序和图形结果。1、短时傅里叶变换：s=xlsread('G:\数字信号处理应用\现代信号处理课件及程序\程序\轴承故障.xls');x1=s;fs=1000;N=length(x1);t=1:N;figure(1);plot(t,x1,'LineWidth',2);xlabel('时间t/s');ylabel('幅值A');y1=fft(x1,N);%对信号进行快速Fourier变换mag1=abs(y1);%求得Fourier变换后的振幅k=0:N-1;f1=k*fs/N;figure(2)gridon%网格开启stem(f1,mag1);%绘制N点DFI的幅频特性图xlabel('f1');ylabel('幅度');axis([0,512,0,1.2*max(abs(y1))]);gridonfigure(2)h=window(1111,'hamming');sig=x1;tfrstft(sig,1:1024,1024,h);xlabel('时间(秒)');ylabel('频率(Hz)');1-1信号的波形图图1-1工程信号的波形图1-2信号的幅频特性图图1-2信号的幅频特性图1-3信号的短时傅里叶变换图图1-3信号的短时傅里叶变换图2、gabor变换：s=xlsread('G:\数字信号处理应用\现代信号处理课件及程序\程序\轴承故障.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);xlabel('时间(秒)');ylabel('频率(Hz)');2-1对信号进行gabor变换图2-1信号的gabor变换图对信号进行gabor变换的过程中，在过抽样的情况下，对信号进行时频分析。gabor变换是短时傅里叶变换加窗后的一种特殊情况。3、Cohen类时频分布s=xlsread('G:\数字信号处理应用\现代信号处理课件及程序\程序\轴承故障.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('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('平滑伪Wigner-ville·分布')axis('xy');xlabel('时间(秒)');ylabel('频率(Hz)');figure(5);tfrcw(hilbert(x));title('cohen时频分布')axis('xy');xlabel('时间(秒)');ylabel('频率(Hz)');3-1Wiger-ville分布图3-1信号的Wiger-Ville分布图3-2伪Wiger-ville分布图3-2信号的伪Wiger-Ville分布图3-3平滑伪Wiger-ville分布图3-3信号的平滑伪Wiger-Ville分布图3-4Cohen时频分布图3-4Cohen时频分布4、小波变换s=xlsread('G:\数字信号处理应用\现代信号处理课件及程序\程序\轴承故障.xls');x1=s;fs=1000;N=length(x1);t=1:N;figure(1);plot(t,x1,'LineWidth',2);xlabel('时间t/s');ylabel('振幅A');[c,l]=wavedec(x,6,'db3');%重构第1-6层逼近系数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);%显示逼近系数figure(2)subp