现代数字信号处理

现代数字信号处

理实验报告

姓名：杨航

学号：2017170707

班级：电子与通信40班

分数：

计算机作业一

题目

在计算机上用如下方法产生随机信号()

u n的观测样本：首先产生一段零均值、方差为2σ的复高斯白噪声序列()

v n；然后在()

v n上叠加三个复正弦信号，它们的归一化频率分别是f1=0.15，f2=0.17和f3=0.25。调整2σ和正弦信号的幅度，使在f1、f2和f3处得信噪比分别为30dB、30dB和27dB。

(a)令信号观测样本长度N=32，试估计出信号自相关函数

^ () 0m r

自相关函数是信号间隔的函数，间隔有正负间隔，所以n个长度的信号，有2n-1个自相关函数值，分别描述的是不同信号间隔的相似程度，并且该2n-1个自相关函数值关于n对称。

Matlab中用于计算自相关函数的指令是xcorr

其指令格式为

c= xcorr(x,'option')

其中选项为

"biased"为有偏的互相关函数估计；

"unbiased"为无偏的互相关函数估计；

"coeff"为0延时的正规化序列的自相关计算；

"none"为原始的互相关计算；

所以利用这个函数我们就可以得到信号的自相关函数

仿真结果为

（b）令信号观测样本长度N=256，试用BT法和周期图法估计()

u n的功率谱，这里设BT法中所用自相关函数的单边长度M=64。

BT法即经典功率谱估计法中的间接法，即先由N个观察值，估计出自相关函数，然后求自相关函数的傅里叶变换，以此变换结果作为对功率谱的估计。

直接法功率谱估计是间接法功率谱估计的一个特例，又称为周期图法，它是把随机信号的N个观察值，直接进行傅里叶变换，然后取其幅值的平方，再除以N，作为对功率谱的估计。

其仿真结果为

(c)令信号观测样本长度N=256，试用Levinson-Durbin迭代算法求解AR模

型的系数并估计()

u n的功率谱，模型的阶数取为p=16。

Levinson-Durbin迭代算法是为用来避免矩阵求逆运算，降低计算量的一种实用方法。

实际运算中LD算法是用自相关的函数估计值来得到AR模型的系数，然后实用FFT就可以计算出P阶AR模型的功率谱

图利用Levinson-Durbin迭代算法计算功率谱流程图

其仿真结果为

（d）利用上述3中方法得到的功率谱，分别估计出3个频率分量具体的归一化频率值，重复上述实验100次，统计频率值估计的均值和均方误差。

利用上诉两个公式我们就可以求出估计值的均值和均方误差

然后是统计实验数据

正弦信号在其频率值处有明显的波峰，所以频率谱图中中波峰处的频率值即为频率分量的具体归一化频率。

统计数据后可以得到数据表

表归一化频率均值、均方误差

（e）逐渐增加噪声功率，给出频率值估计的均方误差随着信噪比（dB）的变化曲线，比较3种算法的性能。

三种图像的均方误差变化曲线为

性能比较

1)周期图法得到的功率谱分辨率最高，但是方差性能最差，功率谱起伏剧烈，容易出现虚假谱峰；

2)BT法由于使用了平滑窗对直接法估计的功率谱进行了平滑，因此方差性能比周期图法好，功率谱比直接法估计的要平滑，但其分辨率比直接法低。

3)由于AR模型是一个有理分式，因而估计出的谱要比经典法的谱平滑。AR 谱估计的一些方法隐含着数据和自相关函数的外推，可获得高的分辨率

代码以及注释

（1）观测样本长度N=32，试估计出信号的自相关函数

%产生零均值、方差为1的复高斯白噪声序列

N=32;

noise=(randn(1,N)+j*randn(1,N))/sqrt(2);

%产生三个复正弦信号并产生观察样本

f1=0.15;%信号的归一化频率

f2=0.17;

f3=0.25;

SNR1=30;%信号的信噪比

SNR2=30;

SNR3=27;

A1=10^(SNR1/20);%信号的幅度

A2=10^(SNR2/20);

A3=10^(SNR3/20);

signal1=A1*exp(j*2*pi*f1*(0:N-1));%产生复正弦信号

signal2=A2*exp(j*2*pi*f2*(0:N-1));

signal3=A3*exp(j*2*pi*f3*(0:N-1));

un=signal1+signal2+signal3+noise;%观测样本信号u（n）

%基于FFT的自相关函数快速计算方法

Uk=fft(un, 2*N);%对u（n）进行2N点的FFT

Sk=(1/N)* abs(Uk).^2;%计算功率谱估计S（k）

r0=ifft(Sk);%对功率谱估计S（k）求逆FFT

r1=[r0(N+2:2*N),r0(1:N)];%根据教材式求得自相关函数

figure(1);

stem(real(r1));%自相关函数的实部

figure(2);

stem(imag(r1))%自相关函数的虚部

（2）令观测样本长度N=256，试分别用BT法和周期图法估计u(n)的功率谱，其中BT法所用自相关函数的单边长度M=64

%产生零均值、方差为1的复高斯白噪声序列

N=256;

noise=(randn(1,N)+j*randn(1,N))/sqrt(2);

%产生三个复正弦信号并产生观察样本

f1=0.15;%信号的归一化频率

f2=0.17;

f3=0.25;

SNR1=30;%信号的信噪比

SNR2=30;

SNR3=27;

A1=10^(SNR1/20);%信号的幅度

A2=10^(SNR2/20);

A3=10^(SNR3/20);

signal1=A1*exp(j*2*pi*f1*(0:N-1));%产生复正弦信号

signal2=A2*exp(j*2*pi*f2*(0:N-1));

signal3=A3*exp(j*2*pi*f3*(0:N-1));

un=signal1+signal2+signal3+noise;%观测样本信号u（n）

%BT法估计u（n）的功率谱

M=64;%自相关函数的单边长度

r2=xcorr(un, M, 'biased');%直接计算自相关函数

NF=1024;%BT法中FFT的点数

BT=fftshift(fft(r2, NF))%BT法计算功率谱

BT=BT/max(BT);%归一化用于画图

figure(1);

plot((-511:512)/1024,10*log10(BT));

%周期图法估计u（n）的功率谱

NF=1024; %周期图法中FFT的点数

Spr1=fftshift((1/NF)*abs(fft(un,NF)).^2)%计算信号的周期图

Spr= Spr1/max(Spr1);%归一化用于画图

t=[-0.5:1/NF:0.5-1/NF];

figure(2);

plot(t,10*log10(Spr));

（3）令信号观测样本长度N=256，试用Levinson-Durbin迭代算法求解AR模型系数并估计信号功率谱，阶数p=16

%产生零均值、方差为1的复高斯白噪声序列

N=256;

noise=(randn(1,N)+j*randn(1,N))/sqrt(2);

%产生三个复正弦信号并产生观察样本

f1=0.15;%信号的归一化频率

f2=0.17;

f3=0.25;

SNR1=30;%信号的信噪比

SNR2=30;

SNR3=27;

A1=10^(SNR1/20);%信号的幅度

A2=10^(SNR2/20);

A3=10^(SNR3/20);

signal1=A1*exp(j*2*pi*f1*(0:N-1));%产生复正弦信号

signal2=A2*exp(j*2*pi*f2*(0:N-1));

signal3=A3*exp(j*2*pi*f3*(0:N-1));

un=signal1+signal2+signal3+noise;%观测样本信号u（n）

%Levinson-Durbin迭代算法求解AR模型系数并估计信号功率谱

p=16;%AR模型的阶数

r0=xcorr(un, p, 'biased');%计算自相关函数r（0）

r=r0(p + 1 : 2 * p + 1);%提取r（0），r（1）...r(p)

%计算一阶AR模型的系数与输入方差

a(1, 1)=-r(2)/r(1);%1阶AR模型系数

sigma(1)=r(1)-(abs(r(2))^2)/r(1);%一阶AR模型的输入方差

%Levinson-Durbin迭代算法

for m=2:p

k(m)= -(r(m+1)+sum(a(m-1,1:m-1).*r(m:-1:2)))/sigma(m-1);%反射系数

a(m,m)=k(m);

for i=1:m-1

a(m, i)=a(m-1, i) + k(m) * conj(a(m-1, m-i)); %m阶AR的系数

end

sigma(m)=sigma(m-1) * (1-abs(k(m))^2); %m阶AR模型的方差

end

%计算十六阶AR模型的功率谱

NF = 1024;%AR方法中FFT的点数

Par = sigma(p) ./ fftshift(abs(fft([1, a(p,:)], NF)) .^2);%p劫AR模型的输入方差

Par=Par/max(Par);%归一化用于画图

plot((-511:512)/1024,10*log10(Par));

（4）利用上述3中方法得到的功率谱，分别估计出3个频率分量具体的归一化频率值，重复上述实验100次，统计频率值估计的均值和均方误差。（是的第k次估计值）function res = test(sig)

%产生零均值、方差为1的复高斯白噪声序列

N=256;

noise=sig * (randn(1,N)+1i*randn(1,N))/sqrt(2);

%产生三个复正弦信号并产生观察样本

f1=0.15;

f2=0.17;

f3=0.25;%信号的归一化频率

SNR1=30;

SNR2=30;

SNR3=27;

%信号的信噪比

A1=10^(SNR1/20);

A2=10^(SNR2/20);

A3=10^(SNR3/20);%信号的幅度

signal1=A1*exp(1i*2*pi*f1*(0:N-1));%产生复正弦信号

signal2=A2*exp(1i*2*pi*f2*(0:N-1));

signal3=A3*exp(1i*2*pi*f3*(0:N-1));

un=signal1+signal2+signal3+noise;%观测样本信号u（n）

p=16;r0=xcorr(un,p,'biased');r=r0(p+1:2*p+1);%计算自相关函数r（0）

a(1,1)=-r(2)/r(1);%1阶AR模型系数

sigma(1)=r(1)-(abs(r(2))^2)/r(1);%一阶AR模型的输入方差

for m=2:p;

k(m)=-(r(m+1)+sum(a(m-1,1:m-1).*r(m:-1:2)))/sigma(m-1);%反射系数

a(m,m)=k(m);

for i=1:m-1;

a(m,i)=a(m-1,i)+k(m)*conj(a(m-1,m-i));%m阶AR的系数

end

sigma(m)=sigma(m-1)*(1-abs(k(m))^2); %m阶AR模型的方差

end

%计算十六阶AR模型的功率谱

NF=1024;%AR方法中FFT的点数

m1=-511:512;

Par=sigma(p)./fftshift(abs(fft([1,a(p,:)],NF)).^2);%p阶AR模型的功率谱%统计十六阶AR模型的频率值估计的均值和误差估计

res = 10*log10(Par)/max(10*log10(Par));

res = [res;m1]';

res = sortrows(res,1);

res = res(end-2:end,:);

res = sort(res(:,2));

end

clear all;

for sig=1:100%重复100次

temp = 0;

temp1 = [];

for i = 1

temp = test(sig)'./1000;%3个功率谱的均值估计

temp1 = [temp1; temp];%3个功率谱的误差估计

end

temp1

avg = mean(temp1);%均值估计

v = sum(var(temp1))/3;%均方误差估计

end

function res = test(sig)

%产生零均值、方差为1的复高斯白噪声序列

N=256;

noise=sig * (randn(1,N)+1i*randn(1,N))/sqrt(2);

%产生三个复正弦信号并产生观察样本

f1=0.15;

f2=0.17;

f3=0.25;%信号的归一化频率

SNR1=30;

SNR2=30;

SNR3=27;

%信号的信噪比

A1=10^(SNR1/20);

A2=10^(SNR2/20);

A3=10^(SNR3/20);%信号的幅度

signal1=A1*exp(1i*2*pi*f1*(0:N-1));%产生复正弦信号

signal2=A2*exp(1i*2*pi*f2*(0:N-1));

signal3=A3*exp(1i*2*pi*f3*(0:N-1));

un=signal1+signal2+signal3+noise;%观测样本信号u（n）

p=16;r0=xcorr(un,p,'biased');r=r0(p+1:2*p+1);%计算自相关函数r（0）

a(1,1)=-r(2)/r(1);%1阶AR模型系数

sigma(1)=r(1)-(abs(r(2))^2)/r(1);%一阶AR模型的输入方差

for m=2:p;

k(m)=-(r(m+1)+sum(a(m-1,1:m-1).*r(m:-1:2)))/sigma(m-1);%反射系数

a(m,m)=k(m);

for i=1:m-1;

a(m,i)=a(m-1,i)+k(m)*conj(a(m-1,m-i));%m阶AR的系数

end

sigma(m)=sigma(m-1)*(1-abs(k(m))^2); %m阶AR模型的方差

end

%计算十六阶AR模型的功率谱

NF=1024;%AR方法中FFT的点数

m1=-511:512;

Par=sigma(p)./fftshift(abs(fft([1,a(p,:)],NF)).^2);%p阶AR模型的功率谱

%统计十六阶AR模型的频率值估计的均值和误差估计

res = 10*log10(Par)/max(10*log10(Par));

res = [res;m1]';

res = sortrows(res,1);

res = res(end-2:end,:);

res = sort(res(:,2));

end

（5）逐渐增加噪声功率，给出频率值估计的均方误差随着噪声功率的变化曲线，比较3中算法的性能。

clear all;

y = [];

for sig=1:100%重复100次

temp = 0;

temp1 = [];

for i = 1:5%增加噪声功率

temp = test(sig)'./1000;%3个功率谱的均值估计

temp1 = [temp1; temp];%3个功率谱的误差估计

end

temp1

avg = mean(temp1);%均值估计

v = sum(var(temp1))/3;%均方误差估计

y = [y,v];

end

y

plot(y);

function res = test(sig)

%产生零均值、方差为1的复高斯白噪声序列

N=256;

noise=sig * (randn(1,N)+1i*randn(1,N))/sqrt(2);

%产生三个复正弦信号并产生观察样本

f1=0.15;

f2=0.17;

f3=0.25;%信号的归一化频率

SNR1=30;

SNR2=30;

SNR3=27;

%信号的信噪比

A1=10^(SNR1/20);

A2=10^(SNR2/20);

A3=10^(SNR3/20);%信号的幅度

signal1=A1*exp(1i*2*pi*f1*(0:N-1));%产生复正弦信号

signal2=A2*exp(1i*2*pi*f2*(0:N-1));

signal3=A3*exp(1i*2*pi*f3*(0:N-1));

un=signal1+signal2+signal3+noise;%观测样本信号u（n）

p=16;r0=xcorr(un,p,'biased');r=r0(p+1:2*p+1);%计算自相关函数r（0）a(1,1)=-r(2)/r(1);%1阶AR模型系数

sigma(1)=r(1)-(abs(r(2))^2)/r(1);%一阶AR模型的输入方差

for m=2:p;

k(m)=-(r(m+1)+sum(a(m-1,1:m-1).*r(m:-1:2)))/sigma(m-1);%反射系数a(m,m)=k(m);

for i=1:m-1;

a(m,i)=a(m-1,i)+k(m)*conj(a(m-1,m-i));%m阶AR的系数

end

sigma(m)=sigma(m-1)*(1-abs(k(m))^2); %m阶AR模型的方差

end

%计算十六阶AR模型的功率谱

NF=1024;%AR方法中FFT的点数

m1=-511:512;

Par=sigma(p)./fftshift(abs(fft([1,a(p,:)],NF)).^2);%p阶AR模型的功率谱%统计十六阶AR模型的频率值估计的均值和误差估计

res = 10*log10(Par)/max(10*log10(Par));

res = [res;m1]';

res = sortrows(res,1);

res = res(end-2:end,:);

res = sort(res(:,2));

end

计算机作业 三

一、题目

考虑一个由如下差分方程定义的AR 过程x(n)

12()(1)(2)()x n a x n a x n v n =----+

其中v(n)是零均值，方差为 2

v σ的加性白噪声，其中参数a1和a2均为实数：

(a) 当x(n)为单位方差时计算噪声方差 2

v σ 。并给出不同实现的x(n)。

(b) 给定输入信号x(n)，用两阶LMS 滤波器估计参数和的值。画出100次试验中和的平均值收敛曲线，迭代步长为0.04。计算时间常数，并与理论值比较。 (c)在一次实现中，计算预测误差()()()f n x n x n =-，以及两个系数的

的误差111()()n a h

n ε=--222()()n a h n ε=--。分别画出法f （n ）的能量谱。

(d) 计算100次不同实现中预测误差平方的平均学习曲线。根据试验结果计算时间常数和残留能量，并与理论值比较。

二、理论分析

二阶线性预测滤波器，其输出为x(n)的预测值并可表示为：

预测误差为： 对于二阶自回归模型，LMS 滤波器的抽头系数，也即待预测参数:

w1(n+1)=w1(n)+2*u*e(n)*x(n-1) w2(n+1)=w2(n)+2*u*e(n)*x(n-2)

收敛准则：线性预测滤波器采用LMS 算法根据一定的准则，按照一定的步长u 调整w1和w2，使均方误差E{e(n)*e(n)}尽可能最小。

LMS 算法步骤

设置变量和参数： X(n):输入的加噪信号 v(n):是零均值的白噪声

W(n):LMS 滤波器的抽头系数，也即待预测参数 d(n):期望信号（即真实输入信号x(n)） e(n ）:误差信号序列

u:步长,u 的选择和收敛速度收敛曲线的稳定性有关，

u 应满足：rho_max = max(eig(x*x.')); % 输入信号相关矩阵的最大特征值

mu = rand()*(1/rho_max) ; % 收敛因子 0 < mu < 1/rho_max

n_times:一次独立实验的迭代次数上限 K:k 次独立实验，求均值。

)2()1()(?21-+-=n x w n x w n x

)(?)()(n x

n x n e -=

三、仿真结果

(a)

当x(n)为单位方差时计算噪声方差 。并给出不同实现的x(n)。 因为x(n)为单位方差，所以 =1，r3=inv([1,-a1,-a2;0,-(1+a2),0;0,-a1,-1])*[1;a1;a2];%inv 是求逆矩阵，x ’是求转置矩阵，当为复数矩阵时，求的是共轭转置

R2=[r3(1),r3(2),r3(3);r3(2),r3(1),r3(3);r3(3),r3(2),r3(1)]; sigma_v=r3(1); r1=r3(2); r2=r3(3);

(b)

假定a1=-1.6966;a2=0.88，产生N=1000的样本序列。

100

200

300

400

500

600

700

800

900

1000

输入信号x

图输入信号

(c)

图u=0.04 a1及a2的预测曲线

图f（n)的预测误差曲线和能量谱

图a1的预测误差曲线和能量谱

图a2的预测误差曲线和能量谱