matlab中周期图功率谱法的实现原理

matlab中的周期图功率谱法原理是通过计算采样信号的FFT，获得离散点的幅度，再根据幅度与功率之间的关系，转换为离散点的功率，再通过坐标变换将离散点的功率图转换为连续功率谱密度。

1、计算采样信号x(n)的DFT，使用FFT方法来计算。如果此处将复频率处的幅度对称到物理实际频率，得到的就是单边谱，否则就是双边谱。



2、根据正余弦信号功率与幅度的关系以及直流功率与幅度的关系，将幅度转换为离散功率谱。



3、对横纵坐标进行转换，横坐标乘以频率分辨率转换为实际连续物理频率，纵坐标除以频率分辨率转换为功率谱密度。

将上述步骤综合后，得：



等式最右边的2倍是双边到单边的转换，乘以Ts是为了满足单位转换，计算出的PSD单位为V^2/Hz。

下面是自己编写的周期图法功率谱和调用matlab中的周期图法函数得到的功率谱的比较。

%本程序旨在分析单边功率谱密度的计算方法，对比matlab中自带的periodogram方法
clc;
clear;
close all;
%% 定义计算参数
F0=16;%信号频率
N_sample=2048;%总采样点数
Fs=512;%采样频率
f_index=(0:N_sample/2)/N_sample*Fs;%计算FFT之后的频率标注
n=0:N_sample-1;%采样序号
t=n*1/Fs;%时间标度
A=1;%信号幅度

%% 具体计算
x = A*sin(2*pi*F0*t)+2*A*cos(2*pi*4*F0*t)+3*A*cos(2*pi*3*F0*t)+1;  %  +randn(size(t))
[Pxx,f]=periodogram(x,[],'onesided',N_sample,Fs);%周期图法计算功率谱

x_FFT=abs(fft(x,N_sample))./N_sample*2;%计算FFT
x_FFT(1)=x_FFT(1)/2;%直流分量的对应
x_FFT(N_sample/2+1)=x_FFT(N_sample/2+1)/2;%奈奎斯特频率点处的对应

FFT_PSD=x_FFT.^2/2;%利用FFT计算功率谱密度
FFT_PSD(1)=x_FFT(1).^2;%直流功率
FFT_PSD(N_sample/2+1)=x_FFT(N_sample/2+1).^2;%奈奎斯特频率处的功率

FFT_PSD=FFT_PSD/Fs*N_sample;

delta=FFT_PSD(1:N_sample/2+1)-Pxx(1:N_sample/2+1)';%求两者的计算差异

%% 开始绘图
figure(1);
subplot(4,1,1)
plot(t,x);
title('时域信号x','FontSize',18);
grid on;
subplot(4,1,2)
plot(f,Pxx);
title('利用periodogram求得的功率谱密度','FontSize',18);
grid on;
subplot(4,1,3)
plot(f_index,x_FFT(1:N_sample/2+1));
title('DFT','FontSize',18);
grid on;
subplot(4,1,4)
plot(f_index,FFT_PSD(1:N_sample/2+1));
title('FFT手动功率谱密度','FontSize',18);
grid on;
figure(2);
plot(f_index,delta);
title('两谱密度的差值','FontSize',18);