講這個(gè)話題，就要先搞清楚頻譜、功率譜的概念，可參考我的另一篇文章

[信號的頻譜 頻譜密度 功率譜密度 能量譜密度]

做信號處理的朋友應(yīng)該都會fft比較熟悉，就是求傅里葉變換。我在這里也不再去講這個(gè)函數(shù)了，但需要注意的一點(diǎn)：實(shí)信號的頻譜關(guān)于0頻對稱，是偶函數(shù)，如果st = cos(2pif0*t)+1; t的長度為4000，那么0頻的位置在第一個(gè)點(diǎn)，做fftshift后，0頻的位置在低2001個(gè)點(diǎn)的位置，fft后的信號關(guān)于第2001個(gè)點(diǎn)對稱，而不是4000個(gè)點(diǎn)左右對稱。

pwelch是用來求功率譜的，采用Welch平均周期法對信號進(jìn)行譜估計(jì)，它通過分段選取數(shù)據(jù)進(jìn)行加窗求功率，再進(jìn)行平均，pwelch函數(shù)的使用方式為：

pxx = pwelch(x,window,noverlap,nfft)
[pxx,f] = pwelch(x,window,noverlap,f,fs)

其中，

• X表示輸入序列；
• window：當(dāng)window是一個(gè)數(shù)值時(shí)，表示窗函數(shù)長度，即分段長度L，默認(rèn)的窗函數(shù)為hamming窗；當(dāng)window是一個(gè)序列時(shí)，表示窗函數(shù)序列；
• NFFT表示FFT的點(diǎn)數(shù)，X為實(shí)數(shù)時(shí)，當(dāng)NFFT是偶數(shù)時(shí)，Pxx的長度是(NFFT/2+1)；當(dāng)NFFT是奇數(shù)時(shí)，Pxx的長度是(NFFT+1)/2；X為復(fù)數(shù)時(shí)，Pxx的長度就是NFFT，如果NFFT沒有指定，則默認(rèn)是256或者比X長度大的2的N次冪
• Fs 繪制功率譜曲線的采樣頻率，默認(rèn)值為1
• Pxx表示功率譜估計(jì)值
• F表示Pxx值所對應(yīng)的頻率點(diǎn)
• NOVERLAP指定分段重疊的樣本數(shù) ，如果NOVERLAP=L/2，則可得到重疊50%的Welch法平均周期圖

下面我們分別用fft和fwelch來求信號的功率譜。

clc;close all;clear all;

fs = 10e6;
N = 4000;
t = (0:N-1)/fs;
f0 = 1e5;

st = cos(2*pi*f0*t) + 1;
st_fft = fft(st);
psdx = abs(st_fft(1:end/2+1)).^2/fs/N;  %功率譜密度為能量譜密度除以時(shí)間，摸值的平方即為能量譜
psdx(2:end) = 2*psdx(2:end);            %乘2是因?yàn)閒ft結(jié)果是對稱的，在計(jì)算功率時(shí)需要把功率加回來；第一個(gè)點(diǎn)是0頻，這個(gè)點(diǎn)并不對稱
freq = linspace(0,fs/2,length(psdx));
[pxx,f] = pwelch(st,rectwin(N),32,N,fs);
figure;plot(freq,psdx);title('fft方法求功率譜密度');grid on
figure;plot(f,pxx);title('fwelch方法求功率譜密度');grid on