西安电子科技大学数字信号处理大作业

这是西安电子科技大学计算机学院数字信号处理的大作业

 

coolwx当时选这门课的时候自以为信号学的很好，所以作为别的方向也“自信满满”选了这门数字信号处理，然而这门课真的硬核而且无聊

 

说实话，coolwx很佩服电子信息他们的课程，我觉得XDU中最简单的课程几乎都在CS系，而EE系的课程真的很难。

 

DSP作为一门硬核课程，他的大作业其实coolwx用了一个上午就搞定了，其实真的不是很难，只要你理解了书本（然而这本书不是很友好）

 

起初coolwx还想入门语音信号这种东西，然而经过一学期DSP的折磨，coolwx现在动摇了，转而去研究graphics

 

本大作业主要的思想如下：

利用傅里叶级数展开模拟信号

采样---数字化

进行DFT

设计低通高通带通滤波器

分别是IIR和FIR

具体的作业要求coolwx做完就删掉了，所以找不到了，不过后来的同学们，这门课肯定是祖传实验，看懂coolwx的代码应该还是很容易的。

coolwx之所以开源这份代码就是要给嵌入式后来的同学们留一个参考（因为coolwx其实是参考书本上写的，没有在网上找到xdu的数字信号处理大作业（很简单，其实coolwx也想找的，然而实在找不到））

 

所以开源以下代码留作纪念，留给后世一个参考。

t=0:0.1:30;%用于模拟域画图
TWOPI = 2*pi;%2pi用于傅里叶级数展开
w1 = TWOPI*1;%1
w2 = TWOPI*2;%2
w3 = TWOPI*3;%3
w4 = TWOPI*4;%1
w5 = TWOPI*5;%1
w6 = TWOPI*6;%2
w7 = TWOPI*7;%3
w8 = TWOPI*8;%3
w9 = TWOPI*9;%1
yt = 1*cos(w1*t)+2*sin(w2*t)+3*cos(w3*t+pi/2)+1*sin(w4*t+1/4*pi)+1*cos(w5*t)+2*cos(w6*t)+3*cos(w7*t)+3*cos(w8*t)+1*cos(w9*t); %显然，最高的f（频率）为9Hz,周期为1
figure(1);
plot(t,yt);
title("模拟域中的信号时域波形");
xlabel("t");
ylabel("y(t)");
%采样，采样角频率为4（小于2*f）
n1 = 0:1/4:1;%此时采样取一个周期，采样频率为4
y1n = 1*cos(w1*n1)+2*sin(w2*n1)+3*cos(w3*n1+pi/2)+1*sin(w4*n1+1/4*pi)+1*cos(w5*n1)+2*cos(w6*n1)+3*cos(w7*n1)+3*cos(w8*n1)+1*cos(w9*n1);
n1_number = 0:1:4;
figure(2);
subplot(2,1,1);
stem(n1_number,y1n);
title("采样信号，采样频率为4，采样一个周期");
xlabel("n");
ylabel("y1(n)");
%采样，采样角频率为20（大于2*f)
n2 = 0:1/20:1;%仍然是一个周期，也就是0到1的范围内
n2_number = 0:1:20;
y2n = 1*cos(w1*n2)+2*sin(w2*n2)+3*cos(w3*n2+pi/2)+1*sin(w4*n2+1/4*pi)+1*cos(w5*n2)+2*cos(w6*n2)+3*cos(w7*n2)+3*cos(w8*n2)+1*cos(w9*n2);
subplot(2,1,2);
stem(n2_number,y2n);
title("采样信号，采样频率为20，采样一个周期");
xlabel("n");
ylabel("y2(n)");
%下面对y1取4点fft，并画出其幅度和相位谱线
f1n = fft(y1n,4);%作fft
A1n = abs(f1n);%幅度谱
figure(3);
subplot(2,1,1);
n_number_draw = 0:3;
stem(n_number_draw,A1n);%画fft的幅频响应
title("对y1n做4点fft,画出幅频响应");
xlabel("n");
ylabel("A(n)幅频响应");
subplot(2,1,2);
fain = angle(f1n);
stem(n_number_draw,fain);%画fft的幅频响应
title("对y1n做4点fft,画出相频响应");
xlabel("n");
ylabel("fi(n)相频响应");
figure(4);
subplot(2,1,1);
stem(n_number_draw/4*2,A1n);
title("对y1n做4点fft,画出幅频响应,归一化到2pi中");
xlabel("n/pi");
ylabel("A(n)幅频响应之后归一化到2pi之中");
subplot(2,1,2);
stem(n_number_draw/4*2,fain);
title("对y1n做4点fft,画出相频响应，归一化到2pi之中");
xlabel("n/pi");
ylabel("fi(n)相频响应之后归一化到2pi之中");
f2n = fft(y2n,20);%作fft
A2n = abs(f2n);%幅度谱
figure(5);
subplot(2,1,1);
n_number_draw = 0:19;
stem(n_number_draw,A2n);%画fft的幅频响应
title("对y2n做20点fft,画出幅频响应");
xlabel("n");
ylabel("A(n)幅频响应");
subplot(2,1,2);
fain2 = angle(f2n);
stem(n_number_draw,fain2);%画fft的幅频响应
title("对y2n做20点fft,画出相频响应");
xlabel("n");
ylabel("fi(n)相频响应");
figure(6);
subplot(2,1,1);
stem(n_number_draw/20*2,A2n);%画fft的幅频响应
title("对y2n做20点fft,画出幅频响应，归一化到2pi之中");
xlabel("n/pi");
ylabel("fi(n)幅频响应");
subplot(2,1,2);
stem(n_number_draw/20*2,fain2);%画fft的幅频响应
title("对y2n做20点fft,画出相频响应，归一化到2pi之中");
xlabel("n/pi");
ylabel("fi(n)相频响应");

%下面进行IIR滤波器设计 （巴特沃斯滤波器，低通/高通/带通）
IIR_wp=0.3;%归一化通带边界频率 0.3*pi
IIR_ws=0.6;%归一化阻带边界频率 0.6*pi
Rp = 2;%通带最大衰减2db
As = 60;%阻带最小衰减60db
[N,wc]=buttord(IIR_wp,IIR_ws,Rp,As);
[B,A]=butter(N,wc);%算出巴特沃斯低通滤波器的B和A
%下面进行滤波
figure(7);
[H,w]=freqz(B,A,[0:2*pi/400:2*pi]);%先画幅频响应
Hf=abs(H); 
subplot(2,1,1);
plot(w/(2*pi)*2,(Hf));
xlabel("数字角频率w/pi");
ylabel("幅频响应");
subplot(2,1,2);
plot(w/(2*pi)*2,20*log10(Hf));
xlabel("数字角频率w/pi");
ylabel("损耗函数");

y = filter(B,A,y2n);
f3n = fft(y,20);
figure(8);
subplot(2,1,1);
stem(0:1:20,y);
title("时域的结果");
subplot(2,1,2);
title("做fft之后的频域的结果");
xlabel("n");
ylabel("幅频特性");
stem(0:19,abs(f3n));

%下面做高通滤波器
IIR_wp=0.6;%归一化通带边界频率 0.6*pi
IIR_ws=0.3;%归一化阻带边界频率 0.3*pi
Rp = 2;%通带最大衰减2db
As = 60;%阻带最小衰减60db
[N,wc]=buttord(IIR_wp,IIR_ws,Rp,As);
[B,A]=butter(N,wc,"high");%算出巴特沃斯低通滤波器的B和A
%下面进行滤波
figure(9);
[H,w]=freqz(B,A,[0:2*pi/400:2*pi]);%先画幅频响应
Hf=abs(H); 
subplot(2,1,1);
plot(w/(2*pi)*2,(Hf));
xlabel("数字角频率w/pi");
ylabel("幅频响应");
subplot(2,1,2);
plot(w/(2*pi)*2,20*log10(Hf));
xlabel("数字角频率w/pi");
ylabel("损耗函数");

y = filter(B,A,y2n);
f3n = fft(y,20);
figure(10);
subplot(2,1,1);
stem(0:1:20,y);
title("时域的结果");
subplot(2,1,2);
title("做fft之后的频域的结果");
xlabel("n");
ylabel("幅频特性");
stem(0:19,abs(f3n));

IIR_wp=[0.3,0.8];%归一化通带边界频率 0.3*pi,0.8pi
IIR_ws=[0.4,0.7];%归一化阻带边界频率 0.6*pi,0.7pi
Rp = 2;%通带最大衰减2db
As = 60;%阻带最小衰减60db
[N,wc]=buttord(IIR_wp,IIR_ws,Rp,As);
[B,A]=butter(N,wc);%算出巴特沃斯低通滤波器的B和A
%下面进行滤波
figure(11);
[H,w]=freqz(B,A,[0:2*pi/400:2*pi]);%先画幅频响应
Hf=abs(H); 
subplot(2,1,1);
plot(w/(2*pi)*2,(Hf));
xlabel("数字角频率w/pi");
ylabel("幅频响应");
subplot(2,1,2);
plot(w/(2*pi)*2,20*log10(Hf));
xlabel("数字角频率w/pi");
ylabel("损耗函数");

y = filter(B,A,y2n);
f3n = fft(y,20);
figure(12);
subplot(2,1,1);
stem(0:1:20,y);
title("时域的结果");
subplot(2,1,2);
stem(0:19,abs(f3n));
xlabel("n");
title("做fft之后的频域的结果");
ylabel("幅频特性");


%下面利用窗函数法设计FIR滤波器
wp = pi/2;
ws = pi/4;
Bt = wp-ws;
N0 = ceil(6.2*pi/Bt);%找出h(n)所需要的长度
N = N0+mod(N0+1,2);%由于是高通滤波器，所以必须保证hn的长度是奇数
wc = (wp+ws)/2/pi;
hn = fir1(N-1,wc,'high',hanning(N));
figure(13);
[H,w]=freqz(hn,1,[0:2*pi/400:2*pi]);%先画幅频响应
Hf=abs(H); 
subplot(3,1,1);
stem(0:1:N-1,hn);
title("h(n)的波形");
subplot(3,1,2);
plot(w/(2*pi)*2,20*log10(Hf));
xlabel("数字角频率w/pi");
ylabel("损耗函数");
title("损耗函数曲线");
subplot(3,1,3);
plot(w/(2*pi)*2,(Hf));
xlabel("数字角频率w/pi");
ylabel("幅频特性");
title("曲线");

figure(14);
y = filter(hn,1,y2n);
f3n= fft(y,20);

subplot(2,1,1);
stem(0:1:20,y);
title("时域的结果");
subplot(2,1,2);
stem(0:19,abs(f3n));
xlabel("n");
title("做fft之后的频域的结果");
ylabel("幅频特性");

%%下面设计低通滤波器
wp = pi/2;
ws = pi/4;
Bt = wp-ws;
N0 = ceil(6.2*pi/Bt);%找出h(n)所需要的长度
N = N0+mod(N0+1,2);
wc = (wp+ws)/2/pi;
hn = fir1(N-1,wc,hanning(N));
figure(15);
[H,w]=freqz(hn,1,[0:2*pi/400:2*pi]);%先画幅频响应
Hf=abs(H); 
subplot(3,1,1);
stem(0:1:N-1,hn);
title("h(n)的波形");
subplot(3,1,2);
plot(w/(2*pi)*2,20*log10(Hf));
xlabel("数字角频率w/pi");
ylabel("损耗函数");
title("损耗函数曲线");
subplot(3,1,3);
plot(w/(2*pi)*2,(Hf));
xlabel("数字角频率w/pi");
ylabel("幅频特性");
title("曲线");

figure(16);
y = filter(hn,1,y2n);
f3n= fft(y,20);

subplot(2,1,1);
stem(0:1:20,y);
title("时域的结果");subplot(2,1,2);stem(0:19,abs(f3n));xlabel("n");
title("做fft之后的频域的结果");
ylabel("幅频特性");