實驗一熟悉matlab一、實驗內容練習matlab的使用，熟悉離散卷積運算，產生復指數序列二、實驗目的1、掌握離散卷積計算機實現。2、進一步對離散信號卷積算法的理解。三、原理及算法概要算法：把沖激響應h(n)與輸入序列x(n)分別輸入到程序中，然后調用離散卷積函數y=conv(x.,h)即可得到所要求的結果。原理：離散卷積定義為當序列為有限長時，則 四、程序與運行結果%實現復指數序列程序n=0:30;x=exp(.05+i*pi/6).^n;xr=real(x);xi=imag(x);xm=abs(x);xa=angle(x);figure;subplot(221);stem(n,xr);title('實部');subplot(222);stem(n,xi);title('虛部');subplot(223);stem(n,xm);title('模');subplot(224);stem(n,xa);title('相角');程序運行結果如圖1所示。從圖中可以看出，復指數序列的實部和虛部都是幅度按指數增長的正弦序列。圖1復指數序列波形x1=[1111];nx1=0:3;%實現卷積程序h1=[10.80.640.8^30.8^4];nh1=0:4;

y1=conv(x1,h1);

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

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

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

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

subplot(2,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(2,3,4);stem(x2);title('序列x2');

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

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

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

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

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

圖2卷積程序波形五、結果分析有限長序列的離散卷積計算結果與理論值一致，而存在無限長序列做卷積時，由于在程序處理時是用比較長有限長序列代替的，所以與理論值基本相同。六、實習心得主要進行的是有matlab實現離散卷積的計算，分為三個類型：沖激響應h(n)與輸入序列x(n)都為有限長，一個序列為有限長一個序列為無限長和兩個序列均為無限長。這部分實習內容不難，是原來做過的，主要是為了拾揀起原來所學的知識。第一種類型是最簡單，也是最基本的，直接調用函數y=conv(x.,h)即可。在第二種和第三種類型的計算中遇到了一些困難，在輸入序列時，由于存在無限長的序列，不知道該怎么輸入，查過了相關的題目才弄明白。通過本專題的實習，讓我對數字信號處理中的一些基礎知識有了一些回顧，對原來所學過的知識也熟悉了一些。實驗二FFT一、實驗內容設有離散序列x(n)=0.51sin(2πn/32)+sin(4πn/32)-0.1sin(6πn/32);采集數據長度N=32，分析32點的頻譜，并畫出幅頻特性。二、實驗目的1、了解DFT及FFT的性質和特點2、利用FFT算法計算信號的頻譜。三、關鍵算法算法：讀入離散序列x(n)=0.51sin(2πn/32)+sin(4πn/32)-0.1sin(6πn/32)，采集長度為N=32的數據，調用matlab中的函數fft（x）對其作離散傅里葉變換得到32點的頻譜。原理概要：當抽樣數N=2M時，以下為算法蝶形圖。一般規律如下：當N=2M時，則要進行M次分解，即進行M級蝶形單元的計算2、按自然順序輸入，輸出是碼位倒置。3、每一級包含N/2個基本蝶形運算4、第L級有2L-1個蝶群，蝶群間隔為N/2如果是Matlab實現的話，可用以下兩種方法計算信號頻譜1、調用庫函數為：fft（），直接計算X（k）2、進行矩陣運算四、程序及運行結果n=0:31;x=0.51*sin(2*pi*n/32)+sin(4*pi*n/32)-0.1*sin(6*pi*n/32);y=fft(x);y1=abs(y);y2=angle(y);subplot(2,1,1);stem(y1);subplot(2,1,2);stem(y2);title('Hzhao');typedefstruct{floatreal;floatimag;}COMPLEX;externvoidfft(COMPLEX*x,intm);externvoidifft(COMPLEX*x,intm);voidfft(COMPLEXS*x,intm){staticCOMPLEX*w;staticintmstore=0;staticintn=1;COMPLEXu,temp,tm;COMPLEX*xi,*xip,*xj,*wptr;inti,j,k,l,le,windex;doublearg,w_real,w_imag,wrecur_real,wrecur_imag,wtemp_reaL;if(m!=mstore){if(mstore!=0)free(w);mstore=m;if(m==0)return;n=1<<m;le=n/2;w=(COMPLEX*)colloc(le-1,sizeof(COMPLEX));if(!w){printf("\nUnabletoallocatecomplexwarray\n");exit(1);}arg=4.0*atan(1.0)/le;wrecur_real=w_real=cos(arg);wrecur_imag=w_imag=-sin(arg);xj=w;for(j=1;j<le;j++){xj->real=(float)wrecur_real;xj->imag=(float)wrecur_imag;xj++;wtemp_real=wrecur_real*w_real_wrecur_imag*w_imag;wrecur_imag=wrecur_real*w_imag+wrecur_imag*w_real;wrecur_real=wtemp_real;}}le=n;windex=1;for(l=0;1<m;1++){le=le/2;for(i=0;i<n;i=i+2*le){xi=x+i;xip=xi+le;temp.real=(xi->real+xip->real);temp,imag=(xi->imag+xip->imag);xip->real=(xi->real-xip->real);xip-imag=(xi->imag-xip->imag);xi=temp;}wptr=w+windex-1;for(j=;j<le;j++){u=*wptrfor(i=j;i<n;i=i+2*le){xi=x+i;xip=xi+le;temp.real=(xi->real+xip->real);temp.imag=(xi->imag+xip->imag);tm.real=xi->real-xip->real;tm.imag=xi->imag-xip->imag;xip->real=(tm.real*u.real-tm.imag*u.imag);xip->imag=(tm.real*u,imag+tm.imag*u.real);*xi=temp;}wptr=wptr+windex;}windex=2*windex;}for(i=0;i<n;++i){j=0for(k=0;k<m;++k)j=(j<<1)|(1&(i>>k));if(i<j){xi=x+i;xj=x+j;temp=*xj;*xj=*xi;*xi=temp;}}}voidifft(COMPLEX*x,intm){staticCOMPLEX*w){staticintmstore=0;staticintn=1;COMPLEXu,temp,tm;COMPLEX*xi,*xip,*xj,*wptr;inti,j,k,l,le,windex;floatscale;doublearg,w_real,w_imag,wrecur_real.wrecur_imag,wtemp_real;五、結果分析N點DFT的頻譜分辨率是2π/N。一節指出可以通過補零觀察到更多的頻點，但是這并不意味著補零能夠提高真正的頻譜分辨率。這是因為x[n]實際上是x(t)采樣的主值序列，而將x[n]補零得到的x'[n]周期延拓之后與原來的序列并不相同，也不是x(t)的采樣。因此是不同離散信號的頻譜。對于補零至M點的x'的DFT，只能說它的分辨率2π/M僅具有計算上的意義，并不是真正的、物理意義上的頻譜。頻譜分辨率的提高只能通過提高采樣頻率實現。六、專題實習心得離散傅里葉變換是一種快速算法，由于有限長序列在其頻域也可離散化為有限長序列，因此離散傅里葉變換在數字信號處理中是非常有用的。DFT是重要的變換，在分析有限長序列的有用工具、信號處理的理論上有重要意義、運算方法上起核心作用，譜分析、卷積、相關都可以通DFT在計算機上實現。DFT解決了兩個問題：一是離散與量化，二是快速運算。通過編程實踐體會到了時域、頻域信號的對應關系，也對采樣頻率的含義有了深刻的認識，同時也加深了對采樣信號頻譜周期性的理解。專題三數字濾波器的設計一、實驗內容1、設計一個Butterworth和Chebyshev數字低通濾波器，并學習模擬濾波器轉換數字濾波器的兩種方法，練習課本169頁、181頁、187頁編程題2、分析不同濾波器的特點和結果。3、編程設計實現FIR濾波器。二、實驗目的1.掌握不同IIR濾波器的性質、特點，和用窗函數法設計FIR濾波器的原理和方法。2.通過實驗學習如何設計各種常用的IIR濾波器和FIR濾波器，以便在實際工作中能根據具體情況使用IIR和FIR濾波器。三、原理及算法概要IIR濾波器算法：輸入通帶截止頻率Wp，阻帶截止頻率Ws，通帶衰減Rp，阻帶衰減Rs，通過這些數值調用[NWn]=buttord(Wp,Ws,Rp,Rs)函數計算巴特沃斯數字濾波器的階數N和截止頻率wn，再根據階數N通過函數[b,a]=butter(N,Wn)，即可得到所要的巴特沃斯濾波器。IIR濾波器原理概要：1、濾波器類型1.1Butterworth濾波器Butterworth濾波器的特點是在通帶內的頻率特性是平坦的,并且隨著頻率的增加而衰減。Butterworth濾波器又是最簡單的濾波器。N階低通Butterworth濾波器的幅度平方函數為:1.2ChebyshevⅠ型濾波器ChebyshevⅠ型濾波器的在通帶內的響應是等紋波的，而在阻帶內是單調下降的，或在通帶內是單調下降的，在阻帶內是等紋波的特性。其幅度平方函數為 其中VN是Chebyshev多項式。1.3橢圓濾波器橢圓濾波器在通帶和阻帶內都是等紋波振蕩。橢圓濾波器的特性函數為：其中UN是N階雅可比函數。2、變換方法2.1沖激響應不變法沖激響應不變法的基本原理是從濾波器的沖激響應出發，對模擬濾波器沖激響應h(t)進行取樣，所得到的離散序列h(nT)作為數字濾波器的單位取樣響應。H(z)是由H(s)通過下式的對應關系得到。2.2雙線性變換是在所得到滿足性能指標要求的模擬濾波器的基礎上,通過變換,從而得到相應的數字濾波器。FIR濾波器算法:通過其通帶截止頻率ωp與阻帶截止頻率ωs算出其過渡帶的寬度與濾波器的長度，從而得到理想濾波器的截止頻率，根據所要求的理想濾波器，得到hd（n）。由于其通帶截止頻率處的衰減不大于3分貝與阻帶衰減不小于40分貝，我選擇最接近的漢寧窗，最后調用函數h=hd.*win及freqz(h,1,512)得到實際漢寧窗的響應和實際濾波器的幅度響應。FIR濾波器原理概要：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濾波器，最簡單的窗函數為矩形窗W（n）=1n<(N-1)/20n>(N-1)/2加窗后的頻譜加窗后使實際頻響偏離理想頻響，影響主要有兩個方面：(1)通帶和阻帶之間存在過渡帶，過渡帶寬度取決于窗函數頻響的主瓣寬度。(2)通帶和阻帶區間有紋波，這是由窗函數的旁瓣引起的，旁瓣越多，紋波越多。增加窗函數的寬度N，其主瓣寬度減小，但不改變旁瓣的相對值。為了改善濾波器的性能，要求窗函數的主瓣寬度盡可能窄，以獲得較窄的過渡帶；旁瓣衰減盡可能大，數量盡可能大，從而改善紋波狀況，使實際頻響H（ejω）更好地逼近理想頻響Hd（ejω）2.除了矩形窗外，一般還可以采用以下幾種窗函數Battlet窗：漢寧窗：海明窗布來克曼窗其中：WR是矩形窗的頻譜四、程序及運行結果1.w1=2*400*tan(2*pi*90/(2*400));w2=2*400*tan(2*pi*110/(2*400));wr=2*400*tan(2*pi*120/(2*400));[N,wn]=buttord([w1w2],[100wr],3,10,'s');[B,A]=butter(N,wn,'s');[num,den]=bilinear(B,A,400);[h,w]=freqz(num,den);f=w/pi*200;plot(f,20*log10(abs(h)));axis([40,160,-30,10]);grid;xlabel('頻率/kHz')ylabel('幅度/dB')2.[B,A]=cheby1(3,0.7,2*pi*1000,'s');[num1,den1]=impinvar(B,A,4000);[h1,w]=freqz(num1,den1);[B,A]=cheby1(3,0.7,2/0.00025,'s');[num2,den2]=bilinear(B,A,4000);[h2,w]=freqz(num2,den2);f=w/pi*2000;plot(f,abs(h1),'-',f,abs(h2),'-');grid;xlabel('頻率/Hz')ylabel('幅值/dB')3.w1=2*400*tan(2*pi*90/(2*400));w2=2*400*tan(2*pi*110/(2*400));wr=2*400*tan(2*pi*120/(2*400));[N,wn]=cheb1ord([w1w2],[100wr],3,10,'s');Rp=0.7;[B,A]=cheby1(N,Rp,wn,'s');[num,den]=bilinear(B,A,400);[h,w]=freqz(num,den);f=w/pi*200;plot(f,20*log10(abs(h)));axis([40,160,-30,10]);grid;xlabel('頻率/kHz')ylabel('幅度/dB')4.[B,A]=butter(3,2*pi*1000,'s');%3階濾波器、3分貝截止頻率2*pi*1000[num1,den1]=impinvar(B,A,4000);%沖激響應不變法[h1,w]=freqz(num1,den1);%freqz是繪制一個數字濾波器的頻率響應的函數，返回值h1是通過濾波器后的響應值向量，w是輸入的頻率向量

[B,A]=butter(3,2/0.00025,'s');[num2,den2]=bilinear(B,A,4000);%雙線性變換[h2,w]=freqz(n