现代信号处理例的题目及matlab代码实现.docx

现代信号处理例的题目及matlab代码实现.docx

《现代信号处理例的题目及matlab代码实现.docx》由会员分享，可在线阅读，更多相关《现代信号处理例的题目及matlab代码实现.docx（15页珍藏版）》请在冰豆网上搜索。

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

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

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

程序代码如下：

clear;

clearall;

N=1024;%数据长度

Nfft=1024;%FFT所采用的数据长度

n=0:

N-1;

wn=randn（1,N）;%产生随机白噪声

subplot（2,2,1）;%绘出白噪声序列

plot（n,wn）;

title（'白噪声'）;

%500次实现的平均功率谱密度

s=zeros（1,N）;

fori=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次实现的平均功率谱密度'）;

gridon;

%1000次实现的平均功率谱密度

s=zeros（1,N）;

fori=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次实现的平均功率谱密度'）;

gridon;

%500次实现的平均功率谱密度

s=zeros（1,N）;

fori=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次实现的平均功率谱密度'）;

gridon;

实验结果图如下：

2仿真如下随机过程：

其中：

Vn是均值为0，方差为1的Gaussian白噪声过程，Φ为随机相位，在[0，2π]间服从均匀分布。

试对其中的正弦波频率进行估计（在不同的数据长度下，N=16，64，128，1024，可使用经典谱估计中的任何一种方法），并讨论数据长度对估计分辨率和平滑特性的影响。

解答：

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

程序代码如下：

clearall;

%********************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'）;

gridon;

%*******************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'）;

gridon;

%*******************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'）;

gridon;

%*******************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'）;

gridon;

实验结果如下：

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

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

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

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

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

①

②

频率范围[－2π，2π]

解答：

程序代码如下：

clearall;

closeall;

N=32;

dt=1;%设置最大点数

n=[0:

1:

N-1];

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的参数值

用两种方法：

1直接用MA方程：

2谱分解方法。

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

解答：

MA方程方法：

由

，

，c0=7.24,c1=3.6

得b0=1.3379，b1=2.6907。

程序代码如下：

b=[1.33792.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方程方法得到的谱估计'）;

谱分解方法：

程序代码如下：

clearall;

N=456;

B1=[12];

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）;

fori=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）;

fori=1:

2

r1（i,1）=-Ry（260+i）;

forj=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（'谱分解方法得到的谱估计'）;

实验结果图如下：