功率谱MATLAB程序

%直接法又称为周期图法，就是把x（n）的N个数据视为已能量有限的序列，直接计算
%x（n）的离散傅里叶变换，得X（k），然后再取其幅值的平方，作为功率谱的估计
t=linspace（011000）;
signal=4*sin（2*pi*50*t）+5*sin（2*pi*200*t）;
noise=randn（size（t））;
symbol=signal+noise;
Y=fft（symbol128）;
f=1000*（0:127）/128;
P1=Y.*conj（Y）/128; %直接法
plot（f10*log10（P1（1:128）））;
xlabel（‘frequency‘）
ylabel（‘power‘）
title（‘直接法‘）
%或者用下面代码
Fs=1000;% 采样频率
n=0:1/Fs:1;
xn=cos（2*pi*40*n）+3*cos（2*pi*100*n）+randn（size（n））;
window=boxcar（length（xn））;%矩形窗
nfft=1024;
[Pxxf]=periodogram（xnwindownfftFs）;
plot（f10*log10（Pxx））;

clear;
fs=500;  %采样率
df=1;    %频率分辨率
N=floor（fs/df）+1;%计算的序列点数
t=0:1/fs:（N-1）/fs;%截取信号的时间段
f=0:df:fs;%功率谱估计的频率分辨率和范围
xt=sin（2*pi*50*5）+2*sin（2*pi*130*t）+randn（1length（t））;
%截取时间段上的离散信号样本序列



%利用周期图法进行功率谱估计，但是其得出的功率谱很不光滑，相应的估计协方差比较大
%增加采样点数也不能使周期图变得更加平滑，这是周期图的缺点，在后面对其改进。
Px=abs（fft（xt））.^2/（N^2）;%功率谱估计
Pav_tm=sum（xt.^2）/N;%在时域计算信号功率
Pav_fn=sum（Px）;%通过功率谱计算信号功率
figure（1）
subplot（221）   %作出功率谱密度图
plot（f10*log10（Px））;
xlable（‘频率（Hz）‘）;
Ylable（‘功率谱（dB）‘）;
title（‘周期图得出的功率谱估计‘）;



%对周期图进行改进的思想是将信号段进行估计，然后再将这些估计结果进行品均，
%从而减小估计的协方差，使估计功率谱图变得平滑。
%本程序是将以上501点的信号分3段，分别作周期图估计，然后平均。
Px=（abs（fft（xt（1:167）））.^2+abs（fft（xt（168:334）））.^2+abs（fft（xt（335:501）））.^2）/3/（（N/3）^2）;%分为三段
Pav_tm=sum（xt.^2）/N;%在时域计算信号功率
Pav_fn=sum（Px）;%通过功率谱计算信号功率
subplot（222）   %作出功率谱密度图
plot（0:3:fs10*log10（Px））;
xlable（‘频率（Hz）‘）;
Ylable（‘功率谱（dB）‘）;
title（‘分段估计得出的功率谱估计‘）;



%增加分段数可以进一步降低估计的协方差，然而每段中的数据点太少，就会使估计的
%频率分辨率下降很多。在样本信号序列总点数一定的条件下，可以采用使分段相互重叠
%的方法来增加分段数，从而保持每段中信号点数不变，这样就在保证频率分辨率的前提下
%进一步降低了估计