《现代信号分析与处理》课终作业(20090709)

《现代信号分析与处理》课终作业

姓名：司爱威

年级：2009级博士研究生

学号：111200911005

专业：车辆工程

2009-7-16

一、线性卷积程序

1.1 子函数程序（matlab语言）

function y=saw_fun_conv(x,h) % x,h分别表示两个要卷积的信号Nx=length(x);

Nh=length(h);

Ny=Nx+Nh-1;

for i=Nx+1:Ny

x(i)=0;

end

for i=Nh+1:Ny

h(i)=0;

end

for i=1:Ny

sum=0;

for j=1:i

sum = sum+ x(j)*h(i-j+1);

end

y(i)=sum;

end

1.2 验证程序

x=[1 2 3 4];

h=[1 1];

y=saw_fun_conv(x,h)

subplot(2,2,1);

stem(x,'bo');xlabel('x(k)');grid on;

subplot(2,2,3);

stem(h,'bo');axis([1 4 0 5]);xlabel('h(k)');grid on;

subplot(1,2,2);

stem(y,'bo');xlabel('y(n)');grid on;

1.3 验证波形

如图1.1所示。其中

x=[1 2 3 4];

h=[1 1];

y=[1 3 5 7 4].

图1.1 线性卷积程序验证波形

二、DFT程序

2.1 子函数程序

function X=saw_fun_dft(x) % x表示原信号

N=length(x);

W=exp(-j*2*pi/N);

for k=1:N

sum=0;

for n=1:N

sum=sum+x(n)*W.^((n-1)*(k-1));

end

Xa(k)=abs(sum);

end

for i=1:N/2

X(i)=Xa(i);

end

2.2 验证程序

N=256;

f1=0.1;f2=0.2;fs=1;

a1=5;a2=3;

w=2*pi/fs;

x=a1*sin(w*f1*(0:N-1))+a2*sin(w*f2*(0:N-1))+randn(1,N);

X=saw_fun_dft(x);

subplot(2,1,1);

plot(x);xlabel({'n','(a)'});ylabel('x(n)');grid on;

f=0:1/N:fs/2-1/N;

subplot(2,1,2);

plot(f,X);xlabel({'f/Hz','(b)'});ylabel('X(f)');grid on;

2.3 验证波形

如图2.1所示。其中

x=5sin(0.2πn)+3sin(0.4πn) +randn(1,256), n=0,2, (255)

图（a）为原信号波形，图（b）为DFT变换后的波形。

图2.1 DFT程序验证波形

三、窗函数法设计FIR滤波器程序

3.1 子函数程序

function h=saw_fun_firwin(N,fc) %低通滤波器，选用汉明窗，N表示滤波器阶数，

fc表示截止频率。

if rem(N,2)==0

n2=N/2-1;

flag=1;

else

n2=N/2;

flag=0;

end

wc1=2*pi*fc;

for i=0:n2

s=i-N/2;

h(i+1)=(sin(wc1*s)/(pi*s))*(0.54-0.46*cos(2*i*pi/(N+1)));

h(N-i+1)=h(i+1);

end

if(flag==1)

h(N/2+1)=wc1/pi;

end

3.2 验证程序

N=60;

fc=0.15;

f1=0.1;f2=0.2;fs=1;

a1=5;a2=3;

w=2*pi/fs;

x=a1*sin(w*f1*(0:N-1))+a2*sin(w*f2*(0:N-1))+randn(1,N);

X=saw_fun_dft(x); %利用第2题DFT程序

N2=20;

h=saw_fun_firwin(N2,fc)

x2=saw_fun_conv(x,h); %利用第1题线性卷积程序

X2=saw_fun_dft(x2); %利用第2题DFT程序

subplot(321);

stem(h,'.');xlabel({'n','(a)'});ylabel('h(n)');grid on;

subplot(322);

f11=0:0.5/N:0.5-0.5/N;

plot(f11,abs(freqz(h,1,N)));xlabel({'f','(b)'});ylabel('H(f)');grid on;

subplot(323);

plot(x);xlabel({'n','(c)'});ylabel('x(n)');grid on;

subplot(324);

f12=0:1/N:0.5-1/N;

plot(f12,X);xlabel({'f','(d)'});ylabel('X(f)');grid on;

subplot(325);

plot(x2);xlabel({'n','(e)'});ylabel('x2(n)');grid on;

subplot(326);

f13=0:1/80:0.5-1/80;

plot(f13,X2);xlabel({'f','(f)'});ylabel('X2(f)');grid on;

3.3 验证波形

如图3.1所示。其中

x=5sin(0.2πn)+3sin(0.4πn) +randn(1,60), n=0,2, (255)

图（a）为滤波器系数；图（b）为滤波器幅频响应曲线；图（c）为信号x的时域波形；图（d）为信号x的频谱；图（e）为信号x通过滤波器后的时域波形；图（f）为信号x通过滤波器后的频谱。

图3.1 窗函数法设计FIR滤波器程序验证波形

四、Welch法功率谱估计程序

4.1 子函数程序

function P=saw_fun_welch(x,M) %选用汉明窗，x表示原信号，M表示分段长度。N=length(x);

K=(N-M/2)/(M/2);

s=0;

m2=M/2;

for k=0:K-1

for i=0:M-1

xa(k+1,i+1)=x(s+i+1);

end

s=s+m2;

end

U=0;

for k=0:K-1

for i=0:M-1

win(i+1)=0.54-0.46*cos(2*pi*i/(M-1));

xin(k+1,i+1)=xa(k+1,i+1)*win(i+1);

U=U+win(i+1)*win(i+1);

end

end

for k=0:K-1

X=saw_fun_dft(xin(k+1,1:M)); %利用第2题DFT程序

for i=0:m2-1

S(k+1,i+1)=(abs(X(i+1)))*(abs(X(i+1)));

end

end

for i=0:m2-1

Sxx(i+1)=0;

for k=0:K-1

Sxx(i+1)=Sxx(i+1)+S(k+1,i+1);

end

P(i+1)=Sxx(i+1)/(K*U);

end

4.2 验证程序

N=256;

M=64;

f1=0.1;f2=0.2;fs=1;

a1=5;a2=3;

w=2*pi/fs;

x=a1*sin(w*f1*(0:N-1))+a2*sin(w*f2*(0:N-1)) +randn(1,N);

P=saw_fun_welch(x,M);

subplot(2,1,1);

plot(x);xlabel({'n','(a)'});ylabel('x(n)');grid on;

f=0:1/M:fs/2-1/M;

subplot(2,1,2);

plot(f,P);xlabel({'f/Hz','(b)'});ylabel('P(f)');grid on;

4.3 验证波形

如图4.1所示。其中

x=5sin(0.2πn)+3sin(0.4πn)+ randn(1,256), n=0,1,2,…,255

图（a）为原信号x(n)的波形，图（b）为信号x(n)的功率谱。

图4.1 Welch法功率谱估计程序验证波形

五、Burg法AR谱估计程序

5.1 子函数程序

function P=saw_fun_arburg(x,order) % x表示原信号，order表示AR模型阶数N=length(x);

a(1)=1;

r0=0;

for i=0:N-1

r0=r0+x(i+1)*x(i+1)/N;

end

v(1)=r0;

for i=1:N-1

ef(i+1)=x(i+1);

eb(i)=x(i);

end

for k=1:order

sumn=0;

sumd=0;

for i=k:N-1

sumn=sumn+ef(i+1)*eb(i);

sumd=sumd+ef(i+1)*ef(i+1)+eb(i)*eb(i);

end

a(k+1)=-2*sumn/sumd;

for i=1:k-1

b(i+1)=a(i+1)+a(k+1)*a(k-i+1);

end

for i=1:k-1

a(i+1)=b(i+1);

end

v(1)=(1-a(k+1)*a(k+1))*v(1);

end

W=exp(-j*2*pi/N);

for k=1:N

sum=0;

for n=2:order

sum=sum+x(n)*W.^((n-1)*(k-1));

end

Xa(k)=abs(sum);

end

for k=N/2+1:N

S(k-N/2)=Xa(k);

P(k-N/2)=v/(abs(1+S(k-N/2))*abs(1+S(k-N/2)));

end

5.2 验证程序

N=256;

f1=0.1;f2=0.2;fs=1;

a1=5;a2=3;

w=2*pi/fs;

x=a1*sin(w*f1*(0:N-1))+a2*sin(w*f2*(0:N-1))+randn(1,N);

order=12;

P=saw_fun_arburg(x,order);

subplot(2,1,1);

plot(x);xlabel({'n','(a)'});ylabel('x(n)');grid on;

f=0:1/N:fs/2-1/N;

subplot(2,1,2);

plot(f,P);xlabel({'f/Hz','(b)'});ylabel('P(f)');grid on;

5.3 验证波形

如图5.1所示。其中

x=5sin(0.2πn)+3sin(0.4πn)+ randn(1,256), n=0,2, (255)

图（a）为原信号x(n)的波形，图（b）为信号x(n)的功率谱。

图5.1 Burg法功率谱估计程序验证波形