專題一離散卷積的計算一、實驗內容設線性時不變（LTI）系統的沖激響應為h(n)，輸入序列為x(n)1、h(n)=(0.8)n，0≤n≤4；x(n)=u(n)-u(n-4)2、h(n)=(0.8)nu(n),x(n)=u(n)-u(n-4)3、h(n)=(0.8)nu(n),x(n)=u(n)求以上三種情況下系統的輸出y(n)。二、實驗目的1、掌握離散卷積計算機實現。2、進一步對離散信號卷積算法的理解。三、原理及算法概要算法：把沖激響應h(n)與輸入序列x(n)分別輸入到程序中，然后調用離散卷積函數y=conv(x.,h)即可得到所要求的結果。原理：離散卷積定義為當序列為有限長時，則 四．理論計算1、h(n)=(0.8)n，0≤n≤4；x(n)=u(n)-u(n-4)(a)當n<0時，y(n)=0(b)當時,(0.8)n(c)當時，(0.8)n(d)當n>7時，y(n)=0理論結果與上圖實驗結果圖中所示吻合。2、h(n)=(0.8)nu(n),x(n)=u(n)-u(n-4)(a)當n<0時，y(n)=0(b)當時,(0.8)n(c)當時，(0.8)n(d)當時，(0.8)n(e)當n>23時，y(n)=0在卷積序列中，因為n趨于無窮，在實驗中無法實現，而0.8的20次冪基本接近于0，所以這里只取到n為20就可以了。由以上可以看出，理論結果與上圖實驗結果圖中所示吻合。3、h(n)=(0.8)nu(n),x(n)=u(n)(a)當n<0時，y(n)=0(b)當時,(0.8)n(c)當時，(0.8)n(d)當n>140時，y(n)=0在卷積序列中，因為n趨于無窮，在實驗中無法實現，而0.8的70次冪的疊加基本接近于最大值，所以這里只取到n為70就可以了。由以上可以看出，理論結果與上圖實驗結果圖中所示吻合。五．程序x1=[1111];nx1=0:3;

h1=[10.80.640.8^30.8^4];nh1=0:4;

y1=conv(x1,h1);

subplot(3,3,1);stem(nx1,x1);title('序列x1');

xlabel('n');ylabel('x1(n)');

subplot(3,3,2);stem(nh1,h1);title('序列h1');

xlabel('n');ylabel('h1(n)');

subplot(3,3,3);stem(y1);title('序列y1');

xlabel('n');ylabel('y1(n)');

x2=[1111];nx2=0:3;

nh2=0:1:20;

h2=(0.8).^nh2;

y2=conv(x2,h2);

subplot(3,3,4);stem(x2);title('序列x2');

xlabel('n');ylabel('x2(n)');

subplot(3,3,5);stem(h2);title('序列h2');

xlabel('n');ylabel('h2(n)');

subplot(3,3,6);stem(y2);title('序列y2');

xlabel('n');ylabel('y2(n)')

nx3=0:1:20;

x3=1.^nx3;

nh3=0:1:20;

h3=(0.8).^nh3;

y3=conv(x3,h3);

subplot(3,3,7);stem(nx3,x3);title('序列x3');

xlabel('n');ylabel('x3(n)');

subplot(3,3,8);stem(nh3,h3);title('序列h3');

xlabel('n');ylabel('h3(n)');

subplot(3,3,9);stem(y3);title('序列y3');

xlabel('n');ylabel('y3(n)'六．程序運行結果七．結果分析有限長序列的離散卷積計算結果與理論值一致，而存在無限長序列做卷積時，由于在程序處理時是用比較長有限長序列代替的，所以與理論值基本相同。八.專題實習心得該專題主要進行的是有matlab實現離散卷積的計算，分為三個類型：沖激響應h(n)與輸入序列x(n)都為有限長，一個序列為有限長一個序列為無限長和兩個序列均為無限長。這部分實習內容不難，是原來做過的，主要是為了拾揀起原來所學的知識。第一種類型是最簡單，也是最基本的，直接調用函數y=conv(x.,h)即可。在第二種和第三種類型的計算中遇到了一些困難，在輸入序列時，由于存在無限長的序列，不知道該怎么輸入，查過了相關的題目才弄明白。通過本專題的實習，讓我對數字信號處理中的一些基礎知識有了一些回顧，對原來所學過的知識也熟悉了一些。專題二離散傅里葉變換及其應用實驗內容設有離散序列x(n)=cos(0.48πn)+cos(0.52πn)分析下列三種情況下的幅頻特性。采集數據長度N=16，分析16點的頻譜，并畫出幅頻特性。采集數據長度N=16，并補零到64點，分析其頻譜，并畫出幅頻特性。采集數據長度N=64，分析46點的頻譜，并畫出幅頻特性。觀察三幅不同的幅頻特性圖，分析和比較它們的特點及形成原因。2、實現序列的內插和抽取所對應的離散傅里葉變換。（選做）求和對應的128點傅里葉變換，比較三個結果所得到的幅頻特性，分析其差別產生的原因。二、實驗目的1、了解DFT及FFT的性質和特點2、利用FFT算法計算信號的頻譜。三、關鍵算法算法1：讀入離散序列x(n)=cos(0.48πn)+cos(0.52πn)，采集長度為N=16的數據，調用matlab中的函數fft（x,16）與fft（x,64）對其作離散傅里葉變換得到16點、64點的頻譜。采集數據長度為N=64，，調用matlab中的函數fft（x,46）對其作離散傅里葉變換得到46點的頻譜。算法2：讀入離散序列信號，每隔N=4采集一個數值，得到新的抽取離散序列，每隔K=0.25采集一個數值得到新的內插離散序列。分別對其做離散傅氏變換，得到對應的頻譜圖。原理概要：當抽樣數N=2M時，以下為算法蝶形圖。 一般規律如下：當N=2M時，則要進行M次分解，即進行M級蝶形單元的計算2、按自然順序輸入，輸出是碼位倒置。3、每一級包含N/2個基本蝶形運算4、第L級有2L-1個蝶群，蝶群間隔為N/2如果是Matlab實現的話，可用以下兩種方法計算信號頻譜1、調用庫函數為：fft（），直接計算X（k）2、進行矩陣運算四．程序程序1:n=0:1:15;

x1=cos(0.48*3.14*n)+cos(0.52*3.14*n);

g1=abs(fft(x1,16));

subplot(3,2,1);stem(x1);title('x1');

subplot(3,2,2);stem(g1);title('g1');

n2=0:1:15;

x2=cos(0.48*3.14*n2)+cos(0.52*3.14*n2);

x2=[x2zeros(1,48)];

g2=abs(fft(x2,64));

subplot(3,2,3);stem(x2);title('x2');

subplot(3,2,4);stem(g2);title('g2');

n3=0:1:64;

x3=cos(0.48*3.14*n3)+cos(0.52*3.14*n3);

g3=abs(fft(x3,46));

subplot(3,2,5);stem(x3);title('x3');

subplot(3,2,6);stem(g3);title('g3');

程序2：x=[];

forn=0:127

x=[x(cos(1/36*pi*n)+cos(1.5/36*pi*n))]

end

X=abs(fft(x));

subplot(3,2,1);stem(x);title('序列x');

xlabel('n');ylabel('x(n)');

subplot(3,2,2);stem(X);title('序列X');

xlabel('\omega');ylabel('X');

x1=[];

forn=0:4:508

x1=[x1(cos(1/36*pi*n)+cos(1.5/36*pi*n))]

end

X1=abs(fft(x1));

subplot(3,2,3);stem(x1);title('序列x1');

xlabel('n');ylabel('x1(n)');

subplot(3,2,4);stem(X1);title('序列X1');

xlabel('\omega');ylabel('X1');

x2=[];

forn=0:0.25:31.75

x2=[x2(cos(1/36*pi*n)+cos(1.5/36*pi*n))]

end

X2=abs(fft(x2));

subplot(3,2,5);stem(x2);title('序列x2');

xlabel('n');ylabel('x2(n)');

subplot(3,2,6);stem(X2);title('序列X2');

xlabel('\omega');ylabel('X2');

五．程序運行結果結果1：結果2：六．結果分析N點DFT的頻譜分辨率是2π/N。一節指出可以通過補零觀察到更多的頻點，但是這并不意味著補零能夠提高真正的頻譜分辨率。這是因為x[n]實際上是x(t)采樣的主值序列，而將x[n]補零得到的x'[n]周期延拓之后與原來的序列并不相同，也不是x(t)的采樣。因此是不同離散信號的頻譜。對于補零至M點的x'的DFT，只能說它的分辨率2π/M僅具有計算上的意義，并不是真正的、物理意義上的頻譜。頻譜分辨率的提高只能通過提高采樣頻率實現。七.專題實習心得離散傅里葉變換是一種快速算法，由于有限長序列在其頻域也可離散化為有限長序列，因此離散傅里葉變換在數字信號處理中是非常有用的。DFT是重要的變換，在分析有限長序列的有用工具、信號處理的理論上有重要意義、運算方法上起核心作用，譜分析、卷積、相關都可以通DFT在計算機上實現。DFT解決了兩個問題：一是離散與量化，二是快速運算。通過編程實踐體會到了時域、頻域信號的對應關系，也對采樣頻率的含義有了深刻的認識，同時也加深了對采樣信號頻譜周期性的理解。八．思考題1、直接計算DFT與用FFT計算，理論上的速度差別應是多少？整個DFT運算總共需要4N2次實數乘法和2N(2N-1)次實數加法。整個FFT運算總共需要(N2/2)+(N/2)=N(N+1)/2≈N2/2次復數乘法和N(N/2-1)+N=N2/2次復數加法。由此可見，理論上用FFT計算運算速度是DFT的二倍。什么是高密度頻譜及高分辨率頻譜？實驗1中至少應采集幾個樣本才能區分兩個頻率分量？DFT的頻譜分辨率一定時，在尾部補零可以得到DFT的高密度頻譜，可以細化當前分辨率下的頻譜，但不能改變DFT的分辨率。高分辨率是指對信號中兩個靠得較近的頻譜分量的識別能力比較高，在采樣頻率不變時，采樣點數N越大，頻譜分辨率越高。在實驗1中至少要采集8個點才能區分兩個頻率分量。專題三IIR濾波器的設計一、實驗內容1、設計一個Butterworth數字低通濾波器，設計指標如下：通帶截止頻率：0.2π，幅度衰減不大于1分貝阻帶截止頻率：0.3π，幅度衰減大于15分貝2、讓不同頻率的正弦波通過濾波器，驗證濾波器性能。3、分析不同濾波器的特點和結果。4、編程設計實現IIR濾波器。二、實驗目的掌握不同IIR濾波器的性質、特點。通過實驗學習如何設計各種常用的IIR濾波器，以便在實際工作中能根據具體情況使用IIR濾波器。三、原理及算法概要算法：輸入通帶截止頻率Wp，阻帶截止頻率Ws，通帶衰減Rp，阻帶衰減Rs，通過這些數值調用[NWn]=buttord(Wp,Ws,Rp,Rs)函數計算巴特沃斯數字濾波器的階數N和截止頻率wn，再根據階數N通過函數[b,a]=butter(N,Wn)，即可得到所要的巴特沃斯濾波器。設計一個正弦波信號，再調用函數A=filter(b,a,I)讓正弦波信號通過濾波器，得到濾波信號。原理概要：1、濾波器類型Butterworth濾波器Butterworth濾波器的特點是在通帶內的頻率特性是平坦的,并且隨著頻率的增加而衰減。Butterworth濾波器又是最簡單的濾波器。N階低通Butterworth濾波器的幅度平方函數為:Chebyshev濾波器Chebyshev濾波器的在通帶內的響應是等紋波的，而在阻帶內是單調下降的，或在通帶內是單調下降的，在阻帶內是等紋波的特性其幅度平方函數為 其中VN是Chebyshev多項式。橢圓濾波器橢圓濾波器在通帶和阻帶內都是等紋波振蕩。橢圓濾波器的特性函數為：其中UN是N階雅可比函數。模擬濾波器設計完成后,就可以進行典型濾波器設計的第二個步驟,通過沖激響應不變法和雙線性變換轉變為相應的數字濾波器。2、變換方法（1）沖激響應不變法沖激響應不變法的基本原理是從濾波器的沖激響應出發，對模擬濾波器沖激響應h(t)進行取樣，所得到的離散序列h(nT)作為數字濾波器的單位取樣響應。H(z)是由H(s)通過下式的對應關系得到。（2）雙線性變換是在所得到滿足性能指標要求的模擬濾波器的基礎上,通過變換,從而得到相應的數字濾波器。四．程序Wp=0.2;

Ws=0.3;

Rp=1;

Rs=15;

[NWn]=buttord(Wp,Ws,Rp,Rs)%用于計算巴特沃斯數字濾波器的階數N和截止頻率wn

[b,a]=butter(N,Wn);%計算N階巴特沃斯數字濾波器系統函數分子、分母多項式的系數向量b、a,設計所需的低通濾波器

[h,omega]=freqz(b,a,512);%返回量h包含了離散系統頻響，調用中若N默認，默認值為512。

plot(omega/pi,20*log10(abs(h)));grid;

xlabel('\omega/\pi');ylabel('Gain,dB');

title('IIRButterworthLowpassFilter');

Wp=0.2;

Ws=0.3;

Rp=1;

Rs=15;

[N1,Wn1]=buttord(Wp,Ws,Rp,Rs);%用于確定階次

[b,a]=butter(N,Wn);%用于直接設計巴特沃茲數字濾波器，即為IIR濾波器

%freqz(b,a);

t=1:300

I=sin(0.1*pi*t)+sin(0.4*pi*t);%設計正弦波

plot(I);

figure;

A=filter(b,a,I);%正弦波通過濾波器

plot(A);

五．程序運行結果N=6Wn=0.2329N1=6Wn1=0.2329六．結果分析Butterworth濾波器在通帶內的頻率特性是平坦的,并且隨著頻率的增加而衰減。正弦信號在經過IIR濾波器濾波后，高頻信號被濾除，低頻信號被保留了下來。七.專題實習心得通過本專題的學習，熟悉和鞏固了Butterworth數字低通濾波器的設計方法和原理，實現濾波器設計的有關經典算法，更重要的是熟練掌握使用MATLAB語言設計各種要求的數字濾波器。這一部分的內容相對來說是有些難度了，做起來花費的精力也多了一些，不過對數字信號處理內容的掌握上又加深了一層。八．思考題1、比較雙線性變換與沖激響應不變法的優缺點。沖激響應不變法的特點：模仿模擬濾波器的單位抽樣相應，時域模仿性好產生頻率響應的混疊失真。保持系統的相位特性換不改變系統的因果性和穩定性;因為存在頻率響應的混疊效應，只適用于限帶的模擬濾波器，高通和帶阻濾波器不宜采用沖激響應不變法，沖激響應不變法是一種時域模仿，缺點是產生頻率響應的混疊失真。沖激響應不變法產生混疊失真的原因是映射是多對一的映射。雙線性變換的特點是一對一的映射，消除了多值變換性，也就消除了頻譜混疊。和頻率變換相比，數字變換有什么優點？數字濾波器的幅頻響應相對于模擬濾波器的幅頻響應有畸變3、實驗中求取相應模擬濾波器原形時，為什么要對模擬頻率Ωc進行歸一化？用歸一化巴特沃思低通濾波器為原型濾波器，一旦歸一化低通濾波器的系統函數確定后，其它巴特沃思低通濾波、高通、帶通、帶阻濾波器的傳遞函數都可以通過變換法從歸一化低通原型的傳遞函數Ha(S)得到。歸一化原型濾波器是指截止頻率Ωc已經歸一化成Ωc’=1的低通濾波器。對于截止頻率為某個Ωc的低通濾波器，則令S/Ωc代替歸一化原型濾波器系統函數中的S∧，即S∧S/Ωc對于其他高通、帶通、帶阻濾波器，可應用后面討論到的頻帶變換法，由其變換得出。專題四用窗函數設計FIR濾波器一、實驗內容選取合適窗函數設計一個線性相位FIR低通濾波器，使它滿足如下性能指標：通帶截止頻率：ωp=0.5π,通帶截止頻率處的衰減不大于3分貝；阻帶截止頻率：ωs=0.66π，阻帶衰減不小于40分貝。二、實驗目的1、掌握用窗函數法設計FIR濾波器的原理和方法。2、熟悉線性相位濾波器特性。3、了解各種窗函數對濾波器特性的影響。三、原理及算法概要算法:通過其通帶截止頻率ωp與阻帶截止頻率ωs算出其過渡帶的寬度與濾波器的長度，從而得到理想濾波器的截止頻率，根據所要求的理想濾波器，得到hd（n）。由于其通帶截止頻率處的衰減不大于3分貝與阻帶衰減不小于40分貝，我選擇最接近的漢寧窗，最后調用函數h=hd.*win及freqz(h,1,512)得到實際漢寧窗的響應和實際濾波器的幅度響應。原理概要：1、利用窗函數法設計FIR濾波器FIR濾波器的最大特點是其相位特性可以設計為嚴格的線性，而其幅值可以任意設置，這樣輸出波形就不會相位失真。理想低通濾波器的單位取樣響應hd（n）是無限長的，所以要用一個有限長的因果序列h（n）進行逼近，最有效的方法是截斷hd（n），即用有限長的窗函數w（n）來截取hd（n），表示為h（n）=hd（n）w（n）為獲得線性相位的FIR濾波器，h（n）必須滿足中心對稱條件，序列h（n）應有一定的延遲α，且α=（N-1）/2頻率響應逼近hd（ejw）的FIR濾波器，最簡單的窗函數為矩形窗1n<(N-1)/2W（n）=0n>(N-1)/2加窗后的頻譜加窗后使實際頻響偏離理想頻響，影響主要有兩個方面：通帶和阻帶之間存在過渡帶，過渡帶寬度取決于窗函數頻響的主瓣寬度。通帶和阻帶區間有紋波，這是由窗函數的旁瓣引起的，旁瓣越多，紋波越多。增加窗函數的寬度N，其主瓣寬度減小，但不改變旁瓣的相對值。為了改善濾波器的性能，要求窗函數的主瓣寬度盡可能窄，以獲得較窄的過渡帶；旁瓣衰減盡可能大，數量盡可能大，從而改善紋波狀況，使實際頻響H（ejω）更好地逼近理想頻響Hd（ejω）除了矩形窗外，一般還可以采用以下幾種窗函數Battlet窗：漢寧窗：海明窗布來克曼窗其中：WR是矩形窗的頻譜四．程序wp=0.5*pi;

ws=0.66*pi;

wdelta=ws-wp;%過渡帶寬度

N=ceil(8*pi/wdelta)%濾波器長度

ifrem(N,2)==0

N=N+1;

end

Nw=N;

wc=(wp+ws)/2;%理想低通濾波器的截止頻率

n=0:N-1;

alpha=(N-1)/2;

m=n-alpha+0.00001;

hd=sin(wc*m)./(pi*m);%一個響應

win=(hanning(Nw))';%漢寧窗

h=hd.*win;%實際漢寧窗的響應

freqz(h,1,512);%實際濾波器的幅度響應

五．程序運行結果N=50六．結果分析上圖中第一個圖為FIR濾波器的相頻特性，第二個圖為幅頻特性；根據濾波器特性：通帶截止頻率處的衰減不大于3分貝；阻帶衰減不小于40分貝。選用漢寧窗。FIR濾波器的最大特點是其相位特性可以設計為嚴格的線性，而其幅值可以任意設置，這樣輸出波形就不會相位失真。理想低通濾波器的單位取樣響應hd（n）是無限長的，所以要用一個有限長的因果序列h（n）進行逼近，最有效的方法是截斷hd（n），即用有限長的窗函數w（n）來截取hd（n），表示為h（n）=hd（n）w（n）為獲得線性相位的FIR濾波器，h（n）必須滿足中心對稱條件，序列h（n）應有一定的延遲α，且α=（N-1）/2頻率響應逼近hd（ejw）的FIR濾波器，最簡單的窗函數為矩形窗1n<(N-1)/2W（n）=0n>(N-1)/2加窗后的頻譜加窗后使實際頻響偏離理想頻響，影響主要有兩個方面：通帶和阻帶之間存在過渡帶，過渡帶寬度取決于窗函數頻響的主瓣寬度。通帶和阻帶區間有紋波，這是由窗函數的旁瓣引起的，旁瓣越多，紋波越多。增加窗函數的寬度N，其主瓣寬度減小，但不改變旁瓣的相對值。為了改善濾波器的性能，要求窗函數的主瓣寬度盡可能窄，以獲得較窄的過渡帶；旁瓣衰減盡可能大，數量盡可能大，從而改善紋波狀況，使實際頻響H（ejω）更好地逼近理想頻響Hd（ejω）七.專題實習心得FIR有限沖擊響應數字濾波器由于設計靈活，在濾波效果和響應時間之間能較好的折衷，因此在數字信號處理領域應用廣泛。用窗函數法設計FIR數字濾波器在截斷點依舊取模擬樣本沖激響應的全值更利于理解和仿真。但是，窗函數法設計FIR數字濾波器的頻率特性與模擬樣本頻率特性不完全一樣，且具有線性相移的特點。窗函數法設計FIR數字濾波器是傅立葉變換的典型運用。通過學習有益于加深對數字濾波器的理解，并提高使用傅立葉變換分析解決實際問題的能力。八：思考題在這個實驗中，你選擇哪一種窗函數？試說明原因。選用漢寧窗。根據所給定濾波器特性：通帶截止頻率處的衰減不大于3分貝；阻帶衰減不小于40分貝。參照書本上的各濾波器的特性。2、用窗函數法設計FIR濾波器時，濾波器的過渡帶寬度和阻帶衰減和哪些因素有關？窗函數的主瓣寬度決定濾波器過渡帶寬度，旁瓣幅值決定阻帶衰減。選擇窗函數應使其頻譜滿足：主瓣寬度盡可能的窄，以使過渡帶盡可能陡；最大旁瓣相對于主瓣盡可能小，使能量盡可能集中在主瓣內，這樣可以使肩峰和波動減小。綜合一、實驗內容錄制一段自己的語音信號，時間為10s左右，并對錄制的信號進行采樣；畫出采樣后語音信號的時域波形和頻譜圖；給定濾波器的性能指標，采用窗函數法或雙線性變換設計濾波器，并畫出濾波器的頻率響應；然后用自己設計的濾波器對采集的語音信號進行濾波，畫出濾波后信號的時域波形和頻譜，并對濾波前后的信號進行對比，分析信號的變化；回放語音信號；最后，用MATLAB設計一信號處理系統界面。二．算法調用函數functionpushbutton1_Callback(hObject,eventdata,handles)實現一個信號處理系統界面。選擇左鍵時，用雙線性變換法設計濾波器來對信號進行處理，選擇右鍵時，用窗函數法設計濾波器來對信號進行處理。讀取語音信號，對語音信號進行f=8000的頻率進行采樣，調用函數y1=fft(x1,2048)對所采集的點做2048點FFT變換。先設計butterworth模擬濾波器，再用雙線性變換法實現模擬濾波器到數字濾波器的轉換。最后調用函數f1=filter(bz,az,x2)對加了噪聲的語音信號進行濾波，得到濾波后的頻譜圖。三．程序functionpushbutton1_Callback(hObject,eventdata,handles)

%hObjecthandletopushbutton1(seeGCBO)

%eventdatareserved-tobedefinedinafutureversionofMATLAB

%handlesstructurewithhandlesanduserdata(seeGUIDATA)

fs=8000;%語音信號采樣頻率為8000

x1=wavread('F:\實習\信息工程實習\第五題\語音.wav');

t=(0:length(x1)-1)/8000;

y1=fft(x1,2048);%對信號做2048點FFT變換

f=fs*(0:1023)/2048;

figure(1)

subplot(2,2,1)

plot(t,x1)%做原始信號的時域波形

gridon;axistight;

title('原始語音信號');

xlabel('time(s)');

ylabel('幅度');

subplot(2,2,2)

plot(f,abs(y1(1:1024)))%做原始信號的FFT頻譜

gridon;axistight;

title('原始語音信號的FFT頻譜')

xlabel('Hz');

ylabel('幅度');

%雙線性變換法設計的巴特沃斯濾波器

A1=0.05;A2=0.10;

d=[A1*cos(2*pi*3800*t)+A2*sin(2*pi*3600*t)]';

x2=x1+d;

wp=0.8*pi;

ws=0.85*pi;

Rp=1;

Rs=15;

Fs=8000;

Ts=1/Fs;

wp1=2/Ts*tan(wp/2);%將模擬指標轉換為數字指標

ws1=2/Ts*tan(ws/2);

[N,Wn]=buttord(wp1,ws1,Rp,Rs,'s');%選擇濾波器最小階數

[Z,P,K]=buttap(N);%創建butterworth模擬濾波器

[Bap,Aap]=zp2tf(Z,P,K);

[b,a]=lp2lp(Bap,Aap,Wn);

[bz,az]=bilinear(b,a,Fs);%用雙線性法實現模擬到數字的轉換

[H,W]=freqz(bz,az);%繪制頻率響應曲線

subplot(2,2,3)

plot(W*Fs/(2*pi),abs(H))

gridon;axistight;

xlabel('頻率（Hz）')

ylabel('頻率響應')

title('Butterworth')

f1=filter(bz,az,x2);

figure(2)

subplot(2,2,1)

plot(t,x2)%畫出濾波前的時域圖

gridon;axistight;

title('濾波前的時域波形');

subplot(2,2,2)

plot(t,f1);%畫出濾波后的時域圖

gridon;axistight;

title('濾波后的時域波形');

y3=fft(f1,2048);

y2=fft(x2,2048);

subplot(2,2,3);

plot(f,abs(y2(1:1024)));%畫出濾波前的頻譜圖

gridon;axistight;

title('濾波前的頻譜')

xlabel('Hz');

ylabel('幅度');

subplot(2,2,4)

plot(f,abs(y3(1:1024)));%畫出濾波后的頻譜圖

gridon;axistight;

title('濾波后的頻譜')

xlabel('Hz');

ylabel('幅度');

functionpushbutton2_Callback(hObject,eventdata,handles)

%hObjecthandletopushbutton2(seeGCBO)

%eventdatareserved-tobedefinedinafutureversionofMATLAB

%handlesstructurewithhandlesanduserdata(seeGUIDATA)

fs=8000;%語音信號采樣頻率為8000

x1=wavread('F:\實習\信息工程實習\第五題\語音.wav');

t=(0:length(x1)-1)/8000;

y1=fft(x1,2048);%對信號做2048點FFT變換

f=fs*(0:1023)/2048;

figure(1)

subplot(2,1,1)

plot(t,x1)%做原始信號的時域波形

gridon;axistight;

title('原始語音信號');

xlabel('time(s)');

ylabel('幅度');

subplot(2,1,2)

plot(f,abs(y1(1:1024)))%做原始信號的FFT頻譜

gridon;axistight;

title('原始語音信號的FFT頻譜')

xlabel('Hz');

ylabel('幅度');

%窗函數設計濾波器

t=(0:length(x1)-1)/8000;

f=fs*(0:2047)/4096;

A1=0.05;A2=0.10;

d=[A1*cos(2*pi*3600*t)+A2*sin(2*pi*3800*t)]';

x2=x1+d;

wp=0.8*pi;

ws=0.85*pi;

wdelta=ws-wp;

N=ceil(6.6*pi/wdelta);%取整

wn=(0.8+0.85)*pi/2;

[bz,az]=fir1(N,wn/pi,hamming(N+1));%選擇窗函數，并歸一化截止頻率

figure(2)

freqz(bz,