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

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

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

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

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

其中，

• X表示輸入序列；
• window：當window是一個數值時，表示窗函數長度，即分段長度L，默認的窗函數為hamming窗；當window是一個序列時，表示窗函數序列；
• NFFT表示FFT的點數，X為實數時，當NFFT是偶數時，Pxx的長度是(NFFT/2+1)；當NFFT是奇數時，Pxx的長度是(NFFT+1)/2；X為復數時，Pxx的長度就是NFFT，如果NFFT沒有指定，則默認是256或者比X長度大的2的N次冪
• Fs 繪制功率譜曲線的采樣頻率，默認值為1
• Pxx表示功率譜估計值
• F表示Pxx值所對應的頻率點
• NOVERLAP指定分段重疊的樣本數 ，如果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;  %功率譜密度為能量譜密度除以時間，摸值的平方即為能量譜
psdx(2:end) = 2*psdx(2:end);            %乘2是因為fft結果是對稱的，在計算功率時需要把功率加回來；第一個點是0頻，這個點并不對稱
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