现代信号处理例题及matlab代码实现

《现代信号处理》期末考核作业

1 MATLAB仿真均值为0，方差为1的白噪声信号，信号长度N=1024，并用周期图法分别求500、1000和1500次实现的平均功率谱密度，画图。

程序代码如下：

clear; clear all; N=1024;%数据长度

Nfft=1024;?T所采用的数据长度 n=0:N-1;

wn=randn(1,N);%产生随机白噪声 subplot(2,2,1);%绘出白噪声序列 plot(n,wn); title('白噪声');

P0次实现的平均功率谱密度 s=zeros(1,N); for i=1:500

wn=randn(1,N);%产生随机白噪声

Pxx=10*log10(abs(fft(wn,Nfft).^2)/N);%Fourier振幅谱平方的平均值，并转换为db s=s+Pxx; end s=s/500;

f=(0:length(Pxx)-1)/length(Pxx);%绘出频率序列 subplot(222); plot(f,s);

xlabel('频率/Hz');ylabel('功率谱/dB'); title('500次实现的平均功率谱密度'); grid on;

00次实现的平均功率谱密度 s=zeros(1,N); for i=1:1000

wn=randn(1,N);%产生随机白噪声

Pxx=10*log10(abs(fft(wn,Nfft).^2)/N);%Fourier振幅谱平方的平均值，并转换为db s=s+Pxx; end s=s/1000;

f=(0:length(Pxx)-1)/length(Pxx);%绘出频率序列 subplot(223); plot(f,s);

xlabel('频率/Hz');ylabel('功率谱/dB'); title('1000次实现的平均功率谱密度'); grid on;

P0次实现的平均功率谱密度 s=zeros(1,N); for i=1:1500

wn=randn(1,N);%产生随机白噪声

Pxx=10*log10(abs(fft(wn,Nfft).^2)/N);%Fourier振幅谱平方的平均值，并转换为db s=s+Pxx; end s=s/1500;

f=(0:length(Pxx)-1)/length(Pxx);%绘出频率序列 subplot(224); plot(f,s);

xlabel('频率/Hz');ylabel('功率谱/dB'); title('1500次实现的平均功率谱密度'); grid on;

实验结果图如下：

2仿真如下随机过程：

xn?sin(n+φ)?Vn?sin(n+φ)

23其中：Vn是均值为0，方差为1的Gaussian白噪声过程，Φ为随机相位，在[0，2π]间服从均匀分布。试对其中的正弦波频率进行估计（在不同的数据长度下，N=16，64，128，1024，可使用经典谱估计中的任何一种方法），并讨论数据长度对估计分辨率和平滑特性的影响。

解答：使用周期图法对不同数据长度的信号进行估计。

程序代码如下：

clear all;

%********************N=16;第一种情况数据长度为16********************** N = 16; Nfft = 16; n = 0:N-1;

xn = sin(0.5*pi*n+2*pi*rand)+sin(0.3333*pi*n+2*pi*rand)+randn(1,N); figure(1); subplot(2,1,1);

??

plot(n,xn); ylabel('幅值(V)'); xlabel('时间(s)'); title('原始信号');

Pxx=10*log10(abs(fft(xn,Nfft).^2)/N); Fourier振幅谱平方的平均值，并转换为db

f = 0:length(Pxx)-1; 绘出频率序列 subplot(212);plot(f,Pxx);

xlabel('频率/Hz');ylabel('功率谱/dB'); title('周期图N=16'); grid on;

%*******************N=64********************* N = 64; Nfft = 64; n = 0:N-1;

xn = sin(0.5*pi*n+2*pi*rand)+sin(0.3333*pi*n+2*pi*rand)+randn(1,N); figure(2); subplot(2,1,1);

plot(n,xn); ylabel('幅值(V)'); xlabel('时间(s)'); title('原始信号');

Pxx=10*log10(abs(fft(xn,Nfft).^2)/N); f = 0:length(Pxx)-1; subplot(212);plot(f,Pxx);

xlabel('频率/Hz');ylabel('功率谱/dB'); title('周期图N=64'); grid on;

%*******************N=64********************* N = 128;

Nfft = 128; n = 0:N-1;

xn = sin(0.5*pi*n+2*pi*rand)+sin(0.3333*pi*n+2*pi*rand)+randn(1,N); figure(3); subplot(2,1,1);

plot(n,xn); ylabel('幅值(V)'); xlabel('时间(s)'); title('原始信号');

Pxx = 10*log10(abs(fft(xn,Nfft).^2)/N); f = 0:length(Pxx)-1; subplot(212);plot(f,Pxx);

xlabel('频率/Hz');ylabel('功率谱/dB'); title('周期图N=128'); grid on;

%*******************N=64********************* N = 1024; Nfft =1024; n = 0:N-1;

xn = sin(0.5*pi*n+2*pi*rand)+sin(0.3333*pi*n+2*pi*rand)+randn(1,N); figure(4); subplot(2,1,1);

plot(n,xn); ylabel('幅值(V)'); xlabel('时间(s)'); title('原始信号');

Pxx = 10*log10(abs(fft(xn,Nfft).^2)/N); f = 0:length(Pxx)-1; subplot(212);plot(f,Pxx);

xlabel('频率/Hz');ylabel('功率谱/dB');

title('周期图N=1024'); grid on;

实验结果如下：

数据长度N=16时的谱估计结果

数据长度N=64时的谱估计结果

数据长度N=128时的谱估计结果

数据长度N=1024时的谱估计结果

3计算并画图描绘下列函数离散时间傅立叶变换的幅度（幅频）特性：

π/5)?exp(in① x[n]??

0?π/5)?sin(nx[n]?② ?0?频率范围[－2π，2π] 解答： 程序代码如下：

clear all; close all; N=32;

dt=1;%设置最大点数 n=[0:1:N-1];

n=0，1…31 n=其它 n=0，1…31 n=其它

k=[0:1:N-1]; WN=exp(-j*2*pi/N); nk=n'*k;

WNnk=WN.^nk;%离散Fourier变换矩阵 t=n*dt;

%**************************xn信号 xn=exp(j*n*pi/5); figure(1); subplot(311); plot(t,xn); title('xn信号');

Xk=xn*WNnk;%对xn进行Fourier变换 magXk=abs(Xk); phaXk=angle(Xk); k=0:length(magXk)-1; subplot(312);

plot(k/(N*dt),magXk*2/N); title('xn信号的振幅谱'); subplot(313);

plot(k/(N*dt),unwrap(phaXk)); ylabel('相位角/rad'); title('xn信号的相位谱');

%**************************yn信号 yn=sin(n*pi/5); figure(2); subplot(311); plot(t,yn); title('yn信号');

Yk=yn*WNnk;%对xn进行Fourier变换 magYk=abs(Yk);

phaYk=angle(Yk); k=0:length(magYk)-1; subplot(312);

plot(k/(N*dt),magYk*2/N); title('yn信号的振幅谱'); subplot(313);

plot(k/(N*dt),unwrap(phaYk)); ylabel('相位角/rad'); title('yn信号的相位谱');

实验结果图如下：

信号①的离散时间傅里叶变换

信号②的离散时间傅里叶变换

4设Rxx(0)=7.24，设Rxx(1)=3.6，试确定如下一阶MA的参数值

H(z)?b0?b1Z?1

用两种方法：

?p?Rxx(0)??c0???bqbq?r??????其中c?① 直接用MA方程：??q?1r?????0??cp?Rxx(p)??????② 谱分解方法。

最后描绘用MA模型得到的谱估计。

解答：

r?pr?p

MA方程方法：

由c0?b12，c1?b0b1，c0=7.24, c1=3.6

得b0=1.3379，b1=2.6907。

程序代码如下：

b=[1.3379 2.6907]; %MA系统系数 w=linspace(0,pi,512);

H=freqz(b, w);%产生信号的频域响应 Ps=abs(H).^2;%计算得到功率谱 plot(w./(2*pi),Ps); xlabel('频率/Hz'); ylabel('功率谱/dB');

title('直接用MA方程方法得到的谱估计');

谱分解方法： 程序代码如下：

clear all; N=456; B1=[1 2]; A1=[1];

w=linspace(0,pi,512);

%采用自协方差法对AR模型参数进行估计%

y1=filter(B1,A1,randn(1,N)).*[zeros(1,200),ones(1,256)]; [Py11,F]=pcov(y1,1,512,1); %AR(1)的估计% [Py13,F]=periodogram(y1,[],512,1); %****MA模型****% y=zeros(1,256); for i=1:256 y(i)=y1(200+i); end

ny=[0:255];

z=fliplr(y);nz=-fliplr(ny);

nb=ny(1)+nz(1);ne=ny(length(y))+nz(length(z)); n=[nb:ne]; Ry=conv(y,z); R1=zeros(2,2); r1=zeros(2,1); for i=1:2

r1(i,1)=-Ry(260+i); for j=1:2

R1(i,j)=Ry(260+i-j); end end R1; r1;

a1=inv(R1'*R1)*R1'*r1; %利用最小二乘法得到的估计参数 %对MA的参数进行估计% A1;

A11=[1,a1']; %AR的参数的估计值

B11=fliplr(conv(fliplr(B1),fliplr(A11))); %MA模型的分子

y21=filter(B11,A1,randn(1,N));%.*[zeros(1,200),ones(1,256)]; %由估计出的MA模型产生数据

[Ama1,Ema1]=arburg(y21,32);

B1; b1=arburg(Ama1,1); %求出MA模型的参数 %---求功率谱---% w=linspace(0,pi,512);

H11=freqz(b1,A11,w);%产生信号的频域响应 figure

plot(w./(2*pi), abs(H11).^2);

xlabel('频率/Hz');ylabel('功率谱/dB'); title('谱分解方法得到的谱估计');

实验结果图如下：

本文来源：https://www.bwwdw.com/article/7ygo.html