最近看了一點(diǎn)FFT和功率譜估計(jì)的書籍，有些方法隨手記下，剛興趣的朋友可以看下。
%用Fourier變換求取信號(hào)的功率譜---周期圖法
clf;
Fs=1000;
N=256;Nfft=256;%數(shù)據(jù)的長(zhǎng)度和FFT所用的數(shù)據(jù)長(zhǎng)度
n=0:N-1;t=n/Fs;%采用的時(shí)間序列
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振幅譜平方的平均值，并轉(zhuǎn)化為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;%數(shù)據(jù)的長(zhǎng)度和FFT所用的數(shù)據(jù)長(zhǎng)度
n=0:N-1;t=n/Fs;%采用的時(shí)間序列
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振幅譜平方的平均值，并轉(zhuǎn)化為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;

%用Fourier變換求取信號(hào)的功率譜---分段周期圖法
%思想：把信號(hào)分為重疊或不重疊的小段，對(duì)每小段信號(hào)序列進(jìn)行功率譜估計(jì)，然后取平均值作為整個(gè)序列的功率譜
clf;
Fs=1000;
N=1024;Nsec=256;%數(shù)據(jù)的長(zhǎng)度和FFT所用的數(shù)據(jù)長(zhǎng)度
n=0:N-1;t=n/Fs;%采用的時(shí)間序列
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振幅譜平方的平均值，并轉(zhuǎn)化為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;

%運(yùn)用信號(hào)重疊分段估計(jì)功率譜
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振幅譜平方的平均值，并轉(zhuǎn)化為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;

%用Fourier變換求取信號(hào)的功率譜---welch方法
%思想：welch法采用信號(hào)重疊分段，加窗函數(shù)和FFT算法等計(jì)算一個(gè)信號(hào)序列的自功率譜（PSD）和兩個(gè)信號(hào)序列的互功率譜（CSD），采用MATLAB自
%帶的函數(shù)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;

%功率譜估計(jì)----多窗口法（multitaper method ,MTM法）
%思想：利用多個(gè)正交窗口獲得各自獨(dú)立的近似功率譜估計(jì)，綜合這些得到一個(gè)序列的功率譜估計(jì)；相對(duì)于普通的周期圖有更大的自由度；MTM法采用一個(gè)參數(shù)：時(shí)間帶
%寬積NW，這個(gè)參數(shù)用以定義計(jì)算功率譜所用窗的數(shù)目為2*NW-1，NW越大，時(shí)間域分辨率越高而頻率分辨率越低，使得功率譜估計(jì)的波動(dòng)減小；隨著NW的增大
%，每次估計(jì)中譜泄露增多，總功率譜估計(jì)的偏差增大
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;

%功率譜估計(jì)----最大熵法（maxmum entmpy method,MEM法）
%思想：假定隨機(jī)序列為平穩(wěn)高斯過程利用已知的自相關(guān)序列rxx(0),rxx(1),rxx(2)...rxx(p)為基礎(chǔ)，外推自相關(guān)序列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方法估計(jì)功率譜
noverlap=128;
dflag='none';
subplot(2,1,2)
psd(xn,Nfft,Fs,window,noverlap,dflag);
xlabel('頻率/Hz');ylabel('功率譜/dB');
title('Welch方法估計(jì)功率譜');grid on;

%功率譜估計(jì)----多信號(hào)分類法（multiple signal classification,music法）
%注：適用于白白噪聲中的多正弦波頻率估計(jì)
%思想：將數(shù)據(jù)自相關(guān)矩陣看成是由信號(hào)自相關(guān)矩陣和噪聲自相關(guān)矩陣兩部分組成，求他們的矩陣特征值向量
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方法估計(jì)功率譜');grid on;