马上就要DSP考试了，DSP实验下周就结束了，课程设计也该提上日程了，啊哈 ，今天失踪许久的小六又回来了，带着新鲜出炉的低通滤波器来了。
首先这个功能大概就是根据音乐的幅频特性自己确定滤波器，应该可以降噪，然后声音变低沉变圆润，话不多说，来感受感受：
这是龙哥2018年留下的一段音频，也是输入音频：
这是输出音频：
确实符合我上面说的那些，然后我们一起来复习复习FIR滤波器的设计步骤：
(1)给定希望逼近的频率响应Hd(jw)，若所给指标为边界角频率和通带、阻带衰减，可选择理想滤波器作逼近函数。比如我们这里想设计一个低通滤波器，就可以使用理想低通滤波器的频率响应Hd(jw)，然后由此求单位冲激响应。(2)求单位冲激响应hd(n)，即作Hd(jw)的逆序列傅里叶变换(IDTFT)(3)根据阻带衰减指标，选择窗函数的形状，这个查表就好了，然后根据过渡带宽求得阶数N，同时得保证N的奇偶性。(4)求FIR滤波器的单位冲激响应h(n)=hd(n)*w(n)(5)计算FIR滤波器的频率响应H(jw)
然后到了代码解析时间，纯手工解析：
clc,clf,clear;%例行清理%% 音频获取[m,fs]=audioread('testaudio.mp3');%音频输入N=length(m);%获取长度N2=floor(N/2);%地板取整这里audioread函数是音频读取函数，里面写文件名，输出m是信号序列，fs是采样频率，我的这个文件比较小，fs是8000Hz，如果用手机或电脑录音的话可能是48000Hz，对matlab来说太大了，不过我觉得你可以再采样，这样刚好满足了课程设计要求。length函数是求长度。floor函数这里我愿称之为地板函数，取不大于本数的整数。
%% 傅里叶分析m_fft=fft(m);%傅里叶变换m_fft_a=abs(m_fft);%取模db=20.*log10(m_fft_a(1:N2));%换算分贝n=((0:N2-1)./N.*2)';%横轴频率换算到[0，1]subplot(311);plot(n,db);%画图title('输入幅频特性');xlabel('w/\pi');ylabel('20lg|H(jw)|(dB)');axis([0 1 -150 70]);这里就是实验里也常用的结构，基本都这个套路，腻了腻了。
%% 滤波器设计logicp=(db>40);%大于40分贝为1ffp=logicp.*((0:N2-1)')./N;%取出频率wp=max(ffp)*2*pi;%找边界logics=(db>15);%大于15分贝为1ffs=logics.*((0:N2-1)')./N;ws=max(ffs)*2*pi;这一部分是最具个人特色的一部分，也是这个设计的一个小亮点
(虽然可能有点偏离课程设计要求)这里先产生了一个逻辑序列，再作乘积找到最大值存在的位置，即边界，比如wp就是最后一个幅度超过40分贝的角频率，ws是最后一个幅度超过15分贝的角频率。这里40和15可以自己设定，它会根据信号本身频谱产生独特的设计要求，对我这种选择困难患者最好了。
dw=ws-wp;%过渡带宽NO=ceil(6.6*pi/dw);%确定阶数N，并取整NO=NO+mod(NO+1,2);%取奇数no=[0:NO-1];wh=(hamming(NO))';%获得窗wc=(wp+ws)/2;%获得截止频率a=(NO-1)/2;%获得ahd=sin(wc*(no-a+eps))./pi./(no-a+eps);%获得理想低通滤波器单位冲激响应h=hd.*wh;%获得数字滤波器单位冲激响应[H,w]=freqz(h,[1],N,'whole');%获得其频率响应A=abs(H);%取模DB=20.*log10((A+eps)/max(A));%换算分贝subplot(312);plot(w/pi,DB);%画图title('系统幅频特性');xlabel('w/\pi');ylabel('20lg|H(jw)|(dB)');axis([0 1 -150 70]);这里基本就是书后的了，这里自作主张选择了海明窗，没有什么别的想法，就是懒。不过这里要注意NO，结束之后你可以在旁边的工作区看到NO的值，可能会很大，太大了你可以调宽一下过渡带宽，即上面的“40”和“15”两个分贝。
%% 滤波M_fft=m_fft.*H;%获得输出频谱M_fft_a=abs(M_fft);%取模dbM=20.*log10(M_fft_a(1:N2));%换算分贝subplot(313);plot(n,dbM);%画图title('输出幅频特性');xlabel('w/\pi');ylabel('20lg|H(jw)|(dB)');axis([0 1 -150 70]);M=ifft(M_fft);%傅里叶反变换这里也没啥好说的，可能就是得到滤波器频率响应之后怎么办，我们知道这个滤波器频率响应就是系统函数，只要在频域作乘积就好了。
%% 输出音频audiowrite('testout.wav',M,fs);%输出这里是输出，之所以选择wav格式，因为m4a和mp4格式采样频率都要求44100和48000，如果是输出48000，输入8000最后得到的就是六倍速，此外它还支持ogg和flac格式，
当然这个程序仅供参考，因为它可能不是非常符合课程设计要求，其实这只是我给我自己的生日礼物，啊哈哈哈哈哈。再给大家看看图
可以看出高频部分被抑制的很厉害。
