準(zhǔn)備前的工作——一段收集來(lái)的代碼
一诞仓、輸入信號(hào)時(shí)域&頻域圖像（非聲音信號(hào)，而是函數(shù)）
f1=10;%第一個(gè)點(diǎn)頻信號(hào)分量頻率
f2=30;%第二個(gè)點(diǎn)頻信號(hào)分量頻率
f3=45;%第三個(gè)點(diǎn)頻信號(hào)分量頻率
fs=100;%采樣率
T=2;%間隔長(zhǎng)度
n=round(T*fs);%采樣點(diǎn)個(gè)數(shù)
t=linspace(0,T,n);
y=cos(2*pi*f1*t)+cos(2*pi*f2*t)+cos(2*pi*f3*t)+randn(size(t));
figure;
subplot(2,1,1);
plot(t,y);
title('輸入信號(hào)時(shí)域圖像');
xlabel('t/s');
ylabel('V');
fft_y=fftshift(fft(y));
f=linspace(-fs/2,fs/2,n);
subplot(2,1,2);
plot(f,abs(fft_y));
title('輸入信號(hào)頻域圖像');
xlabel('f/Hz');
ylabel('V');
https://blog.csdn.net/Qinlong_Stm32/article/details/84570524
文章介紹的簡(jiǎn)單易懂↑↑↑
——————————————————————————————
思考流程：
接下來(lái)，我需要對(duì)該信號(hào)進(jìn)行處理：
首先，如何加指定頻段加入固定幅值的噪聲偷溺？頻域圖像的橫坐標(biāo)又該如何確定跨释？
其次，如何設(shè)計(jì)IIR濾除指定頻段的噪聲吹截？該頻譜圖又該如何確定橫坐標(biāo)瘦陈？
最后，濾波前和濾波后的誤差曲線如何畫波俄？誤差均方值又該如何計(jì)算晨逝？
一、輸入聲音信號(hào)的時(shí)域圖像和頻域圖像
clear;
clc;
%輸入原始信號(hào)
[y,fs]=audioread('E:\son\dudu.wav'); %將聲音放于matlab中
info=audioinfo('E:\son\dudu.wav')
%sound(y,fs);
T=1/fs; %采樣時(shí)間
t=(0:length(y)-1)*T;%時(shí)間
f=(0:length(y)-1)*fs/length(y);
figure(1);
yz=y(:,1);%左聲道
subplot(2,1,1);
plot(t,yz);%輸入信號(hào)時(shí)域曲線
title('原始信號(hào)時(shí)域');
xlabel('時(shí)間');
ylabel('振幅');
subplot(2,1,2);
n=length(yz);%進(jìn)行變換的點(diǎn)數(shù)
y1=fft(yz,n); %對(duì)n點(diǎn)進(jìn)行傅里葉變換到頻域
F=fs/length(yz); %譜分辨率懦铺，頻譜間隔
plot(f,abs(y1));%左聲道頻譜圖
axis([0,1000,0,10000]);
title('原始信號(hào)頻譜');
xlabel('F(Hz)');
ylabel('H(jw)');
grid on
Note: 在plot頻譜圖的時(shí)候捉貌，我取的范圍是（-fs/2,fs/2-F）,為什么不取（0冬念，fs/2-F）呢趁窃？
因?yàn)槲医酉聛?lái)要在100,500,900處加噪聲，前者的取法會(huì)使中間頻率集中在0附近急前，這樣加噪聲會(huì)對(duì)聲音有影響醒陆；而后者會(huì)使頻率集中在2×10^4處，前面的部分幾乎振幅為0裆针，加噪聲后也沒(méi)有影響刨摩。如何選取坐標(biāo)范圍是根據(jù)需求或者要求來(lái)的。但需要注意的是据块，如果y需要選取0~1000hz的码邻，那么對(duì)應(yīng)的f也應(yīng)該改變。
二另假、加100像屋、500、900Hz的噪聲
加上噪聲以后會(huì)出現(xiàn)边篮，這樣的錯(cuò)誤顯示己莺。原因是因?yàn)橄嗉拥臅r(shí)候維度搞錯(cuò)了奏甫。yz是1×n的列向量，x是n×1的行向量凌受，無(wú)法相加阵子。我一開(kāi)始并沒(méi)有想到這上面，這是一種基礎(chǔ)的語(yǔ)法錯(cuò)誤胜蛉，我需要的是經(jīng)驗(yàn)挠进。
#加噪聲
noise1=sin(2*pi*100*t);
noise2=sin(2*pi*500*t);
noise3=sin(2*pi*900*t);
x1=noise1+noise2+noise3;
x1=yz+x1';
figure(2);
subplot(2,1,1);
plot(t,x1);
title('加噪信號(hào)時(shí)域');
xlabel('時(shí)間');
ylabel('振幅');
subplot(2,1,2);
nx=length(x1);
x=fft(x1,nx);
plot(f,abs(x));%左聲道頻譜圖
axis([0,1000,0,10000]);
title('加噪信號(hào)頻譜');
xlabel('F(Hz)');
ylabel('H(jw)');
三、設(shè)計(jì)IIR濾波器
數(shù)字濾波器的設(shè)計(jì)方法誊册，三個(gè)步驟：
1领突、給出所需要的濾波器的技術(shù)指標(biāo)；
2案怯、設(shè)計(jì)一個(gè)使其逼近所需要的技術(shù)指標(biāo)君旦；
實(shí)現(xiàn)所設(shè)計(jì)的。
3嘲碱、IIR濾波器的設(shè)計(jì)主要借助模擬濾波器z做轉(zhuǎn)換來(lái)進(jìn)行金砍；FIR的設(shè)計(jì)主要建立在對(duì)理想濾波器頻率特性作某種近似的基礎(chǔ)上，有窗函數(shù)法麦锯、頻率抽樣法等恕稠；
線性相位濾波器通常采用FIR型，其單位脈沖響滿足一定條件時(shí)扶欣，其相位特性是嚴(yán)格線性的谱俭。
一、IIR濾波器的設(shè)計(jì)步驟如下：
1)宵蛀、按一定規(guī)則把給定的數(shù)字濾波器技術(shù)指標(biāo)轉(zhuǎn)換為模擬低通濾波器的技術(shù)指標(biāo)昆著；
2)、根據(jù)轉(zhuǎn)換后的技術(shù)指標(biāo)設(shè)計(jì)模擬低通濾波器系統(tǒng)函數(shù)H(s)术陶；
3)凑懂、再按一定規(guī)則將轉(zhuǎn)換成；若所設(shè)計(jì)的數(shù)字濾波器是低通的梧宫，那么完成整個(gè)設(shè)計(jì)工作接谨；
4)、將高通塘匣，帶通或帶阻數(shù)字濾波器的技術(shù)指標(biāo)轉(zhuǎn)化為低通模擬濾波器的技術(shù)指標(biāo)脓豪，然后從第(2)步開(kāi)始設(shè)計(jì)低通濾波器H(s)，再將轉(zhuǎn)換為所需的H(z)忌卤。
我現(xiàn)在需要確定基本參數(shù)扫夜，根據(jù)該參數(shù)設(shè)計(jì)濾波器。如何確定參數(shù)？——fdatool笤闯；確定參數(shù)以后如何寫代碼堕阔？
1.參數(shù)可以通過(guò)fdatool確定，其實(shí)我覺(jué)得直接自定義一個(gè)即可颗味。衰減一個(gè)是0.1超陆，一個(gè)是60dB。濾波器貌似都是這個(gè)標(biāo)準(zhǔn)(不知道為什么)
2.如何將濾波器圖形轉(zhuǎn)化成代碼語(yǔ)言——可以參考書本浦马，很簡(jiǎn)單的定義
3.我現(xiàn)在的問(wèn)題是时呀，如何將濾波器加到聲音信號(hào)里實(shí)現(xiàn)濾波。
######我有點(diǎn)累了晶默。退唠。。荤胁。∈赫現(xiàn)在搞出來(lái)一個(gè)半成品（因?yàn)椴惶_定這個(gè)圖像和聲音對(duì)不對(duì)總感覺(jué)很怪+對(duì)于代碼有一些知識(shí)盲點(diǎn)未完全吃透）仅政，這些知識(shí)點(diǎn)需要攻破，比如每個(gè)濾波器的公式這么多選哪個(gè)才好盆驹？圖像的輸出的橫坐標(biāo)怎么確定圆丹？
######算了，先寫FIR的代碼吧躯喇。IIR晚點(diǎn)再繼續(xù)寫辫封。
四、待濾波信號(hào)加濾波器
%%%%%%%###帶阻濾波器100Hz####%%%%%
x2=filter(ww100,x1);
figure(3);
subplot(2,1,1);
plot(t,x2);
title('100HZ被濾掉后的時(shí)域');
xlabel('時(shí)間');
ylabel('振幅');
subplot(2,1,2);
plot(f,abs(fft(x2)));
axis([0 1000 0 10000]);
title('100Hz被濾后的頻譜');
xlabel('F(Hz)');
ylabel('H(jw)');
%%%%%%%#####帶阻濾波器500Hz#####%%%%%%%%%
x3=filter(ww500,x2);
figure(4);
subplot(2,1,1);
plot(t,x3);
title('100HZ廉丽、500HZ被濾掉后的時(shí)域');
xlabel('時(shí)間');
ylabel('振幅');
subplot(2,1,2);
plot(f,abs(fft(x3)));
axis([0 1000 0 10000]);
title('100Hz倦微、500Hz被濾后的頻譜');
xlabel('F(Hz)');
ylabel('H(jw)');
%%%%%%%%%########帶阻濾波器900Hz######%%%%%%%%%
x4=filter(ww900,x3);
figure(5);
subplot(2,1,1);
plot(t,x4);
title('全被濾掉后的時(shí)域');
xlabel('時(shí)間');
ylabel('振幅');
subplot(2,1,2);
plot(f,abs(fft(x4)));
axis([0 1000 0 10000]);
title('全部被濾后的頻譜');
xlabel('F(Hz)');
ylabel('H(jw)');
sound(x4,fs);
五、濾波前和濾波后的誤差曲線
yend=yz-x4; %誤差
figure(6);
subplot(2,1,1);
plot(t,yend); %誤差時(shí)域曲線
title('誤差時(shí)域');
xlabel('時(shí)間');
ylabel('振幅');
subplot(2,1,2);
yp=abs(y1)-abs(fft(x4));%此處要取絕對(duì)值正压，不然會(huì)將相位也算進(jìn)去
plot(f,abs(yp)); %誤差頻譜分布
axis([0 1000 0 10000]);
title('誤差頻譜');
xlabel('F(Hz)');
ylabel('H(jw)');
S=std(yend) %均方值
六欣福、不同濾波器的定義
100Hz帶阻濾波器
function Hd = ww100
Fs = 48000; % Sampling Frequency
Fpass1 = 95; % First Passband Frequency
Fstop1 = 96; % First Stopband Frequency
Fstop2 = 102; % Second Stopband Frequency
Fpass2 = 103; % Second Passband Frequency
Apass1 = 0.1; % First Passband Ripple (dB)
Astop = 70; % Stopband Attenuation (dB)
Apass2 = 0.1; % Second Passband Ripple (dB)
match = 'stopband'; % Band to match exactly
% Construct an FDESIGN object and call its ELLIP method.
h = fdesign.bandstop(Fpass1, Fstop1, Fstop2, Fpass2, Apass1, Astop, ...
Apass2, Fs);
Hd = design(h, 'ellip', 'MatchExactly', match);
500Hz帶阻濾波器
function Hd = ww500
Fs = 48000; % Sampling Frequency
Fpass1 = 497; % First Passband Frequency
Fstop1 = 498; % First Stopband Frequency
Fstop2 = 503; % Second Stopband Frequency
Fpass2 = 504; % Second Passband Frequency
Apass1 = 0.1; % First Passband Ripple (dB)
Astop = 70; % Stopband Attenuation (dB)
Apass2 = 0.1; % Second Passband Ripple (dB)
match = 'stopband'; % Band to match exactly
% Construct an FDESIGN object and call its ELLIP method.
h = fdesign.bandstop(Fpass1, Fstop1, Fstop2, Fpass2, Apass1, Astop, ...
Apass2, Fs);
Hd = design(h, 'ellip', 'MatchExactly', match);
900Hz帶阻濾波器
function Hd = ww900
Fs = 48000; % Sampling Frequency
Fpass1 = 895; % First Passband Frequency
Fstop1 = 897; % First Stopband Frequency
Fstop2 = 903; % Second Stopband Frequency
Fpass2 = 905; % Second Passband Frequency
Apass1 = 0.1; % First Passband Ripple (dB)
Astop = 60; % Stopband Attenuation (dB)
Apass2 = 0.1; % Second Passband Ripple (dB)
match = 'stopband'; % Band to match exactly
% Construct an FDESIGN object and call its ELLIP method.
h = fdesign.bandstop(Fpass1, Fstop1, Fstop2, Fpass2, Apass1, Astop, ...
Apass2, Fs);
Hd = design(h, 'ellip', 'MatchExactly', match);
心得
1、如果題目中給出的是3*pi rad焦履，可以直接使用公式拓劝，也可以換算成Hz.如果換算成HZ,公式如下：
2、其實(shí)實(shí)現(xiàn)帶通和帶阻嘉裤，不過(guò)是將通帶頻率和阻帶頻率換成了坐標(biāo)范圍
思考題
1）本實(shí)驗(yàn)中Butterworth郑临，Chebyshev，Elliptic三種濾波器有什么差別屑宠？對(duì)濾波效果有什么影響厢洞？
相同階數(shù)時(shí)：
巴特沃斯濾波器通帶最平坦，阻帶下降慢。
切比雪夫?yàn)V波器通帶等紋波犀变，阻帶下降較快妹孙。
橢圓濾波器，橢圓濾波器在通帶等紋波（阻帶平坦或等紋波）获枝，阻帶下降最快蠢正。
2）濾波器的阻帶寬度可以設(shè)到多小省店？通帶最小衰減可以設(shè)到多邢浮？阻帶最大衰減可以設(shè)到多大懦傍？
1雹舀、通帶紋波是指在濾波器的頻響中通帶的最大幅值和最小幅值之間的差值，正常的紋波一般小于1db粗俱。不過(guò)也視情況而言说榆，通帶紋波會(huì)導(dǎo)致通帶內(nèi)的幅值大小有變化，一般要求越高寸认，紋波越小越好签财。通帶紋波和濾波器的階數(shù)有關(guān)系，階數(shù)越大紋波越小偏塞。
2唱蒸、 阻帶衰減：在通帶中，有部分信號(hào)通灸叼，部分信號(hào)阻神汹，而阻的部分不能不能全部阻斷，只有部分衰減古今，部分留了下來(lái)屁魏，最小衰減描述了阻礙受阻信號(hào)的能力，衰減越大捉腥，則能力越好蚁堤。
3）在本實(shí)驗(yàn)中，雙線性變換與脈沖不變響應(yīng)對(duì)濾波效果有什么影響但狭？
雙線性變換實(shí)現(xiàn)低通濾波器
%雙線性變換前預(yù)畸
Fs=500;
wp=(100/Fs)*2*pi;
ws=(200/Fs)*2*pi;
Rp=2;
Rs=15;
wp2=2*Fs*tan(wp/2);
ws2=2*Fs*tan(ws/2);
%選擇濾波器的最小階數(shù)
[N,wn]=buttord(wp2,ws2,Rp,Rs,'s');%注意此處輸入的是畸變后的指標(biāo)披诗，輸出N為符合要求的模擬濾波器的最小階數(shù)，wn為3dB帶寬
%創(chuàng)建butterworth模擬濾波器
[Z,P,K]=buttap(N);
%把濾波器零極點(diǎn)模型轉(zhuǎn)化為傳遞函數(shù)模型
[Bap,Aap]=zp2tf(Z,P,K);
%把模擬濾波器原型轉(zhuǎn)換為截止頻率為wn的模擬低通濾波器
[b,a]=lp2lp(Bap,Aap,wn);
%用雙線性法實(shí)現(xiàn)模擬濾波器到數(shù)字濾波器的轉(zhuǎn)換
[bz,az]=bilinear(b,a,Fs);
%繪制頻率響應(yīng)曲線
[H,W]=freqz(bz,az);
subplot(2,1,1);
plot(W/pi,abs(H));
grid;
xlabel('頻率w/pi');
ylabel('幅度絕對(duì)值');
subplot(2,1,2);
plot(W/pi,20*log10((abs(H)+eps)/max(abs(H))));
grid;
xlabel('頻率w/pi');
ylabel('幅度dB');
————————————————
版權(quán)聲明：本文為CSDN博主「Zhang_P_Y」的原創(chuàng)文章立磁，遵循 CC 4.0 BY-SA 版權(quán)協(xié)議呈队，轉(zhuǎn)載請(qǐng)附上原文出處鏈接及本聲明。
原文鏈接：https://blog.csdn.net/lg1259156776/article/details/48603575/
雙線性變換實(shí)現(xiàn)帶通濾波器
https://wenku.baidu.com/view/3f4689221b37f111f18583d049649b6648d70994.html
擴(kuò)
一唱歧、IIR設(shè)計(jì)雙線性變換高通濾波器
https://blog.csdn.net/weixin_30415801/article/details/95458429
f_p = 15000 ; %阻帶上截止頻率
f_st = 30000 ; %通帶下截止頻率
R_p = 0.4 ; %通帶允許的最大衰減
R_st = 30 ; %阻帶允許的最小衰減
f_s = 100000 ; %采樣頻率
T_s = 1 / f_s ; %采樣間隔
%1.將數(shù)字高通濾波器的頻率參數(shù)變換為歸一化的數(shù)字角頻率參數(shù)
omega_p = 2 * pi * f_p / f_s;
omega_st = 2 * pi * f_st / f_s;
%2.預(yù)畸變處理,將歸一化數(shù)字角頻率參數(shù)變換成模擬高通濾波器的角頻率參數(shù)
C = 2 / T_s ;
Omega_p = C * tan( omega_p / 2 );
Omega_st = C * tan( omega_st / 2 );
%3.對(duì)模擬高通濾波器的角頻率參數(shù)做歸一化處理
lamda_p = Omega_p / Omega_p;
lamda_st = Omega_st / Omega_p;
%4.設(shè)計(jì)歸一化模擬濾波器,得到歸一化模擬高通濾波器的轉(zhuǎn)移函數(shù)
[ N , Wn ] = buttord( Omega_p , Omega_st , R_p , R_st , 's' ); %選擇模擬巴特沃斯濾波器的最小階數(shù)
[ z , p , k ] = buttap(N); %創(chuàng)建巴特沃斯模擬低通濾波器
[ Bp , Ap ] = zp2tf(z,p,k); %由零點(diǎn)宪摧、極點(diǎn)粒竖、增益確定傳輸函數(shù)的分子與分母系數(shù)
%5.利用歸一化模擬高通濾波器的轉(zhuǎn)移函數(shù)確定模擬高通濾波器的轉(zhuǎn)移函數(shù)
[ b , a ] = lp2hp(Bp,Ap,Wn);
%6.利用模擬高通濾波器的轉(zhuǎn)移函數(shù)確定IIR數(shù)字濾波器的轉(zhuǎn)移函數(shù)
[ bz , az ] = bilinear(b,a,f_s);
figure(1);
freqz(bz,az);
title('高通濾波器幅度譜和相位譜特性');
設(shè)計(jì)步驟
帶通濾波器
https://wenku.baidu.com/view/3f4689221b37f111f18583d049649b6648d70994.html
https://blog.csdn.net/weixin_30415801/article/details/95458429
f_s1 = 150 ; %阻帶上截止頻率
f_s2 = 600;
f_p1 = 200 ; %通帶下截止頻率
f_p2 = 500;
R_p = 0.5 ; %通帶允許的最大衰減
R_s = 40 ; %阻帶允許的最小衰減
f_s = 2000 ; %采樣頻率
T_s = 1 / f_s ; %采樣間隔
%1.將數(shù)字高通濾波器的頻率參數(shù)變換為歸一化的數(shù)字角頻率參數(shù)
omega_p1 = 2 * pi * f_p1 /f_s;
omega_s1 = 2 * pi * f_s1 /f_s;
omega_p2 = 2 * pi * f_p2 /f_s;
omega_s2 = 2 * pi * f_s2 /f_s;
%2.預(yù)畸變處理,將歸一化數(shù)字角頻率參數(shù)變換成模擬高通濾波器的角頻率參數(shù)
C = 2 / T_s ;
Omega_p1 =C*tan( omega_p1 / 2 );
Omega_s1 = C*tan( omega_s1 / 2 );
Omega_p2 = C*tan( omega_p2 / 2 );
Omega_s2 = C*tan( omega_s2 / 2 );
W0=sqrt(Omega_p1*Omega_p2);
Omega_BW=Omega_p1-Omega_p2;
%3.對(duì)模擬高通濾波器的角頻率參數(shù)做歸一化處理
eta_sl = Omega_sl / Omega_BW;
eta_p1 = Omega_p1 / Omega_BW;
eta_p2 = Omega_p2 / Omega_BW;
eta_s2 = Omega_s2 / Omega_BW;
%4.設(shè)計(jì)歸一化模擬濾波器,得到歸一化模擬高通濾波器的轉(zhuǎn)移函數(shù)
wph=[Omega_p1,Omega_p2 ];
wsh=[Omega_s1,Omega_s2];
[ N , Wn ] = buttord( wph , wsh , R_p , R_st , 's' ); %選擇模擬巴特沃斯濾波器的最小階數(shù)
[ z , p , k ] = buttap(N); %創(chuàng)建巴特沃斯模擬低通濾波器
%[ Bp , Ap ] = zp2tf(z,p,k); %由零點(diǎn)、極點(diǎn)几于、增益確定傳輸函數(shù)的分子與分母系數(shù)
%5.利用歸一化模擬高通濾波器的轉(zhuǎn)移函數(shù)確定模擬高通濾波器的轉(zhuǎn)移函數(shù)
BB=real(poly(z));
AA=real(poly(p));
[ b , a ] = lp2bp(BB,AA,W0,Omega_BW);
%6.利用模擬高通濾波器的轉(zhuǎn)移函數(shù)確定IIR數(shù)字濾波器的轉(zhuǎn)移函數(shù)
[ bz , az ] = bilinear(b,a,f_s);
figure(1);
freqz(bz,az);
title('帶通濾波器幅度譜和相位譜特性');
