MatLab实现FFT与功率谱

FFT和功率谱估计

1. 用Fourier变换求取信号的功率谱---周期图法

clf;
Fs=1000;
N=256;Nfft=256;%数据的长度和FFT所用的数据长度
n=0:N-1;t=n/Fs;%采用的时间序列
xn=sin(2*pi*50*t)+2*sin(2*pi*120*t)+randn(1,N);
Pxx=10*log10(abs(fft(xn,Nfft).^2)/N);%Fourier振幅谱平方的平均值，并转化为dB
f=(0:length(Pxx)-1)*Fs/length(Pxx);%给出频率序列
subplot(2,1,1),plot(f,Pxx);%绘制功率谱曲线
xlabel('频率/Hz');ylabel('功率谱/dB');
title('周期图 N=256');grid on;
Fs=1000;
N=1024;Nfft=1024;%数据的长度和FFT所用的数据长度
n=0:N-1;t=n/Fs;%采用的时间序列
xn=sin(2*pi*50*t)+2*sin(2*pi*120*t)+randn(1,N);
Pxx=10*log10(abs(fft(xn,Nfft).^2)/N);%Fourier振幅谱平方的平均值，并转化为dB
f=(0:length(Pxx)-1)*Fs/length(Pxx);%给出频率序列
subplot(2,1,2),plot(f,Pxx);%绘制功率谱曲线
xlabel('频率/Hz');ylabel('功率谱/dB');
title('周期图 N=256');grid on;

1. 用Fourier变换求取信号的功率谱---分段周期图法
   %思想：把信号分为重叠或不重叠的小段，对每小段信号序列进行功率谱估计，然后取平均值作为整个序列的功率谱
   clf;
   Fs=1000;
   N=1024;Nsec=256;%数据的长度和FFT所用的数据长度
   n=0:N-1;t=n/Fs;%采用的时间序列
   randn('state',0);
   xn=sin(2*pi*50*t)+2*sin(2*pi*120*t)+randn(1,N);
   Pxx1=abs(fft(xn(1:256),Nsec).^2)/Nsec; %第一段功率谱
   Pxx2=abs(fft(xn(257:512),Nsec).^2)/Nsec;%第二段功率谱
   Pxx3=abs(fft(xn(513:768),Nsec).^2)/Nsec;%第三段功率谱
   Pxx4=abs(fft(xn(769:1024),Nsec).^2)/Nsec;%第四段功率谱
   Pxx=10*log10(Pxx1+Pxx2+Pxx3+Pxx4/4);%Fourier振幅谱平方的平均值，并转化为dB
   f=(0:length(Pxx)-1)*Fs/length(Pxx);%给出频率序列
   subplot(2,1,1),plot(f(1:Nsec/2),Pxx(1:Nsec/2));%绘制功率谱曲线
   xlabel('频率/Hz');ylabel('功率谱/dB');
   title('平均周期图（无重叠） N=4*256');grid on;
   %运用信号重叠分段估计功率谱
   Pxx1=abs(fft(xn(1:256),Nsec).^2)/Nsec; %第一段功率谱
   Pxx2=abs(fft(xn(129:384),Nsec).^2)/Nsec;%第二段功率谱
   Pxx3=abs(fft(xn(257:512),Nsec).^2)/Nsec;%第三段功率谱
   Pxx4=abs(fft(xn(385:640),Nsec).^2)/Nsec;%第四段功率谱
   Pxx5=abs(fft(xn(513:768),Nsec).^2)/Nsec;%第四段功率谱
   Pxx6=abs(fft(xn(641:896),Nsec).^2)/Nsec;%第四段功率谱
   Pxx7=abs(fft(xn(769:1024),Nsec).^2)/Nsec;%第四段功率谱
   Pxx=10*log10(Pxx1+Pxx2+Pxx3+Pxx4+Pxx5+Pxx6+Pxx7/7);%Fourier振幅谱平方的平均值，并转化为dB
   f=(0:length(Pxx)-1)*Fs/length(Pxx);%给出频率序列
   subplot(2,1,2),plot(f(1:Nsec/2),Pxx(1:Nsec/2));%绘制功率谱曲线
   xlabel('频率/Hz');ylabel('功率谱/dB');
   title('平均周期图（重叠1/2） N=1024');grid on;
   
   
2. 用Fourier变换求取信号的功率谱---welch方法
   %思想：welch法采用信号重叠分段，加窗函数和FFT算法等计算一个信号序列的自功率谱（PSD）和两个信号序列的互功率谱（CSD），采用MATLAB自
   %带的函数psd
   clf;
   Fs=1000;
   N=1024;Nfft=256;n=0:N-1;t=n/Fs;
   window=hanning(256);
   noverlap=128;
   dflag='none';
   randn('state',0);
   xn=sin(2*pi*50*t)+2*sin(2*pi*120*t)+randn(1,N);
   Pxx=psd(xn,Nfft,Fs,window,noverlap,dflag);
   f=(0:Nfft/2)*Fs/Nfft;
   plot(f,10*log10(Pxx));
   xlabel('频率/Hz');ylabel('功率谱/dB');
   title('PSD--Welch方法');grid on;
   
   
3. 功率谱估计----多窗口法（multitaper method ,MTM法）
   %思想：利用多个正交窗口获得各自独立的近似功率谱估计，综合这些得到一个序列的功率谱估计；相对于普通的周期图有更大的自由度；MTM法采用一个参数：时间带
   %宽积NW，这个参数用以定义计算功率谱所用窗的数目为2*NW-1，NW越大，时间域分辨率越高而频率分辨率越低，使得功率谱估计的波动减小；随着NW的增大
   %，每次估计中谱泄露增多，总功率谱估计的偏差增大
   clf;
   Fs=1000;
   N=1024;Nfft=256;n=0:N-1;t=n/Fs;
   randn('state',0);
   xn=sin(2*pi*50*t)+2*sin(2*pi*120*t)+randn(1,N);
   [Pxx1,f]=pmtm(xn,4,Nfft,Fs); %此处有问题
   subplot(2,1,1),plot(f,10*log10(Pxx1));
   xlabel('频率/Hz');ylabel('功率谱/dB');
   title('多窗口法（MTM）NW=4');grid on;
   [Pxx,f]=pmtm(xn,2,Nfft,Fs);
   subplot(2,1,2),plot(f,10*log10(Pxx));
   xlabel('频率/Hz');ylabel('功率谱/dB');
   title('多窗口法（MTM）NW=2');grid on;
   
   
4. 功率谱估计----最大熵法（maxmum entmpy method,MEM法）
   %思想：假定随机序列为平稳高斯过程利用已知的自相关序列rxx(0),rxx(1),rxx(2)...rxx(p)为基础，外推自相关序列rxx(p+1),rxx(p+2)...保证信息熵最大
   clf;
   Fs=1000;
   N=1024;Nfft=256;n=0:N-1;t=n/Fs;
   window=hanning(256);
   randn('state',0);
   xn=sin(2*pi*50*t)+2*sin(2*pi*120*t)+randn(1,N);
   [Pxx1,f]=pmem(xn,14,Nfft,Fs); %此处有问题
   subplot(2,1,1),plot(f,10*log10(Pxx1));
   xlabel('频率/Hz');ylabel('功率谱/dB');
   title('最大熵法（MEM）Order=14');grid on;
   %采用Welch方法估计功率谱
   noverlap=128;
   dflag='none';
   subplot(2,1,2)
   psd(xn,Nfft,Fs,window,noverlap,dflag);
   xlabel('频率/Hz');ylabel('功率谱/dB');
   title('Welch方法估计功率谱');grid on;
   
   
5. 功率谱估计----多信号分类法（multiple signal classification,music法）
   %注：适用于白白噪声中的多正弦波频率估计
   %思想：将数据自相关矩阵看成是由信号自相关矩阵和噪声自相关矩阵两部分组成，求他们的矩阵特征值向量
   clf;
   Fs=1000;
   N=1024;Nfft=256;n=0:N-1;t=n/Fs;
   randn('state',0);
   xn=sin(2*pi*100*t)+2*sin(2*pi*200*t)+randn(1,N);
   pmusic(xn,[7,1.1],Nfft,Fs,32,16);
   xlabel('频率/KHz');ylabel('功率谱/dB');
   title('Welch方法估计功率谱');grid on;