中国科学院刘艳老师现代数字信号处理第四章上机作业.pdf

1、假设一平稳随机信号为 w n n x 1 + ? 0.8 x n ) ( = ( ) ( )，其中 w(n) 是均值为 0，方差 为 1 的白噪声，数据长度为 1024。 （1）、产生符合要求的 w(n) 和 x(n) ； （2）、给出信号 x(n)的理想功率谱； （3）、编写周期图谱估计函数，估计数据长度 N=1024 及 256 时信号功率谱， 分析估计效果。 （4）、编写 Bartlett 平均周期图函数，估计当数据长度 N=1024 及 256 时，分 段数 L 分别为 2 和 8 时信号 x(n) 的功率谱，分析估计效果。 一、解题思路 w(n)可以通过随机序列randn(1,N)来产生，x(n)可以通过对w(n)滤波产生（由 递推式可得系统的传递函数），也可以直接由递推式迭代产生。 由于线性系统的输出功率谱等于输入功率谱乘以传递函数模的平方，X(n)可 以看做w(n)通过一线性系统的输出，H(z)=1/(1-0.8z)。所以x(n)的理想功率谱 P(ejw)=σw2|H(ejw)|2。 周期图方法：直接对观测数据做FFT变换，变换的结果取模的平方再除以数据 长度，作为估计的功率谱。256个观测点时可以对原观测数据以4为间隔提取得到。 Bartlett法：将L组独立的观测数据分别求周期图，再将L个周期图求平均作 为信号的功率谱估计。L组数据可以通过对原观测数据以L为间隔提取得到。 二、MATLAB实现程序及注解 clear all; 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; 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