现代信号处理例题及matlab代码实现.docx
《现代信号处理例题及matlab代码实现.docx》由会员分享，可在线阅读，更多相关《现代信号处理例题及matlab代码实现.docx（16页珍藏版）》请在冰豆网上搜索。
现代信号处理例题及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（'谱分解方法得到的谱估计'）;
实验结果图如下：
升级会员 