DSP参考程序.doc

参考程序

%－－－－－实验1－－－－－
%－－－－－ex1_1－－－－－

clear all; close all; clc;
n = [0:20];
x1 = n.*(stepseq(0,0,20)-stepseq(10,0,20));
x2 = 10*exp(-0.3*(n-10)).*(stepseq(10,0,20)-stepseq(20,0,20));
x = x1+x2;
subplot(2,1,1);stem(n,x);
xlabel('n');ylabel('x(n)');axis([0,20,-1,11])
n=[-10:9];
x=[5,4,32,1];
xtilde=x' * ones(1,4);
xtilde=(xtilde(:))';
subplot(2,1,2);stem(n,xtilde);
xlabel('n');ylabel('xtilde(n)');axis([-10,9,-1,6])
set(gcf,'color','w');

%－－－－－实验1－－－－－
%－－－－－ex1_2－－－－－

clear all; close all; clc;
n=-4:2;x=[1,-2,4,6,-5,8,10];
[x21,n21]=sigshift(5*x,n,-5);
[x22,n22]=sigshift(4*x,n,-4);
[x23,n23]=sigadd(x21,n21,x22,n22);
[x2,n2]=sigadd(x23,n23,3*x,n);
subplot(2,1,1);stem(n2,x2);
xlabel('n');ylabel('x_1(n)');
n=-4:2;x=[1,-2,4,6,-5,8,10];n4=-10:10;
x41=2*exp(0.5*n4);x412=cos(0.1*pi*n4);
[x42,n42]=sigmult(x41,n4,x,n);
[x43,n43]=sigshift(x,n,-2);
[x44,n44]=sigmult(x412,n4,x43,n43);
[x4,n4]=sigadd(x42,n42,x44,n44);
subplot(2,1,2);stem(n4,x4);
xlabel('n');ylabel('x_2(n)');
set(gcf,'color','w');

%－－－－－实验1－－－－－
%－－－－－ex1_3－－－－－

clear all; close all; clc;
n=-5:30;b=[1,-1];a=[1];
x1=n.*(stepseq(0,-5,30)-stepseq(10,-5,30));
x2=(20-n).*[stepseq(105,30)-stepseq(20,-5,30)];
[x2,n2]=sigadd(x1,n,x2,n);
subplot(2,1,1);stem(n2,x2);
xlabel('n');ylabel('x_2(n)');
h=filter(b,a,x2);
subplot(2,1,2);stem(n,h);
xlabel('n');ylabel('h(n)');
set(gcf,'color','w');

%－－－－－实验2－－－－－
%－－－－－ex2_1_1－－－－－

clear all; close all; clc;
n=-5:30;b=[1,-1];a=[1];
x1=n.*(stepseq(0,-5,30)-stepseq(10,-5,30));
x2=(20-n).*[stepseq(10,-5,30)-stepseq(20,-5,30)];
[x2,n2]=sigadd(x1,n,x2,n);
subplot(2,1,1);stem(n2,x2);
xlabel('n');ylabel('x_2(n)');
h=filter(b,a,x2);
subplot(2,1,2);stem(n,h);
xlabel('n');ylabel('h(n)');
set(gcf,'color','w');

clear all; close all; clc;
t=0:0.001:1;xa=cos(20*pi*t);
Ts=0.01;N1=round(1/Ts);n1=0:N1;x1=cos(20*pi*n1*Ts);
subplot(3,1,1);plot(t,xa,n1*Ts,x1,'o');
ylabel('x_1(n)');title('Sampling of x_a(t) using Ts=0.01');
Ts=0.05;N2=round(1/Ts);n2=0:N2;x2=cos(20*pi*n2*Ts);
subplot(3,1,2);plot(t,xa,n2*Ts,x2,'o');
ylabel('x_2(n)');title('Sampling of x_a(t) using Ts=0.05');
Ts=0.1;N3=round(1/Ts);n3=0:N3;x3=cos(20*pi*n3*Ts);
subplot(3,1,3);plot(t,xa,n3*Ts,x3,'o');
ylabel('x_3(n)');title('Sampling of x_a(t) using Ts=0.1');
set(gcf,'color','w');

%－－－－－实验2－－－－－
%－－－－－ex2_1_2－－－－－

clear all; close all; clc;
t=0:0.001:1;xa=cos(20*pi*t);
xa=cos(20*pi*t);
Ts=0.01;N1=round(1/Ts);n1=0:N1;x1=cos(20*pi*n1*Ts);
Ts=0.01;Fs=1/Ts;
xa1=x1*sinc(Fs*(ones(length(n1),1)*t-(n1*Ts)'*ones(1,length(t))));
subplot(3,1,1);plot(t,xa1);axis([0,1,-1.1,1.1]);
ylabel('x_a(t)');title('Reconstruction of x_a(t) when Ts=0.01');
Ts=0.05;Fs=1/Ts;
Ts=0.05;N2=round(1/Ts);n2=0:N2;x2=cos(20*pi*n2*Ts);
xa2=x2*sinc(Fs*(ones(length(n2),1)*t-(n2*Ts)'*ones(1,length(t))));
subplot(3,1,2);plot(t,xa2);axis([0,1,-1.1,1.1]);
ylabel('x_a(t)');title('Reconstruction of x_a(t) when Ts=0.05');
Ts=0.1;Fs=1/Ts;
Ts=0.1;N3=round(1/Ts)n3=0:N3;x3=cos(20*pi*n3*Ts);
xa3=x3*sinc(Fs*(ones(length(n3),1)*t-(n3*Ts)'*ones(1,length(t))));
subplot(3,1,3);plot(t,xa3);axis([0,1,-1.1,1.1]);
ylabel('x_a(t)');title('Reconstruction of x_a(t) when Ts=0.1');
set(gcf,'color','w');

%－－－－－实验2－－－－－
%－－－－－ex2_1_3－－－－－

clear all; close all; clc;
t=0:0.001:1;xa=cos(20*pi*t);
Ts=0.01;N1=round(1/Ts);n1=0:N1;x1=cos(20*pi*n1*Ts);
Ts=0.01;Fs=1/Ts;
xa1=spline(Ts*n1,x1,t);
subplot(3,1,1);plot(t,xa1);axis([0,1,-1.1,1.1]);
ylabel('x_a(t)');title('Reconstruction of x_a(t) when Ts=0.01');
Ts=0.05;N2=round(1/Ts);n2=0:N2;x2=cos(20*pi*n2*Ts);
Ts=0.05;Fs=1/Ts;
xa2=spline(Ts*n2,x2,t);
subplot(3,1,2);plot(t,xa2);axis([0,1,-1.1,1.1]);
ylabel('x_a(t)');title('Reconstruction of x_a(t) when Ts=0.05');
Ts=0.1;N3=round(1/Ts);n3=0:N3;x3=cos(20*pi*n3*Ts);
Ts=0.1;Fs=1/Ts;
xa3=spline(Ts*n3,x3,t);
subplot(3,1,3);plot(t,xa3);axis([0,1,-1.1,1.1])
ylabel('x_a(t)');title('Reconstruction of x_a(t) when Ts=0.1');
set(gcf,'color','w');

%－－－－－实验2－－－－－
%－－－－－ex2_2－－－－－

clear all; close all; clc;
nx=0:9;x1=0.8.^nx;
ny=0:2;x2=[1 0 1];
[ng,g]=SeqLConv(nx,x1,ny,x2);
subplot(221),
stem(nx,x1,'.');grid;axis([-0.5 9.5 -0.2 1.2]);%stem-Plot discrete sequence data
xlabel('n');ylabel('x1(n)');title('序列x1(n)');
subplot(222),
stem(ny,x2,'.');grid;axis([-0.5 3.5 -0.2 1.2]);
xlabel('n');ylabel('x2(n)');title('序列x2(n)');
subplot(212),
stem(ng,g,'.');grid;axis([-1 12 -0.5 2]);
xlabel('n');ylabel('g(n)');title('线性卷积序列g(n)');
set(gcf,'color','w');

%－－－－－实验2－－－－－
%－－－－－ex2_3－－－－－

clear all; close all; clc;
nx=0:4;x=[1 2 3 4 5];
ny=0:5;y=[-1 3 0 -2 -2 1];
[nt1,t1]=SeqCirConv(nx,x,ny,y,6);
[nt2,t2]=SeqCirConv(nx,x,ny,y,10);
[nt3,t3]=SeqCirConv(nx,x,ny,y,12);

subplot(231),
stem(nx,x,'.');grid;axis([-1 5 -0.2 5.2]);
xlabel('n');ylabel('x(n)');title('序列x(n)');
subplot(232),
stem(ny,y,'.');grid;axis([-1 8 -2.2 3.2]);
xlabel('n');ylabel('y(n)');title('序列y(n)');
subplot(233),
stem(nt1,t1,'.');grid;axis([-1 8 -15.2 10]);
xlabel('n');ylabel('t5(n)');title('6点循环卷积');
subplot(234),
stem(nt2,t2,'.');grid;axis([-1 10 -15.2 10]);
xlabel('n');ylabel('t9(n)');title('10点循环卷积');
subplot(235),
stem(nt3,t3,'.');grid;axis([-1 13 -15.2 10]);
xlabel('n');ylabel('t12(n)');title('12点循环卷积');

subplot(236),
[ng,g]=SeqLConv(nx,x,ny,y);
stem(ng,g,'.');grid;axis([-1 13 -15.2 10]);
xlabel('n');ylabel('g(n)');title('线性卷积序列g(n)');

set(gcf,'color','w');

%－－－－－实验3－－－－－
%－－－－－ex3_1－－－－－

clear all; close all; clc;
x=[1,3,5,3,1];nx=0:4;T=0.5; % 给定原始数据
N=length(x);D=2*pi/(N*T); % 求出序列长度及频率分辨率
k=floor((-(N-1)/2):((N-1)/2)); % 求对称于零频率的FFT位置向量
X=fftshift(fft(x,N)); % 求对称于零频率的FFT序列值
subplot(21,1),plot(k*D,abs(X),'o:') % 画幅频特性图
title('幅频特性图');

subplot(2,1,2),plot(k*D,angle(X),'o:') % 画相频特性图
title('相频特性图');

set(gcf,'color','w');

%－－－－－实验3－－－－－
%－－－－－ex3_2－－－－－

clear all;close all;clc;
Fs=80;%抽样频率
SignalNum = 512;%信号点数
n=0:3000;
x = cos(2*20*pi*n/Fs)+2*cos(2*21*pi*n/Fs);
x1=x(1:SignalNum);
X1=fft(x1,2048);
subplot(211),plot(0:SignalNum-1,x1);xlabel('n');ylabel('x(n)');title('时域波形');grid;
subplot(212),plot(abs(X1));xlabel('k');ylabel('|X(k)|');title('幅频特性');grid;
set(gcf,'color','w');

+++++++++++++++++++++++以下是调用的子函数++++++++++++++++++++++++++++
function [nt,t]=SeqCirConv(nx,x,ny,y,N)
nt=0:N-1;
x1=[x zeros(1,N-length(x))];
y1=[y zeros(1,N-length(y))];
i=0:N-1;
y1=y1(mod(-i,N)+1);
Y=zeros(N,N);
for k=1:N
[nyy,yy]=SeqCShift(ny,y1,N,k-1);
Y(k,:)=yy;
end
t=x1*Y';

function [nxc,xc]=SeqCShift(nx,x,N,n0)
[nx1,x1]=SeqZQYT(nx,x,N,2);
if n0>N
n0=mod(n0,N);
end
if n0<-N
n0=n0+N;
end
nxc=min(nx):min(nx)+N-1;
if n0>0
xc=x1(N-n0+1:2*N-n0);
else
xc=x1(abs(n0)+1:abs(n0)+N);
end

function [ng,g]=SeqLConv(nx,x,ny,y)
%%% 函数说明 %%%%
%%% x和y代表参与卷积运算的序列，对应的位置向量用nx和ny表示 %%%%
%%% g代表线性卷积序列，对应的位置向量用ng表示 %%%%
ngmin(nx)+min(ny):max(nx)+max(ny);
g=conv(x,y);

function [ny,y]=SeqZQYT(nx,x,L,k)
%%% 函数功能 %%%%
%%% 将序列x(n)以L为周期进行周期延拓，自x(n)的起点始给出k个周期 %%%%
N=length(nx);
ny=min(nx):k*L+min(nx)-1;y=zeros(1,k*L);
for i=0:k*L-1
if L==N
y(i+1)=x(mod(i,L)+1);
end
if L>N
x1=[x zeros(1,L-N)];
y(i+1)=x1(mod(i,L)+1);
end
if (LN/2)
x2=[x(1:N-L)+x(L+1:N) x(N-L+1:L)];
y(i+1)=x2(mod(i,L)+1);
end
if L<=N/2
if mod(N,2)==1
xb=[x zeros(1,fix(N/2)-1)];
else
xb=[x zeros(1,N/2)];
end
x3=xb(1:L);
if mod(N,L)==0
for t=1:N/L-1
x3=x3+xb(1+t*L:(t+1)*L);
end
else
for t=1:fix(N/L)
x3=x3+xb(1+t*L:(t+1)*L);
end
end
y(i+1)=x3(mod(i,L)+1);
end
end

function [y,n]=sigadd(x1,n1,x2,n2)
% implements y(n)=x1(n)+x2(n)
% ——————————————————
% [y,n]=sigadd(x1,n1,x2,n2)
% y=sum sequence over n, which includes n1 and n2
% x1=first sequence over n1
% x2=second sequence over n2 (n2 can be different from n1)
%
n = min(min(n1),min(n2)):max(max(n1),max(n2)); % duration of y(n)
y1 = zeros(1,length(n));y2 = y1; % initialization
y1(find((n>=min(n1))&(n<=max(n1))==1))=x1; % x1 with duration of y
y2(find((n>=min(n2))&(n<=max(n2))==1))=x2; % x2 with duration of y
y = y1+y2; % sequence addition

function [y,n]=sigmult(x1,n1,x2,n2)
% implements y(n)=x1(n)*x2(n)
% ——————————————————
% [y,n]= sigmult(x1,n1,x2,n2)
% y=product sequence over n, which includes n1 and n2
% x1=first sequence over n1
% x2=second sequence over n2 (n2 can be different from n1)
%
n=min(min(n1),min(n2)):max(max(n1),max(n2)); % duration of y(n)
y1=zeros(1,length(n));y2=y1; % initialization
y1(find((n>=min(n1))&(n<=max(n1))==1))=x1; % x1 with duration of y
y2(find((n>=min(n2))&(n<=max(n2))==1))=x2; % x2 with duration of y
y=y1.*y2;

function [y,n]=sigshift(x,m,n0)
% implements yn)=x(n-n0)
% ——————————————————
% [y,n]=sigshift(x,m,n0)
%
n=m+n0;y=x;

function [x,n]=stepseq(n0,n1,n2)
% Generates x(n)=u(n-n0);n1<=n<=n2
% ——————————————————
% [x,n]=stepseq(n0,n1,n2)
%
n=[n1:n2];x=[(n-n0)>=0];

%---------实验四-------
%---------Ex4_1--------

%---------脉冲响应不变法-----
wp=0.2*pi; ws=0.3*pi;
Fs=4000; T=1/Fs
Wp=wp/T; Ws=ws/T;
alpha_p=2; alpha_s=40;
[c,d]=aft_butt(Wp,Ws,alpha_p,alpha_s);
[b,a]=impinvar(c,d,Fs);
[H,W]=freqz(b,a);
%---------画出幅度响应曲线-------------
subplot(2,1,1);plot(W*Fs/2/pi,abs(H));
xlabel('频率/Hz');ylabel('幅值');grid;
%---------画出相位响应曲线-------------
subplot(2,1,2);plot(W*Fs/2/pi,angle(H));
xlabel('频率/Hz');ylabel('相位');grid;

%---------实验四-------
%---------Ex4_2--------

%---------双线性变换法----------------
wp=0.2*pi; ws=0.3*pi;
Fs=4000; T=1/Fs
Wp=(2/T)*tan(wp/2);
Ws=(2/T)*tan(ws/2);
alpha_p=2; alpha_s=40;
[c,d]=aft_butt(Wp,Ws,alpha_p,alpha_s);
[b,a]=bilinear(c,d,Fs);
[H,W]=freqz(b,a);
%---------画出幅度响应曲线-------------
subplot(2,1,1);plot(W*Fs/2/pi,abs(H));
xlabel('频率/Hz');ylabel('幅值');grid;
%---------画出相位响应曲线-------------
subplot(2,1,2);plot(W*Fs/2/pi,angle(H));
xlabel('频率/Hz');ylabel('相位');grid;

%设计Butterworth低通滤波器-------------
function[b,a]=aft_butt(Wp,Ws,alpha_p,alpha_s)
N=ceil(log10((10^(alpha_p/10)-1)/(10^(alpha_s/10)-1))/(2*log10(Wp/Ws)));
fprintf('\n Butterworth滤波器阶数=%2.0f\n',N);
Omegac=Wp/((10^(alpha_p/10)-1)^(1/(2*N)));
[b,a]=uw_buttap(N,Omegac);

function [b,a]=uw_buttap(N,Omegac)
%未归一化的Butterworth模拟滤波器设计函数

[z,p,k]=buttap(N);
p=p*Omegac;
k=k*Omegac^N;
B=real(poly(z));
b0=k;
b=k*B;
a=real(poly(p));

%----------------实验五----------------------------------
%----------------ex5_1-----------------------------------
N=25;wc=1;
n=[0:1:N-1];
hd=ideal_lp(wc,N); %理想低通滤波器
%----------------使用矩形窗-------------------------------
w_bl1=(rectwin(N))'; %选择矩形窗
h1=hd.*w_bl1; %单位脉冲响应
[db1,mag1,pha1,grd1,w1,H1]=freqz_m(h1,[1]); %幅度响应
%----------------使用汉宁窗------------------------
w_bl2=(hann(N))'; %选择汉宁窗
h2=hd.*w_bl2; %单位脉冲响应
[db2,mag2,pha2,grd2,w2,H2]=freqz_m(h2,[1]); %幅度响应
%----------------画曲线-----------------------------------
subplot(2,1,1); plot(w1,mag1);
title('使用矩形窗设计FIR滤波器'); xlabel('频率'); ylabel('幅值');grid;
subplot(2,1,2); plot(w2,mag2);
title('使用汉宁窗设计FIR滤波器'); xlabel('频率'); ylabel('幅值');grid;

function hd=ideal_lp(wc,N)
alpha=(N-1)/2;
n=[0:1:N-1];
m=n-alpha+eps;
hd=sin(wc*m)./(pi*m);

function [db,mag,pha,grd,w,H]=freqz_m(b,a)
[H,w]=freqz(b,a,1000,'whole');
H=(H(1:1:501))';w=(w(1:501))';
mag=abs(H);
db=(mag+eps)/max(mag);
pha=angle(H);
grd=grpdelay(b,a,w)