matlab经典、现代功率谱估计

上机作业：

1、假设一平稳随机信号为()()()0.81x n x n w n =−+，其中 是均值为0，方差为1的白噪声，数据长度为1024。

（1）、产生符合要求的)(n w 和)(n x ；

（2）、给出信号)(n x 的理想功率谱；

（3）、编写周期图谱估计函数，估计数据长度N=1024及256时信号功率谱，分析估计效果。

（4）、编写Bartlett 平均周期图函数，估计当数据长度N=1024及256时，分段数L 分别为2和8时信号 的功率谱，分析估计效果。

一、解题思路

w(n)可以通过随机序列randn(1,N)来产生，x(n)可以通过对w(n)滤波产生（由递推式可

得系统的传递函数），也可以直接由递推式迭代产生。

由于线性系统的输出功率谱等于输入功率谱乘以传递函数模的平方，X(n)可以看做w(n)通过一线性系统的输出，H(z)=1/(1-0.8z)。所以x(n)的理想功率谱P(e jw )=σw 2|H(e jw )|2。

周期图方法：直接对观测数据做FFT 变换，变换的结果取模的平方再除以数据长度，作为估计的功率谱。256个观测点时可以对原观测数据以4为间隔提取得到。

Bartlett 法：将L 组独立的观测数据分别求周期图，再将L 个周期图求平均作为信号的功率谱估计。L 组数据可以通过对原观测数据以L 为间隔提取得到。 二、MATLAB 实现程序及注解 clc;

clear;close all;

Fs=500; %采样率

N=1024; %观测数据

w=sqrt(1)+randn(1,N); %0均值,方差为1的白噪声,长度1024

x=[w(1) zeros(1,N-1)]; %初始化x(n),长度1024,x(1)=w(1)

for i=2:N

x(i)=0.8*x(i-1)+w(i); %迭代产生观测数据x(n)

end

%% 理想功率谱

[h,w1]=freqz(x);

figure,plot(w1*500/(2*pi),10*log10(abs(h).^2));grid on;

title('理想功率谱');

xlabel('频率'); ylabel('功率db');

%% 周期图法

%1024个观测点

Pxx=abs(fft(x)).^2/N; %周期图公式

Pxx=10*log10(Pxx(index+1)); %化为db

figure;plot(k,Pxx);grid on;

title('周期图1024点');

xlabel('频率'); ylabel('功率db');

% 周期图256个观测点

x1=x(1:4:N);

Pxx1=abs(fft(x1,1024)).^2/N;

Pxx1=10*log10(Pxx1(index+1)); %化为db

figure;plot(k,Pxx1);grid on;

title('周期图256点');

xlabel('频率'); ylabel('功率db');

%% Bartlett平均周期图,N=1024

%分段L=2

L=2;

x_21=x(1:L:N);

x_22=x(2:L:N);

Pxx_21=abs(fft(x_21,1024)).^2/length(x_21);

Pxx_22=abs(fft(x_22,1024)).^2/length(x_22);

Pxx_2=(Pxx_21+Pxx_22)/L;

figure;

subplot(2,2,1),plot(k,10*log10(Pxx_2(index+1)));grid on;

title('N=1024,L=2');

xlabel('频率'); ylabel('功率db');

%分段L=8

L1=8;

x3=zeros(L1,N/L1); %产生L1行,N/L1列的矩阵用以存储分组的数据

for i=1:L1

x3(i,:)=x(i:L1:N); %将原始数据分为8组

end

Pxx3=zeros(L1,1024); %产生L1行,1024列矩阵用以存储分组的周期图

for i=1:L1

Pxx3(i,:)=abs(fft(x3(i,:),1024)).^2/length(x3(i,:)); %分别求周期图,结果保存在Pxx3中,FFT长度为1024

end

for i=1:1024

Pxx3_m(i)=sum(Pxx3(:,i))/L1; %求平均

end

subplot(2,2,2),plot(k,10*log10(Pxx3_m(index+1)));grid on;

title('N=1024,L=8');

xlabel('频率'); ylabel('功率db');

%% Bartlett平均周期图,N=256,求法同上

%分段L=2,分别计算周期图,再取平均

x=x(1:4:N);

L2=2;

x_31=x(1:L2:length(x));

x_32=x(2:L2:length(x));

Pxx_31=abs(fft(x_31,1024)).^2/length(x_31);

Pxx_32=abs(fft(x_32,1024)).^2/length(x_32);

Pxx_3=(Pxx_31+Pxx_32)/L2;

subplot(2,2,3),plot(k,10*log10(Pxx_3(index+1)));grid on;

title('N=256,L=2');

xlabel('频率'); ylabel('功率db');