数字信号处理上机作业

第一题

题目描述

以自己的学号建立序列$x[n]$,$x[n]={5,1,6,0,2,1,9,1,0,0,3,9},n={0,1,…,11}$。

  1. 画出$x[n]$的傅里叶变换$X(e^{jω})$的幅度和相位曲线;
  2. 画出$x[n]$的32点DFT $X[k]$,要求与(1)画在同一幅图上,验证频域取样关系;
  3. 画出$X[k]$经过32点IDFT得到的$x[n]$,验证DFT与IDFT的唯一性。

MATLAB代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
clear,clc;
N = 32; %点数
x = zeros(1,N);
x(1:12) = [5,1,6,0,2,1,9,1,0,0,3,9];
n = 0:N-1;
cnt = 1;
y = zeros(201);

%用fft和ifft函数求DFT和IDFT
f_y = fft(x,N);
f_x = 2*pi/N.*n;
n_x = ifft(f_y,N);

%求傅里叶变换
syms w;
y = sum(x.*exp(-1i*n*w));

%绘图
figure;
subplot(2,1,1)
fplot(abs(y),[0,2*pi])
hold on
stem(f_x,abs(f_y))

%设置坐标轴刻度与标题
xlabel('ω');
ylabel('幅值');
set(gca,'XTick',0:pi/2:2*pi)
set(gca,'xtickLabel',{'0','π/2','π','3π/2','2π'})
axis([0 pi*2 0 40])

subplot(2,1,2)
fplot(angle(y),[0,2*pi])

%设置坐标轴刻度与标题
xlabel('ω');
ylabel('相位');
set(gca,'XTick',0:pi/2:2*pi)
set(gca,'xtickLabel',{'0','π/2','π','3π/2','2π'})
set(gca,'YTick',-pi:pi/2:pi)
set(gca,'ytickLabel',{'-π','π/2','0','π/2','π'})
axis([0 pi*2 -pi pi])

figure;
stem(0:N-1,n_x); %32点IDFT
axis([0 31 0 10])

图像

第二题

题目描述

采用脉冲响应不变法设计一个巴特沃斯离散时间低通IIR滤波器,要求通带截止频率$ω_p=0.3π rad$,阻带截止频率$ω_s=0.4π rad$,通带最大衰减$α_p=1 dB$,阻带最小衰减$α_s=30 dB$。取$T_d=1$,给出滤波器阶数、连续时间和离散时间滤波器的系统函数,并画出连续时间和离散时间系统的对数幅度响应曲线和单位脉冲响应曲线。

阶数与系统函数

滤波器阶数$N=15$

连续时间滤波器系统函数为
$$
H_c(s) = \frac{\sum^{N}{k=0} b_k s^{N-k}} {\sum^{N}{k=0} a_k s^{N-k}},其中
\
a_k = [1,9.5497,45.5983,144.0774,336.2895,613.3351,\902.0572,1088.6924,1086.7497,897.2369,607.8823,\332.1113,141.7800,44.7113,9.3305,0.9735]
\
b_k = [0,0,0,0,0,0,0,0,\0,0,0,0,0,0,0,\0.973565591260273]
$$
离散时间滤波器系统函数
$$
H(z) = \frac{\sum^{N}{k=0} b_k z^{-k}} {\sum^{N}{k=0} a_k z^{-k}},其中
\
a_k = [1,-5.9705,17.9738,-35.5026,
\50.9018,-55.7017,47.8036,-32.6279,
\17.7990,-7.7395,2.6547,-0.7041,
\0.1395,-0.0194,0.0017,-7.122210^{-5}]
\
b_k=[-2.4429
10^{-8},1.180610^{-7},6.558710^{-6},0.0001,0.0009,
\0.0024,0.0027,0.0014,0.0003,4.0372e-05,1.726610^{-6},
\2.2229
10^{-8},2.824810^{-11},6.998710^{-13},0]
$$

MATLAB代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
clear,clc;
wp = 0.3*pi; ap = 1;
ws = 0.4*pi; as = 30;
Td = 1;
Wp = wp/Td;
Ws = ws/Td;
%求连续时间滤波器阶数和截止频率
[N,Wc] = buttord(Wp,Ws,ap,as,'s')
%连续时间滤波器系统函数
[Bs,As] = butter(N,Wc,'s');
%脉冲响应不变法得到离散滤波器系统函数
[Bz,Az] = impinvar(Bs,As,1/Td);
%绘图
width = 1.5;
[H,W] = freqs(Bs,As);
plot(W*Td/pi,20*(log10(abs(H))),'LineWidth',width);
hold on;
[H,w]=freqz(Bz,Az);
plot(w/pi,20*(log10(abs(H))),'LineWidth',width);
axis([0,1,-200,2]);
grid;
legend('离散','连续');
xlabel('ω(*π rad) Ω(*π rad/s)');
ylabel('幅值响应(dB)');
figure;
[h,n]=impz(Bz,Az);
stem(n,h,'LineWidth',width);
hold on;
g=tf(Bs,As);
[h1,t]=impulse(g,0:0.1:40);
plot(t,h1,'LineWidth',width);
legend('离散','连续');
axis([0,40,-0.3,0.5]);
xlabel('t(s)或n(样本)');
ylabel('单位脉冲响应');

图像