實驗五(信號抽樣與恢復(fù))和實驗六(FFT算法的應(yīng)用)

實驗五信號抽樣與恢復(fù)一、實驗?zāi)康膶W(xué)會用MATLAB實現(xiàn)連續(xù)信號的采樣和重建二、實驗原理1．抽樣定理若是帶限信號，帶寬為,經(jīng)采樣后的頻譜就是將的頻譜在頻率軸上以采樣頻率為間隔進(jìn)行周期延拓。因此，當(dāng)時，不會發(fā)生頻率混疊；而當(dāng)<時將發(fā)生頻率混疊。2．信號重建經(jīng)采樣后得到信號經(jīng)理想低通則可得到重建信號，即：=*其中：==所以：=*=*=上式表明，連續(xù)信號可以展開成抽樣函數(shù)的無窮級數(shù)。利用MATLAB中的來表示，有，所以可以得到在MATLAB中信號由重建的表達(dá)式如下：=我們選取信號=作為被采樣信號，當(dāng)采樣頻率=2時，稱為臨界采樣。我們?nèi)±硐氲屯ǖ慕刂诡l率=。下面程序?qū)崿F(xiàn)對信號=的采樣及由該采樣信號恢復(fù)重建：例5-1Sa(t)的臨界采樣及信號重構(gòu)；wm=1;%信號帶寬wc=wm;%濾波器截止頻率Ts=pi/wm;%采樣間隔ws=2*pi/Ts;%采樣角頻率n=-100:100;%時域采樣電數(shù)nTs=n*Ts%時域采樣點f=sinc(nTs/pi);Dt=0.005;t=-15:Dt:15;fa=f*Ts*wc/pi*sinc((wc/pi)*(ones(length(nTs),1)*t-nTs'*ones(1,length(t))));%信號重構(gòu)t1=-15:0.5:15;f1=sinc(t1/pi);subplot(211);stem(t1,f1);xlabel('kTs');ylabel('f(kTs)');title('sa(t)=sinc(t/pi)的臨界采樣信號');subplot(212);plot(t,fa)xlabel('t');ylabel('fa(t)');title('由sa(t)=sinc(t/pi)的臨界采樣信號重構(gòu)sa(t)');grid;例5-2Sa(t)的過采樣及信號重構(gòu)和絕對誤差分析程序和例4-1類似，將采樣間隔改成Ts=0.7*pi/wm,濾波器截止頻率該成wc=1.1*wm，添加一個誤差函數(shù)wm=1;wc=1.1*wm;Ts=0.7*pi/wm;ws=2*pi/Ts;n=-100:100;nTs=n*Tsf=sinc(nTs/pi);Dt=0.005;t=-15:Dt:15;fa=f*Ts*wc/pi*sinc((wc/pi)*(ones(length(nTs),1)*t-nTs'*ones(1,length(t))));error=abs(fa-sinc(t/pi));%重構(gòu)信號與原信號誤差t1=-15:0.5:15;f1=sinc(t1/pi);subplot(311);stem(t1,f1);xlabel('kTs');ylabel('f(kTs)');title('sa(t)=sinc(t/pi)的采樣信號');subplot(312);plot(t,fa)xlabel('t');ylabel('fa(t)');title('由sa(t)=sinc(t/pi)的過采樣信號重構(gòu)sa(t)');grid;subplot(313);plot(t,error);xlabel('t');ylabel('error(t)');title('過采樣信號與原信號的誤差error(t)');例5-3Sa(t)的欠采樣及信號重構(gòu)和絕對誤差分析程序和例4-2類似，將采樣間隔改成Ts=1.5*pi/wm,濾波器截止頻率該成wc=wm=1三、上機(jī)實驗內(nèi)容1．驗證實驗原理中所述的相關(guān)程序；2．設(shè)f(t)=0.5*(1+cost)*(u(t+pi)-u(t-pi))，由于不是嚴(yán)格的頻帶有限信號，但其頻譜大部分集中在[0，2]之間，帶寬wm可根據(jù)一定的精度要求做一些近似。試根據(jù)以下兩種情況用MATLAB實現(xiàn)由f(t)的抽樣信號fs(t)重建f(t)并求兩者誤差，分析兩種情況下的結(jié)果。(1)wm=2,wc=1.2wm,Ts=1;(2)wm=2,wc=2,Ts=2.5實驗六FFT算法的應(yīng)用一、實驗?zāi)康模杭由顚﹄x散信號的DFT的理解及其FFT算法的運(yùn)用。二、實驗原理和例子：N點序列的DFT和IDFT變換定義式如下： ，利用旋轉(zhuǎn)因子具有周期性，可以得到快速算法（FFT）。在MATLAB中，可以用函數(shù)X=fft（x，N）和x=ifft（X，N）計算N點序列的DFT正、反變換。例1對連續(xù)的單一頻率周期信號按采樣頻率采樣，截取長度N分別選N=20和N=16，觀察其DFT結(jié)果的幅度譜。解此時離散序列，即k=8。用MATLAB計算并作圖，函數(shù)fft用于計算離散傅里葉變換DFT，程序如下：k=8;n1=[0:1:19];xa1=sin(2*pi*n1/k);subplot(2,2,1)plot(n1,xa1)xlabel('t/T');ylabel('x(n)');xk1=fft(xa1);xk1=abs(xk1);subplot(2,2,2)stem(n1,xk1)xlabel('k');ylabel('X(k)');n2=[0:1:15];xa2=sin(2*pi*n2/k);subplot(2,2,3)plot(n2,xa2)xlabel('t/T');ylabel('x(n)');xk2=fft(xa2);xk2=abs(xk2);subplot(2,2,4)stem(n2,xk2)xlabel('k');ylabel('X(k)');

計算結(jié)果示于圖2.1，(a)和(b)分別是N=20時的截取信號和DFT結(jié)果，由于截取了兩個半周期，頻譜出現(xiàn)泄漏；(c)和(d)分別是N=16時的截取信號和DFT結(jié)果，由于截取了兩個整周期，得到單一譜線的頻譜。上述頻譜的誤差主要是由于時域中對信號的非整周期截斷產(chǎn)生的頻譜泄漏。實驗內(nèi)容：2N點實數(shù)序列N=64。用一個64點的復(fù)數(shù)FFT程序，一次算出，并繪出。實驗要求：利用MATLAB編程完成計算，繪出相應(yīng)圖形。并與理論計算相比較，說明實驗結(jié)果的原因。用以下代碼實現(xiàn)可得圖6-1所示的DFT圖。>>N=64;>>n=0:2*N-1;>>x=cos(2*pi*7*n/N)+1/2*cos(2*pi*19*n/N);>>X=fft(x,128);>>k=n;>>stem(k,abs(X))>>grid>>xlabel('k');ylabel('|X[k]|');圖6-1理論分析如下：由歐拉公式得：對，其2N點的DFT變換為：當(dāng)時，=0當(dāng)時，即由此可得當(dāng)k=14,38,90,114時有值其余為0），與圖6-1有相同的結(jié)論。例2考慮取時，求的DFT：X（k）；將（1）中的以補(bǔ)零方式使加長到，求X（k）；（3）取，求X（k）要求畫出和X（k），并比較（1）~（3）的結(jié)果。解：首先定義DFT和IDFTfunction[XK]=dft(xn,N)n=[0:1:N-1];k=[0:1:N-1];WN=exp(-j*2*pi/N);nk=n'*k;WNnk=WN.^nk;XK=xn*WNnk;function[xn]=idft(Xk,N)n=[0:1:N-1];k=[0;1:N-1];WN=exp9(-j*pi/N);nk=n'*k;WNnk=WN.^(-nk);xn=(Xk*WNnk)/N;（1）x（n）的10點DFTn=[0:1:99];x=cos(0.48*pi*n)+cos(0.52*pi*n);n1=[0:1:9];y1=x(1:1:10);subplot(211);stem(n1,y1);title('signx(n),n[0,9]');axis([0,10,-2.5,2.5]);text(10.2,-2.5,'n')Y1=dft(y1,10);magY1=abs(Y1(1:1:6));k1=0:1:5;w0=pi*2/10*k1;subplot(212);stem(w0/pi,magY1);title('在頻域dtft分析');xlabel('單位pi')由于采樣的頻率太小，在上圖中無法確定x（n）由兩個頻率成分組成的。以補(bǔ)零方式將x（n）加長到100個樣本n2=[0:1:99];y2=[x(1:1:10)zeros(1,90)];subplot(211);stem(n2,y2);title('信號x(n),0<=n<=9+90zeros');Y2=dft(y2,100);magY2=abs(Y2(1:1:51));k2=0:1:50;w2=2*pi/100*k2;subplot(212);plot(w2/pi,magY2);title('DTFT幅度');xlabel('頻率(單位:pi)')實驗內(nèi)容：用fft函數(shù)代替例2中的dft函數(shù)后，看看運(yùn)行結(jié)果如何？例3:用FFT分析信號頻率成分一被噪聲污染的信號,很難看出它所包含的頻率分量,如一個由50Hz和150Hz正弦信號構(gòu)成的信號,受到均值為零、均方差為0.5的高斯隨機(jī)信號的于擾,數(shù)據(jù)采樣率fs=500Hz.通過FFT來分析其信號頻率成分,用matlab實現(xiàn)如下:fs=500;%采樣頻率fs=500Hz.t=0:1/fs:1;%采樣周期為1/fs.f=sin(2*pi*50*t)+sin(2*pi*150*t);%產(chǎn)生信號f(t)subplot(3,1,1);plot(t,f);title('原始信號');y=f+0.5*randn(1,length(t));%加噪subplot(3,1,2);plot(t,y);title('受噪聲污染的信號');N=256;Y=fft(y,N);%對加噪信號進(jìn)行FFTk=0:N-1;f=fs*k/N;subplot(3,1,3);plot(f,abs(Y));title('FFT(幅度譜)');(由頻譜圖可見,在50Hz和150Hz各出現(xiàn)很長的譜線,表明含噪信號y中含有這二個頻率的信號.在350Hz和450Hz處也出現(xiàn)很長的譜線,這并不是說y中也含350Hz和450Hz的信號,這是由于采樣信號的頻譜是以采樣頻率fs為間隔周期出現(xiàn)而造成的)注意:當(dāng)采樣頻率fs>2fm=2*150=300Hz時,滿足奈奎斯特抽樣定理條件,不會產(chǎn)生頻譜混迭現(xiàn)象.當(dāng)fs<300Hz時則會產(chǎn)生頻譜混迭現(xiàn)象.s>2不發(fā)生混疊s<2發(fā)生混疊s=2臨界例4傅里葉變換的頻移特性

若F(f(t)]=F(ω),則例:設(shè)f(t)=sin(400t),ω0=200.fs=1000;%采樣頻率fs=1000Hz.t=0:1/fs:1;y1=sin(400*pi*t);y2=sin(400*pi*t).*exp(j*200*pi*t);N=512;Y1=fft(y1,N);Y2=fft(y2,N);subplot(3,1,1);plot(t,y1);k=0:N-1;f=fs*k/N;subplot(3,1,2);plot(f,abs(Y1));subplot(3,1,3