数字信号处理实验参考程序

此程序仅供参考非标准答案。

实验一
Q1-4
clear, close all;
N = 40;
num = [0.9 -0.45 0.35 0.002];
den = [1 0.71 -0.46 -0.62];
y = impz(num,den,N);
% Plot the impulse response
stem(y);
xlabel('Time index n'); ylabel('Amplitude');
title('Impulse Response'); grid;
Q1-5
clear,close all,
n=0:39;
x=delta(n);
num = [0.9 -0.45 0.35 0.002];
den = [1 0.71 -0.46 -0.62];
y=filter(num,den,x);
stem(n,y);
xlabel('Time index n');ylabel('Amplitute');
title('Impulse Response');grid;
Q1-6
clear, close all;
N = 40;
num = [0.9 -0.45 0.35 0.002];
den = [1 0.7 -0.46 -0.62];
y = stepz(num,den,N);
% Plot the impulse response
stem(y);
xlabel('Time index n'); ylabel('Amplitude');
title('Impulse Response'); grid;
Q1-8
% Program_P1_10
clear, close all;
h = [3 2 1 -2 1 0 -4 0 3 3 2 1 -2 0 -4]; % impulse response x = [1 -2 3 -4 3 2 1 1 -2 3]; % input sequence
y = conv(h,x); % convolution
n = 0:23;
subplot(2,1,1);
stem(n,y, '. ');
title('Output Obtained by Convolution'); grid;
x1 = [x zeros(1,14)];
y1 = filter(h,1,x1);
subplot(2,1,2);
stem(n,y1, '. ');
xlabel('Time index n'); ylabel('Amplitude');
title('Output Generated by Filtering'); grid;
% Program_P1_4
% Stability test based on the sum of the absolute values of the impulse response samples clear,close all;
num = [1 -4 3]; den = [1 -1.7 1]; % Define a system
N = 200;
h = impz(num,den,N+1);
parsum = 0;
for k = 1:N+1;
parsum = parsum + abs(h(k));
if abs(h(k)) < 10^(-6), break, end
end
% Plot the impulse response
n = 0:N;
stem(n,h)
xlabel('Time index n'); ylabel('Amplitude');
% Print the value of abs(h(k))
disp('V alue =');disp(abs(h(k)));
实验二
% Program P2_3
% Evaluation of the DTFT
clear, close all;
% Compute the frequency samples of the DTFT
w = -4*pi:8*pi/511:4*pi;
num = [0.7 -0.5 0.3 1]; den = [1 0.3 -0.5 0.7];
h = freqz(num, den, w);
% Plot the DTFT
subplot(2,2,1)
plot(w/pi,real(h)); grid; title('Real part of H(e^{j\omega})')
xlabel('\omega /\pi'); ylabel('Amplitude');
subplot(2,2,2)
plot(w/pi,imag(h)); grid;title('Imaginary part of H(e^{j\omega})')
xlabel('\omega /\pi'); ylabel('Amplitude');
subplot(2,2,3)
plot(w/pi,abs(h)); grid;title('Magnitude Spectrum |H(e^{j\omega})|')
subplot(2,2,4)
plot(w/pi,angle(h)); grid;title('Phase Spectrum arg[H(e^{j\omega})]')
xlabel('\omega /\pi'); ylabel('Phase in radians');
% Program P2_4
% Evaluation of the DTFT
clear, close all;
% Compute the frequency samples of the DTFT
w = -4*pi:8*pi/511:4*pi;
num = [1 3 5 7 9 11 13 15 17]; den = [1 0 0 0 0 0 0 0 0];
h = freqz(num, den, w);
% Plot the DTFT
subplot(2,2,1)
plot(w/pi,real(h)); grid; title('Real part of H(e^{j\omega})')
xlabel('\omega /\pi'); ylabel('Amplitude');
subplot(2,2,2)
plot(w/pi,imag(h)); grid;title('Imaginary part of H(e^{j\omega})')
xlabel('\omega /\pi'); ylabel('Amplitude');
subplot(2,2,3)
plot(w/pi,abs(h)); grid;title('Magnitude Spectrum |H(e^{j\omega})|')
xlabel('\omega /\pi'); ylabel('Amplitude');
subplot(2,2,4)
plot(w/pi,angle(h)); grid;title('Phase Spectrum arg[H(e^{j\omega})]')
xlabel('\omega /\pi'); ylabel('Phase in radians');
2.5
% Program2.5
clear, close all;
w = -4*pi:0.02:4*pi;
n = 0:8;
x = (n+1);
X = freqz(x,1,w);
subplot(221)
stem(n,x,'.'), title('Sequence x[n]'), xlabel('Time index n')
subplot(222)
plot(w/pi,abs(X)), title('M agnitude of DTFT of x[n]'), axis([-4,4,0,max(abs(X))]) subplot(224)
plot(w/pi,angle(X)), title('Phase of DTFT of x[n]'), xlabel('Frequency*pi in ratians') axis([-4,4,min(angle(X)),max(angle(X))])
% Program_2_12
clear, close all;
L = input('L=');
k = 0:L-1;
x = input('x=');
X = fft(x,L);
y=ifft(X,L);
subplot(221)
stem(x,'.'), title('Sequence x[n]'), xlabel('Time index n'), axis([0,15,min(x),max(x)]) subplot(222)
stem(k,abs(X),'.'), title('M agnitude of DFT of x[n]'), axis([0,15,0,max(abs(X))]) subplot(224)
stem(k,angle(X),'.'), title('Phase of DFT of x[n]'), xlabel('Frequency index k')
axis([0,15,min(angle(X)),max(angle(X))])
subplot(223)
stem(x,'.'), title('Sequence y[n]'), xlabel('Time index n'), axis([0,15,min(x),max(x)])
% Program_2_15
clear,close all;
g=input('g=');
h=input('h=');
N=length(g);
x=g+j*h;
X=fft(x,N);
for k=0:N-1
G(k+1)=(X(k+1)+conj(X(mod(-k,N)+1)))/2,
H(k+1)=(X(k+1)-conj(X(mod(-k,N)+1)))./2j
end
subplot(221)
stem(g,'.'), title('Sequence g[n]'), xlabel('Time index n'), axis([0,N,0,max(g)]) subplot(222)
stem(h,'.'), title('Sequence h[n]'), xlabel('Time index n'), axis([0,N,0,max(h)]) subplot(223)
stem(abs(G),'.'), title('M agnitude of DFT of g[n]'), axis([0,N,0,max(abs(G))]) subplot(224)
stem(abs(H),'.'), title('M agnitude of DFT of h[n]'), axis([0,N,0,max(abs(H))])
% Program_2_17
clear,close all;
x=0:9;
M=input('M=');
y=circshift(x,M);
subplot(221)
stem(x,'.'), title('Sequence x[n]'), xlabel('Time index n'), axis([0,9,0,max(x)]) subplot(222)
stem(y,'.'), title('Sequence y[n]'), xlabel('Time index n'), axis([0,9,0,max(y)])
% Program_2_18
clear,close all;
g=[1 2 0 1];
h=[2 2 1 1];
y=circonv(g,h);
subplot(221)
stem(g,'.'), title('Sequence g[n]'), xlabel('Time index n'), axis([0,4,0,3])
subplot(222)
stem(h,'.'), title('Sequence h[n]'), xlabel('Time index n'), axis([0,4,0,3])
subplot(223)
stem(y,'.'), title('Sequence y[n]'), xlabel('Time index n'), axis([0,4,0,max(y)])
% Program_2_19
clear, close all;
g=input('g=');
h=input('h=');
M=length(h);
N=length(g);
ge=[g zeros(1,M-1)];
he=[h zeros(1,N-1)];
Ge=fft(ge,M+N-1);
He=fft(he,M+N-1);
y=ifft(Ge.*He,M+N-1);
subplot(221)
stem(g,'.'), title('Sequence g[n]'), xlabel('Time index n'), axis([0,M+N-1,0,max(g)]) subplot(222)
stem(h,'.'), title('Sequence h[n]'), xlabel('Time index n'), axis([0,M+N-1,0,max(h)]) subplot(223)
stem(y,'.'), title('Sequence y[n]'), xlabel('Time index n'), axis([0,M+N-1,0,max(y)])
% Program_2_20
clear,close all,
N=input('Type in the DFT length N=');
n=0:N-1;
x=input('Type in the sequence x(n)=');
k=0:N-1;
X=zeros(1,N);
for m=0:N-1;
for j=0:N-1;
X(m+1)=X(m+1)+x(j+1)*exp(-i*2*pi*m*j/N);
end
end
subplot(211)
stem(n,x,'.')
axis([0,N,1.2*min(x),1.2*max(x)])
title('The original sequence x[n]')
xlabel('Index n')
subplot(212)
stem(k,abs(X),'r.')
axis([0,N,-0.2*max(abs(X)),1.2*max(abs(X))]) title('The DFT of the seqence x[n]')
xlabel('Index k')
% Program_2_21
clear,close all;
N=input('Type in the DFT length N=');
M=log2(N);
n=0:N-1;
x=input('Type in the sequence x(n)=');
x1=zeros(1,N);
x1(1)=x(1);
x1(N)=x(N);
J=0;
for I=1:N-2;
k=1;c=0;
while k~=M+1&c~=1;
if J>=N/(2.^k);
J=J-N/(2.^k);c=0;
else J=J+N/(2.^k);c=1;
end
k=k+1;
end
x1(I+1)=x(J+1);
end
m = 0:N/2-1;
W = exp(-i*2*pi*m/N);
for r=1:M;
B=2^(r-1);
for J=0:B-1;
p=2^(M-r)*J;
for k=J+1:2^r:N-1;
q=x1(k);
x1(k)=q+x1(k+B)*W(p+1);
x1(k+B)=q-x1(k+B)*W(p+1);
end
end
end
k=0:N-1;
subplot(211)
stem(k,x,'.'),title('The sequence x[n]')
xlabel('Index k'),axis([0,max(k),1.1*min(x),1.1*max(x)])
subplot(212)
stem(k,abs(x1),'r.'),title('The DFT of the sequence x[n]')
xlabel('Index k'),axis([0,max(k),-0.2*max(abs(x1)),1.2*max(abs(x1))])
实验四
%Q4_1
clear,close all;
n=0:40;
wp = input('wp=');
ws = input('ws=');
Rp = input('Rp=');
Rs = input('Rs=');
[N,wn] = buttord(wp,ws,Rp,Rs,'s');
[b,a]= butter(N,wn,'s');
[H,w] = freqs(b,a);
subplot(221)
plot(w/pi,abs(H));grid on;title('Magnitude response')
subplot(222)
plot(w/pi,angle(H));grid on;title('Phase response')
subplot(223)
plot(w/(2*pi),20*log(abs(H)));
axis([0,6000,-150,10])
grid on
title('The gain of the butterworth filter')
xlabel('Frequency index w/2*pi');
subplot(224)
h = impz(b,a,12);
n = 0:11;
stem(n,h/max(abs(h)),'.');grid on
title('The impulse response--h=h/max(abs(h))')
%Q4_3
clear,close all;
wp = input('wp=');
ws = input('ws=');
Wp = tan(wp/2);
Ws = tan(ws/2);
Rp = input('Rp=');
Rs = input('Rs=');
[N,wn] = buttord(Wp,Ws,Rp,Rs,'s');
[b,a]= butter(N,wn,'s');
[b,a] = bilinear(b,a,0.5);
[H,w] = freqz(b,a);
subplot(221)
plot(w/pi,abs(H));grid on;title('Magnitude response')
subplot(222)
plot(w/pi,angle(H));grid on;title('Phase response')
subplot(223)
plot(w/(pi),20*log(abs(H)));
axis([0,1,-750,5])
grid on
title('The gain of the IIR filter')
xlabel('Frequency index w/pi')
subplot(224)
h = impz(b,a,50);
n = 0:49;
stem(n,h,'.');grid on
title('The impulse response')
%Q4_5
clear,close all;
x = [-4, -2, 0, -4, -6, -4, -2, -4, -6, -6, -4, -4, -6, -6, -2, 6, 12, 8, 0, -16, -38, -60, -84, -90, -66, -32, -4, -2, -4, 8, 12, 12, 10, 6, 6, 6, 4, 0, 0, 0, 0, 0, -2, -4, 0, 0, 0, -2, -2, 0, 0, -2, -2, -2, -2, 0];
dw = 0.01;
w = -pi:dw:pi;
L = length(x);
n = 0:L-1;
subplot(222)
stem(n,x,'.')
title('The original signal')
grid
wp = 0.2*pi;
ws = 0.3*pi;
Wp = tan(wp/2);
Ws = tan(ws/2);
Rp = 1;
Rs = 30;
[N,wn] = buttord(Wp,Ws,Rp,Rs,'s');
[b,a]= butter(N,wn,'s');
[b,a] = bilinear(b,a,0.5);
H = freqz(b,a,w);
subplot(221)
plot(w/(pi),20*log(abs(H)));
axis([0,0.5,-150,0])
grid on
title('The gain of the IIR filter')
xlabel('Frequency index w/pi')
subplot(223)
h = impz(b,a,L);
stem(n,h,'.')
title('The impulse response')
grid
xlabel('Time index n')
X = freqz(x,1,w);
Y = X.*H;
y = dw*Y*exp(j*w'*n)/(2*pi);
subplot(224)
stem(n,y,'.')
title('The output signal of the filter')
grid
xlabel('Time index n')
% Q4_6
clear,close all;
wp=input('Type in the passband edge frequency (in Radian) wp='); ws=input('T ype in the passband edge frequency (in Radian) ws='); deltaw = ws-wp;
wc=(wp+ws)/2;
L=ceil(8*pi/(ws-wp));
if L/2==ceil(L/2);
L=L+1;
else
L=L;
end
a=(L-1)/2;
n=0:L-1;
w=0.5*(1-cos(2*pi*n/(L-1))).*(u(n)-u(n-L));
h=(wc/pi)*sin(wc*(n-a+eps))./(wc*(n-a+eps)).*w;
N=2000;
H=fft(h,N);
H=H/max(abs(H));
W=0:2*pi/N:2*pi-2*pi/N;
subplot(222)
plot(W/pi,abs(H)),title('The magnitude response of the FIR filter') axis([0,1,0,max(abs(H))]),grid on;
subplot(223)
plot(W/pi,angle(H))
axis([0,1,-pi,pi]),grid on;
title('The phase response of the FIR filter'),xlabel('Frequency (*pi)')
subplot(224)
plot(W/pi,20*log10(abs(H)));
axis([0,1,-100,0]),grid on;
title('The frequency response in dB'),
xlabel('Frequency (*pi)')
subplot(221)
stem(n,h,'.'),
axis([0,L-1,min(h),max(h)]),grid on;
title('The impulse response of the FIR filter')
xlabel('Index n')
% Q4_8
clear,close all;
x = [-4, -2, 0, -4, -6, -4, -2, -4, -6, -6, -4, -4, -6, -6, -2, 6, 12, 8, 0, -16, -38, -60, -84, -90, -66, -32, -4, -2, -4, 8, 12, 12, 10, 6, 6, 6, 4, 0, 0, 0, 0, 0, -2, -4, 0, 0, 0, -2, -2, 0, 0, -2, -2, -2, -2, 0];
wp=input('Type in the passband edge frequency (in Radian) wp=');
ws=input('T ype in the passband edge frequency (in Radian) ws=');
deltaw = ws-wp;
wc=(wp+ws)/2;
L=ceil(8*pi/(ws-wp));
if L/2==ceil(L/2);
L=L+1;
else
L=L;
end
a=(L-1)/2;
n=0:L-1;
w=0.5*(1-cos(2*pi*n/(L-1))).*(u(n)-u(n-L));
h=(wc/pi)*sin(wc*(n-a+eps))./(wc*(n-a+eps)).*w;
y = conv(x,h);
subplot(211)
n = 0:length(x)-1;
stem(n,x,'.')
axis([0,100,-100,50])
grid
title('The original signal')
subplot(212)
n = 0:length(y)-1;
stem(n,y,'.')
grid
axis([0,100,-100,50])
title('The output signal of the filter')
xlabel('Time index n')。