Matlab2019b中常用的音频信号分析，快速傅里叶（FFT），功率谱密度（PSD），以及通过FFT求取功率谱密度的问题

本应该为基础知识的维纳辛钦定理的证明，在matlab调试中遇到的大坑。为了验证，本人编写了一个从求取FFT，到多种方法求取PSD的matlab程序。

问题主要集中在x轴和y轴应该如何取值的问题上，单纯计算出数值很容易，但是想把图画成一样的却很麻烦，网上基本搜不到相关的程序，都只是定理。

这个程序调试了三天才搞定，如果对你有帮助，请关注我，顺手点个赞，感谢。

%% 这个test主要是验证信号的FFT，PSD，以及使用FFT求取PSD的正确性
clear;
clc;
Fs = 1000;            % Sampling frequency
T = 1/Fs;             % Sampling period
L = 1500;             % Length of signal
t = (0:L-1)*T;        % Time vector
S = 0.7*sin(2*pi*50*t) ;%幅值为0.7频率为50的正弦信号


%% 信号的傅里叶FFT变换
Y = fft(S);
P2 = abs(Y/L);
P1 = P2(1:L/2+1);% 实信号的功率谱是对称的,只用取一半即可
P1(2:end-1) = 2*P1(2:end-1);
f = Fs*(0:(L/2))/L;

%% FFT获得功率谱密度PSD
Y = Y(1:L/2+1);% 实信号的功率谱是对称的,只用取一半即可
psdx = (1/(Fs*L)) * abs(Y).^2;% 注意使用采样频率和fft点数归一化
psdx(2:end-1) = 2*psdx(2:end-1);% 为了使得总功率不变,要乘以两倍
freq = 0:Fs/L:Fs/2;

%% Periodogram法测量谱密度PSD
[pxx,fp] = periodogram(S,rectwin(L),L,Fs);

%% pwelch法测量功率谱密度（PSD）
% noverlap=300;%重叠长度，
% Lc=500;%窗函数长度
% nfft=500;%nfft点数
% [pxx,fp]=pwelch(S,...
%     hamming(Lc),...
%     noverlap,...
%     nfft,...
%     Fs,...
%     'onesided',...
%     'power');
%% 自相关系数求功率谱密度PSD
[r,lags]=xcorr(S,500,'coeff');%信号S的自相关
Pr=fft(r);%自相关的傅里叶变换
L1=size(Pr,2);
Lr=size(r,2);%自相关变换后的数据长度
L3=abs(Pr/Lr);
L4=L3(1:Lr/2+1);% 实信号的功率谱是对称的,只用取一半即可
L4(2:end-1)=2*L4(2:end-1);
f1=Fs*(0:(Lr/2))/Lr;%套用单边变换FFT横坐标变换

%% 绘图部分
subplot(4,1,1)
plot(f1,L4)
title('自相关函数傅里叶正变换')
grid on
xlabel('Frequency (Hz)')
ylabel('Amplitude')
%stem(lags,r)

subplot(4,1,2)
plot(freq,pow2db(psdx))
grid on
title('使用FFT法测量功率谱密度')
xlabel('Frequency (Hz)')
ylabel('Power/Frequency (dB/Hz)')

subplot(4,1,3)
plot(f,P1)
title('原始信号直接进行傅里叶分析')
grid on
xlabel('Frequency (Hz)')
ylabel('Amplitude')

subplot(4,1,4)
plot(fp,pow2db(pxx));
title('Periodogram法测量谱密度PSD')
grid on
xlabel('Frequency (Hz)')
ylabel('Power/Frequency (dB/Hz)')

Copyright © 2020 by RichardYang. All rights reserved.
仅供参考，严禁转载，感谢。