西安电子科技大学数字信号处理大作业
这是西安电子科技大学计算机学院数字信号处理的大作业
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("幅频特性");