《信號(hào)與系統(tǒng)實(shí)訓(xùn)指導(dǎo)》課件第4章

第4章信號(hào)和系統(tǒng)的頻域分析4.1相關(guān)MATLAB函數(shù)

4.2周期信號(hào)的傅里葉級(jí)數(shù)

4.3連續(xù)時(shí)間信號(hào)的傅里葉變換

4.4離散時(shí)間信號(hào)的傅里葉變換

4.5系統(tǒng)的頻域分析實(shí)例習(xí)題本章介紹信號(hào)的傅里葉變換和系統(tǒng)的頻域分析以及如何用MATLAB進(jìn)行信號(hào)和系統(tǒng)的頻域分析。4.1節(jié)簡(jiǎn)要介紹頻域分析用到的MATLAB函數(shù)。4.2節(jié)介紹用MATLAB求周期信號(hào)的傅里葉級(jí)數(shù)。4.3節(jié)介紹用MATLAB求連續(xù)時(shí)間信號(hào)的傅里葉變換。4.4節(jié)介紹用MATLAB求離散時(shí)間信號(hào)的傅里葉變換。4.5節(jié)為系統(tǒng)的頻域分析實(shí)例。

信號(hào)與系統(tǒng)的頻域分析涉及到的MATLAB函數(shù)主要有：

fft

ifft

freqs

freqz

下面簡(jiǎn)要介紹這些函數(shù)。4.1相關(guān)MATLAB函數(shù)

1.fft

函數(shù)fft用來計(jì)算離散傅里葉變換,這個(gè)函數(shù)能夠進(jìn)行序列的快速傅里葉變換。如果把序列x(n)的快速傅里葉變換的結(jié)果存在向量X中,可以表示為

X＝fft(x)為了提高計(jì)算效率,向量x的長(zhǎng)度取2的冪次,如64、128、512、1024、2048等。為了得到適當(dāng)?shù)拈L(zhǎng)度,可以將向量x補(bǔ)零。MATLAB可以用下面的命令自動(dòng)完成這個(gè)工作,其中N要定義為2的冪次：

X＝fft(x,N)

向量x的長(zhǎng)度越長(zhǎng),快速傅里葉變換的坐標(biāo)網(wǎng)格就越好。

2.ifft

ifft函數(shù)用于計(jì)算傅里葉反變換：

x＝ifft(X)

3.freqs

函數(shù)freqs可以直接計(jì)算系統(tǒng)的頻率響應(yīng)。要求系統(tǒng)的頻率響應(yīng),可以將系統(tǒng)的頻率響應(yīng)函數(shù)H(jω)的分子和分母的系數(shù)保存到向量num和den中?？梢远x一個(gè)向量w包含需要計(jì)算頻率響應(yīng)的頻率范圍,例如w＝a∶b∶c,其中a是最低頻率,c是最高頻率,b是頻率增量。如果沒有輸出參數(shù),直接調(diào)用

freqs(num,den,w)則MATLAB會(huì)在當(dāng)前繪圖窗口自動(dòng)畫出幅頻特性和相頻特性的曲線。要是不定義w,命令

freqs(num,den)

自動(dòng)選擇畫出H(jω)的頻率范圍。如果命令是

H＝freqs(num,den,w)

向量H為返回值,包含的是w中每個(gè)頻率點(diǎn)對(duì)應(yīng)的頻率響應(yīng)值,它是一個(gè)復(fù)向量,如果需要得到系統(tǒng)的幅頻特性和相頻特性,還需要用abs和angle函數(shù)分別求得。

4.freqz

函數(shù)freqz可以用來計(jì)算離散系統(tǒng)的頻率響應(yīng),即求得離散系統(tǒng)的沖激響應(yīng)的離散時(shí)間傅里葉變換(DTFT)。將系統(tǒng)的頻率響應(yīng)函數(shù)H(ejω)的分子和分母的系數(shù)保存到向量num和den中,命令

[Homega］＝freqz(num,den,N)；計(jì)算出在0和π之間的N個(gè)等分頻率(即ωk＝(π/N)k,0≤k≤N－1)上的頻率響應(yīng)H(ejω)。這里系數(shù)向量num和den的格式與函數(shù)filter的要求相同。輸出向量H保存的是H(ejωk),向量omega保存的是頻率ωk。如果需要得到系統(tǒng)的幅頻特性和相頻特性,可以用abs(H)和angle(H)分別求得。命令

[Homega］＝freqz(num,den,N,′whole′)；能計(jì)算出在0和2π之間的N個(gè)等分頻率(即ωk＝(2π/N)k,0≤k≤N－1)上的頻率響應(yīng)H(ejω)。如果需要自己定義所需的頻率范圍,可以用命令

H＝freqz(num,den,omega)；

例如omega＝－pi∶2*pi/100∶pi。

如果連續(xù)信號(hào)x(t)是周期信號(hào),那么對(duì)于一切t,存在某個(gè)正值T,有

(4.1)

滿足上式的最小正值T稱為信號(hào)的基本周期,而ω0＝2π/T稱為基本頻率。4.2周期信號(hào)的傅里葉級(jí)數(shù)周期信號(hào)可以表示為三角函數(shù)或復(fù)指數(shù)函數(shù)的線性組合,稱為周期信號(hào)x(t)的傅里葉級(jí)數(shù)。周期信號(hào)x(t)表示為復(fù)指數(shù)信號(hào)和的形式為

(4.2)

式中,c0是實(shí)數(shù),ck(k≠0)一般是復(fù)數(shù)。復(fù)指數(shù)項(xiàng)系數(shù)ck可以由下式計(jì)算：

(4.3)

c0項(xiàng)是常數(shù)項(xiàng),即x(t)的直流分量,可以把式(4.3)簡(jiǎn)化為下式求得：

(4.4)

周期信號(hào)x(t)表示為三角函數(shù)和的形式為

(4.5)式中,c0是信號(hào)的直流分量, 是信號(hào)的諧波分量,|ck|和φk分別為ck的模和輻角。

從式(4.2)至式(4.5)可以看出,周期信號(hào)的傅里葉級(jí)數(shù)是把信號(hào)表示為不同頻率的正弦信號(hào)或復(fù)指數(shù)信號(hào)的疊加,這些分量的振幅和相位都是頻率的函數(shù),隨著頻率的變化而變化。以ω作為變量,可以畫出振幅|ck|或2|ck|隨ω變化的曲線,以及相位隨ω變化的曲線。由于x(t)的頻率分量是離散的,在每個(gè)頻率點(diǎn)上可以畫一條對(duì)應(yīng)ωk的垂線,線的長(zhǎng)度就代表振幅或相位的大小,這種圖稱為線譜圖。振幅隨ω變化的線譜圖稱為振幅頻譜；相位隨ω變化的線譜圖稱為相位頻譜。通過周期信號(hào)的譜圖,可以從頻率分量的角度來分析信號(hào)。

例4-1繪制周期鋸齒信號(hào)的頻譜圖。

設(shè)有一個(gè)周期鋸齒信號(hào)x(t),信號(hào)的基本周期為T＝2(ω0＝2π/T＝π),信號(hào)的波形如圖4-1所示。畫出信號(hào)的振幅頻譜和相位頻譜。

圖4-1周期鋸齒信號(hào)解:根據(jù)式(4.4)得直流分量c0為

根據(jù)式(4.3)得ck為

所以信號(hào)的振幅頻譜為

ck的輻角φk可以表示為

所以信號(hào)的相位頻譜為

根據(jù)以上分析結(jié)果,畫出信號(hào)振幅頻譜的MATLAB程序如下:

%振幅頻譜

N＝input(′級(jí)數(shù)諧波項(xiàng)數(shù)N的值:′)；

c0＝0；

w0＝pi；

stem(0,c0)

holdon

fork＝1∶N,

ck＝1/k/pi；

c_k＝ck；

stem([－k*w0］,[c_k］)

stem([k*w0］,[ck］)

end

ylabel(′|c_k|′)

xlabel(′Frequency(rad/sec)′)根據(jù)定義周期信號(hào)的級(jí)數(shù)項(xiàng)數(shù)是無窮多項(xiàng),但是MATLAB程序只能表示出有限多項(xiàng),所以在程序中加入了設(shè)置項(xiàng)數(shù)的input命令:N＝input(′級(jí)數(shù)諧波項(xiàng)數(shù)N的值:′),運(yùn)行程序后,會(huì)在命令窗口顯示提示語句:

級(jí)數(shù)諧波項(xiàng)數(shù)N的值:

在提示語句后面輸入一個(gè)正整數(shù),回車即可運(yùn)行程序。如輸入項(xiàng)數(shù)為“6”,則結(jié)果如圖4-2所示。

圖4-2振幅頻譜畫出信號(hào)相位頻譜的MATLAB程序如下:

%相位頻譜

N＝input(′級(jí)數(shù)諧波項(xiàng)數(shù)N的值:′)；

w0＝pi；

stem(0,0)

holdon

fork＝1∶N,

theta＝－pi/2；

theta＝theta*180/pi；

theta_a＝－theta；

stem([－k*w0］,[theta_a］)

stem([k*w0］,[theta］)

end

xlabel(′Frequency(rad/sec)′)

ylabel(′＼phi_k(degrees)′)

該程序中也加入了設(shè)置項(xiàng)數(shù)的命令,輸入項(xiàng)數(shù)為“6”,運(yùn)行結(jié)果如圖4-3所示。

圖4-3相位頻譜請(qǐng)讀者自己畫出周期函數(shù)三角級(jí)數(shù)的頻譜。

例4-2用級(jí)數(shù)表示周期鋸齒信號(hào)并觀察吉布斯現(xiàn)象。

用級(jí)數(shù)表示例4-1中給出的周期鋸齒信號(hào),畫出N＝5,15,45時(shí)的波形,并觀察吉布斯現(xiàn)象。

解:用傅里葉級(jí)數(shù)表示周期信號(hào)時(shí),需要無限多項(xiàng)才能完全逼近原函數(shù)。但是實(shí)際應(yīng)用中,可以用有限項(xiàng)級(jí)數(shù)來近似表示,所選項(xiàng)數(shù)越多,級(jí)數(shù)就越逼近原函數(shù),當(dāng)N→∞時(shí),級(jí)數(shù)xN(t)應(yīng)趨近信號(hào)x(t)。但是,對(duì)于周期信號(hào)x(t),在x(t)的不連續(xù)處,x(t)的傅里葉級(jí)數(shù)在任何一點(diǎn)實(shí)際上都不等于x(t)。如x(t)在t＝t1處不連續(xù),那么其傅里葉級(jí)數(shù)在t－1和t＋1處都大約偏離9%。這種特性稱為吉布斯(Gibbs)現(xiàn)象。根據(jù)例4-1中的分析,可以寫出周期鋸齒信號(hào)的傅里葉級(jí)數(shù)表達(dá)式:

MATLAB程序如下:

t＝－3∶6/1000∶3；

N＝input(′級(jí)數(shù)諧波項(xiàng)數(shù)N的值:′)；

c0＝0；

w0＝pi；

xN＝c0*ones(1,length(t；%直流分量

fork＝1∶N,

ck＝j(luò)/k/pi*(－1)^k；

c_k＝－ck；

xN＝xN＋ck*exp(j*k*w0*t)＋c_k*exp(－j*k*w0*t)；%疊加諧波分量

end

plot(t,xN)

title([′N＝′,num2str(N)］)

xlabel(′Time(sec)′)

ylabel([′x′,num2str(N),′(t)′］)

N＝5時(shí),信號(hào)波形如圖4-4所示。

圖4-4N＝5時(shí)用級(jí)數(shù)表示的信號(hào)的波形

圖4-5N＝15時(shí)用級(jí)數(shù)表示的信號(hào)的波形

圖4-6N＝45時(shí)用級(jí)數(shù)表示的信號(hào)的波形

連續(xù)時(shí)間信號(hào)x(t)的傅里葉變換(ContinuousTimeFourierTransform,CTFT)定義為

(4.6)

4.3連續(xù)時(shí)間信號(hào)的傅里葉變換傅里葉反變換的定義為

(4.7)

一般情況下,傅里葉變換X(ω)是頻率變量ω的復(fù)數(shù)函數(shù)。為了確切說明X(ω),通常需要求出它的模函數(shù)|X(ω)|和輻角函數(shù)φ(ω)。模函數(shù)|X(ω)|隨ω變化的曲線稱為振幅頻譜；

輻角函數(shù)φ(ω)隨ω變化的曲線稱為相位頻譜。用MATLAB可以計(jì)算連續(xù)時(shí)間傅里葉變換積分的數(shù)值近似。先對(duì)連續(xù)時(shí)間信號(hào)進(jìn)行等間隔采樣,再利用fft函數(shù)計(jì)算連續(xù)時(shí)間傅里葉變換的近似值。根據(jù)傅里葉變換的定義,近似計(jì)算可以用公式表示為

(4.8)在τ足夠小的情況下,式(4.8)等號(hào)右邊的和式就是連續(xù)時(shí)間傅里葉變換很好的近似。若信號(hào)x(t)存在于0≤t<T的范圍內(nèi),那么式(4.8)可以寫成

(4.9)

其中,T＝Nτ,N為一個(gè)整數(shù)。用函數(shù)fft可以對(duì)一組離散的頻率ωk計(jì)算式(4.9)中的和式,如果信號(hào)x(t)的采樣信號(hào)x(nτ)的N個(gè)樣本存在向量x中,那么,向量X＝tau*fft(x)(式中tau代表符號(hào)τ)就表示用函數(shù)fft計(jì)算出式(4.9)中連續(xù)時(shí)間傅里葉變換X(ω)的N個(gè)樣本X(ωk),用向量表示為

(4.10)

其中,

(4.11)

設(shè)N為偶數(shù)。為了提高計(jì)算效率,函數(shù)fft先產(chǎn)生正頻率的樣本。要按頻率上升順序排列,可以使用函數(shù)fftshift。如要將存入向量X中的樣本排列成在 上的傅里葉變換的結(jié)果,就可以用命令X＝fftshift(tau*fft(x。例4-3求信號(hào)x(t)＝e－10tu(t)的連續(xù)時(shí)間傅里葉變換,畫出振幅頻譜,并且與用解析式法求得的頻譜進(jìn)行比較。

解：信號(hào)x(t)＝e－10tu(t)是無限長(zhǎng)的信號(hào),但是當(dāng)t>5時(shí)信號(hào)基本上為零,所以取信號(hào)范圍為0≤t<5,即T＝5。取τ＝0.01,則N＝T/τ＝500。

創(chuàng)建向量x存儲(chǔ)x(nτ)的樣本,對(duì)應(yīng)的時(shí)間向量樣本區(qū)間為t＝[0∶tau∶(T－tau)］；

創(chuàng)建向量X存儲(chǔ)X(ωk)的樣本,對(duì)應(yīng)的頻率向量樣本區(qū)間為w＝－(pi/tau)＋(0∶N－1)*(2*pi/(N*tau))。用解析式求得信號(hào)的頻譜函數(shù)為

可以用函數(shù)abs求出X(ω)的模。

MATLAB程序如下：

tau＝0.01；T＝5；N＝T/tau；

t＝[0∶tau∶(T－tau)］；x＝exp(－10*t)；

X＝fftshift(tau*fft(x；X＝abs(X)；

w＝－(pi/tau)＋(0∶N－1)*(2*pi/(N*tau；

Xe＝1./(10＋j*w)；Xe＝abs(Xe)；

plot(w,X,w,Xe,′ro′)

legend(′X-fft′,′Xe-解析式′)

xlabel(′Frequency(rad/sec)′)

ylabel(′|X|′)

結(jié)果如圖4-7所示。

圖4-7信號(hào)x(t)＝e－10t振幅頻譜圖圖中“—”表示的是用函數(shù)fft根據(jù)式(4.8)至式(4.11)求出的信號(hào)的振幅頻譜,“o”表示的是用解析式法求出的信號(hào)的振幅頻譜,二者基本一致。

離散時(shí)間傅里葉變換(DiscreteTimeFourierTransform,DTFT)也稱為序列的傅里葉變換,定義為

(4.12)

4.4離散時(shí)間信號(hào)的傅里葉變換其傅里葉反變換的定義為

(4.13)

X(ejω)是頻率變量ω的復(fù)數(shù)函數(shù),與連續(xù)時(shí)間傅里葉變換不一樣的是,離散時(shí)間傅里葉變換X(ejω)是一個(gè)周期為2π的周期函數(shù)。|X(ejω)|是它的模函數(shù),φ(ω)是它的輻角函數(shù)。模函數(shù)|X(ejω)|隨ω變化的曲線稱為振幅頻譜；輻角函數(shù)φ(ω)隨ω變化的曲線稱為相位頻譜。離散時(shí)間傅里葉變換X(ejω)是頻率ω的連續(xù)函數(shù),因?yàn)樗且粋€(gè)周期為2π的函數(shù),所以其振幅頻譜和相位頻譜曲線一般都畫在－π≤ω≤π或0≤ω≤2π的區(qū)間內(nèi)。

計(jì)算一個(gè)信號(hào)的離散時(shí)間傅里葉變換也可以用MATLAB中的fft函數(shù),但是在計(jì)算過程中存在兩個(gè)問題。首先,MATLAB的向量只能表示有限長(zhǎng)的信號(hào),如果x(n)是無限長(zhǎng)的,要將其截?cái)?變成一個(gè)有限長(zhǎng)的信號(hào)。第二,X(ejω)是定義在連續(xù)變量ω上的函數(shù),但是MATLAB只能在一組離散的頻率樣本上求值,不過如果將頻率樣本選得足夠多,那么由這些樣本得到的結(jié)論是真正DTFT很好的近似。為了使計(jì)算高效,最好選取在0≤ω<2π區(qū)間內(nèi)的一組頻率樣本,由ωk＝2πk/N,k＝0,1,…,N－1給出N個(gè)等分頻率點(diǎn)。對(duì)于一個(gè)存在于0≤n<M范圍內(nèi)的非零信號(hào)x(n),這組頻率樣本對(duì)應(yīng)的離散時(shí)間傅里葉變換為

(4.14)這個(gè)計(jì)算過程可以用函數(shù)fft實(shí)現(xiàn),如果信號(hào)x(n)的M個(gè)樣本存在向量x中,那么,向量X＝fft(x,N)(N≥M)就表示用函數(shù)fft計(jì)算出式(4.14)中離散時(shí)間傅里葉變換X(ejω)的N個(gè)樣本X(ejωk)。由于N<M的情況比較復(fù)雜,這里就不詳細(xì)介紹了。

由于離散時(shí)間傅里葉變換具有周期性,k≥N/2的向量X的值就是－π≤ω<0區(qū)間上的X(ejωk)的樣本值,因此,如果想在－π≤ωk<π區(qū)間畫出曲線,可以用fftshift函數(shù)將向量X的后一半樣本值與前一半樣本值交換。

用解析式求得信號(hào)的頻譜函數(shù)為例4-4求信號(hào)x(n)＝u(n)－u(n－5)的離散時(shí)間傅里葉變換,畫出振幅頻譜,并且與用解析式法求得的頻譜進(jìn)行比較。

解：創(chuàng)建向量x存儲(chǔ)x(n)的樣本,對(duì)應(yīng)的時(shí)間向量樣本區(qū)間為n＝[0∶4］,所以M＝5；

創(chuàng)建向量X存儲(chǔ)X(ejωk)的樣本,對(duì)應(yīng)的頻率向量樣本區(qū)間為w＝2*pi*[0∶N－1］/N,取N＝100。

可以用函數(shù)abs求出X(ejωk)的模。

MATLAB程序如下：

n＝0∶4；

x＝[ones(1,5)］；

w＝2*pi*[0∶99］/100；%0≤ω<2π

X＝fft(x,100)；

X1＝zeros(size(w；

M＝5；

den＝sin(w/2)；

num＝sin(M*w/2)；

X1(den~＝0)＝num(den~＝0)./den(den~＝0)；

X1(den＝＝0)＝M；

w＝w－pi；%－π≤ω<π

X＝fftshift(fft(x,100；

den＝sin(w/2)；

num＝sin(M*w/2)；

X1(den~＝0)＝num(den~＝0)./den(den~＝0)；

X1(den＝＝0)＝M；

figure

plot(w,abs(X),′r＋′,w,abs(X1

ylabel(′|X|′)

xlabel(′＼omega′)

gridon

legend(′X-fft′,′X1-解析式′)

程序中計(jì)算了0≤ω<2π和－π≤ω<π兩個(gè)區(qū)間的離散時(shí)間傅里葉變換，畫出的是在－π≤ω≤π區(qū)間的振幅頻譜,結(jié)果如圖4-8所示。

圖4-8信號(hào)x(n)＝u(n)－u(n－5)振幅頻譜圖圖中“＋”表示的是用函數(shù)fft根據(jù)式(4.14)求出的信號(hào)的振幅頻譜,“—”表示的是用解析式法求出的信號(hào)的振幅頻譜,二者基本一致。

4.5系統(tǒng)的頻域分析實(shí)例

系統(tǒng)的頻域分析主要研究傅里葉級(jí)數(shù)和傅里葉變換在線性時(shí)不變系統(tǒng)中的應(yīng)用。下面通過幾個(gè)實(shí)例介紹MATLAB在系統(tǒng)頻域分析中的應(yīng)用。例4-5RC電路的頻率響應(yīng)。

RC電路如圖4-9所示,電壓信號(hào)u(t)為輸入信號(hào)x(t),電容兩端的電壓信號(hào)uC(t)為輸出信號(hào)y(t),按下面的要求分析這個(gè)電路。

(1)求出系統(tǒng)的頻率響應(yīng)函數(shù)H(jω)的表達(dá)式,設(shè)RC＝0.001,畫出波形,分析系統(tǒng)的頻率特性；

(2)用函數(shù)freqs求系統(tǒng)的頻率響應(yīng),與(1)中用系統(tǒng)的頻率響應(yīng)函數(shù)表達(dá)式畫出的波形進(jìn)行比較；

(3)設(shè)輸入信號(hào)為x(t)＝cos(10t)＋cos(3000t),用頻域分析法求系統(tǒng)對(duì)x(t)的響應(yīng)y(t),并在－1≤t≤1區(qū)間內(nèi)畫出輸入、輸出信號(hào)的波形。

圖4-9RC電路解：

(1)系統(tǒng)的頻率響應(yīng)函數(shù)表達(dá)式為

MATLAB程序如下：

RC＝0.001；

w＝0∶50∶5000；

H＝(1/RC)./(j*w＋1/RC)；

magH＝abs(H)；

angH＝180*angle(H)/pi；

subplot(2,1,1),plot(w,magH)

subplot(2,1,2),plot(w,angH)；

畫出的波形如圖4-10所示。

圖4-10RC電路的頻率響應(yīng)圖4-10中,圖(a)為系統(tǒng)的幅頻特性,圖(b)為系統(tǒng)的相頻特性。

從系統(tǒng)的頻率響應(yīng)函數(shù)表達(dá)式可以求得:ω＝0時(shí),|H(ω)|＝1,φ(ω)＝0°；ω＝1000rad/s時(shí),|H(ω)|＝0.707,φ(ω)＝－45°；ω＝3000rad/s時(shí),|H(ω)|＝0.316,φ(ω)＝－71.6°；ω→∞時(shí),|H(ω)|→0,φ(ω)→－90°。所以,RC電路為一個(gè)一階低通濾波器,本題中RC電路的截止頻率為ω0＝1000rad/s。

(2)用函數(shù)freqs也可以求出系統(tǒng)的頻率響應(yīng),MATLAB程序如下：

w＝linspace(0,5000)；

a＝[11000］；

b＝1000；

H＝freqs(b,a,w)；

magH＝abs(H)；

angH＝180*angle(H)/pi；

subplot(2,1,1),plot(w,magH)；

subplot(2,1,2),plot(w,angH)；

畫出的波形如圖4-11所示。顯然,圖4-11和圖4-10的結(jié)果是一致的。

圖4-11用函數(shù)freqs得到的RC電路的頻率響應(yīng)

(3)系統(tǒng)的輸入信號(hào)是兩個(gè)正弦信號(hào)的和,一個(gè)信號(hào)的頻率較低,ω1＝10rad/s,一個(gè)信號(hào)的頻率較高,ω2＝3000rad/s,根據(jù)系統(tǒng)的頻率特性,系統(tǒng)的輸出可以表示為

當(dāng)信號(hào)通過低通濾波器時(shí),頻率較高的正弦信號(hào)會(huì)被衰減掉。MATLAB程序如下：

RC＝0.001；

t＝－1∶.2/100∶1；

w1＝10；w2＝3000；

Hw1＝(1/RC)/(j*w1＋1/RC)；

Hw2＝(1/RC)/(j*w2＋1/RC)；

x＝cos(w1*t)＋cos(w2*t)；

y＝abs(Hw1)*cos(w1*t＋angle(Hw1＋abs(Hw2)*cos(w2*t＋angle(Hw2；

subplot(2,1,1),plot(t,x)

ylabel(′x(t)′)

xlabel(′Time(sec)′)

subplot(2,1,2),plot(t,y)；

xlabel(′Time(sec)′)

ylabel(′y(t)′)

畫出的波形如圖4-12所示。

圖4-12輸入波形和輸出波形圖4-12中,圖(a)為輸入信號(hào)的波形,圖(b)為輸出信號(hào)的波形。從圖中可以看出,兩圖中低頻信號(hào)的振幅幾乎相同,而輸出信號(hào)中高頻信號(hào)的幅值有了一定的衰減,大約為原幅度值的0.316倍。例4-6破譯密碼(幅度調(diào)制)。

幅度調(diào)制是一種常用的信號(hào)調(diào)制方法,一種簡(jiǎn)單的幅度調(diào)制可以表示為

其中,m(t)為調(diào)制信號(hào),cos(ω0t)為載波信號(hào),f0是載波的頻率,s(t)為已調(diào)信號(hào)。根據(jù)傅里葉變換的性質(zhì),容易求出已調(diào)信號(hào)的頻譜S(jω)為

其中,M(jω)為調(diào)制信號(hào)m(t)的傅里葉變換。從已調(diào)信號(hào)的頻譜可以看出,幅度調(diào)制把調(diào)制信號(hào)的頻譜搬移到載波的頻率附近。

假設(shè)有一個(gè)代碼被編寫在信號(hào)x(t)中,它是“技術(shù)的未來在于…”這句話的最后幾個(gè)字。破譯這個(gè)代碼所需的信息包含在MATLAB文件ctftmod.m中,文件中包括一個(gè)低通濾波器的頻率響應(yīng)的參數(shù)、典型的“dot”信號(hào)和“dash”信號(hào)的樣本、載波的頻率f1和f2,還有信號(hào)x(t)的樣本。信號(hào)x(t)表示為

信號(hào)m1(t)、m2(t)和m3(t)對(duì)應(yīng)于字母表中的一個(gè)字母,字母表用國(guó)際莫爾斯碼進(jìn)行編碼,編碼見表4-1。表4-1莫爾斯碼編碼表

每個(gè)字母都是由“dot”信號(hào)和“dash”信號(hào)組成的,例如,用莫爾斯碼組成對(duì)應(yīng)于字母“X”的信號(hào)為“-··-”,可以用向量x表示為

x＝[dashdotdotdash］

信號(hào)“dot”和“dash”都是由低頻分量組成的,它們的傅里葉變換的頻率范圍都大致位于給出的低通濾波器的通帶之內(nèi)。當(dāng)信號(hào)“dot”和“dash”被cos(2πf1t)調(diào)制時(shí),傅里葉變換中的大部分能量都會(huì)移到該濾波器的通帶以外。試根據(jù)以上說明解決下面的問題。

(1)用命令whos得到MATLAB文件ctftmod.m中的內(nèi)容。

(2)利用信號(hào)“dot”和“dash”,以莫爾斯碼組成對(duì)應(yīng)于字母“Z”的信號(hào),并對(duì)t畫出這個(gè)信號(hào)的波形。

(3)用freqs畫出文件中給出的濾波器的頻率響應(yīng)。

(4)用lsim求“dot”和“dash”信號(hào)通過濾波器的響應(yīng),畫出原信號(hào)和輸出信號(hào)的波形,說明這兩個(gè)信號(hào)都是由低頻分量組成的。

(5)當(dāng)“dash”信號(hào)被cos(2πf1t)調(diào)制以后,用信號(hào)y(t)表示已調(diào)信號(hào),用lsim求y(t)通過濾波器的響應(yīng),畫出信號(hào)y(t)和輸出信號(hào)yo(t)的波形,說明結(jié)果。

(6)用解析法求下列信號(hào)的傅里葉變換：

m1(t)cos(2πf1t)cos(2πf1t)

m2(t)cos(2πf1t)sin(2πf2t)

m3(t)cos(2πf1t)sin(2πf1t)

(7)利用(5)、(6)中的結(jié)論,設(shè)計(jì)一個(gè)從信號(hào)x(t)中提取信號(hào)m1(t)的方案,畫出信號(hào)m1(t)的波形,并確定該信號(hào)代表的是哪個(gè)字母的莫爾斯碼。

(8)根據(jù)(7)的結(jié)果,分別從信號(hào)x(t)中提取信號(hào)m2(t)和m3(t),畫出m2(t)和m3(t)的波形,并確定這兩個(gè)信號(hào)代表的是哪些字母的莫爾斯碼。

(9)如果你是一個(gè)破譯密碼的偵探，你能根據(jù)以上的介紹得出技術(shù)的未來在于什么嗎？

解：

(1)在MATLAB命令窗口輸入以下命令,可以得到結(jié)果：

clear

loadctftmod

whos

NameSizeBytesClass

af1x648doublearray

bf1x648doublearray

dash1x200016000doublearray

dot1x200016000doublearray

f11x18doublearray

f21x18doublearray

t1x800064000doublearray

x1x800064000doublearray

這些就是文件ctftmod.m中的內(nèi)容。其中af和bf是一個(gè)低通濾波器的頻率響應(yīng)的參數(shù)：

af＝1.0e＋009*[0.00000.00000.00000.00220.10242.3730］

bf＝1.0e＋009*[00.00000.00000.00000.00002.3730］

f1和f2是載波的頻率：

f1＝200

f2＝400

用plot命令畫出典型的“dash”信號(hào)和“dot”信號(hào)的波形,如圖4-13所示。

圖4-13典型的“dash”信號(hào)和“dot”信號(hào)的波形

(2)字母“Z”的莫爾斯碼組成為“--··”,用向量z表示為：z＝[dashdashdotdot］,用命令plot(z)畫出“Z”信號(hào)的波形如圖4-14所示。

(3)用命令freqs(bf,af)可以畫出濾波器在對(duì)數(shù)坐標(biāo)系下的頻率響應(yīng)如圖4-15所示。

圖4-14字母“Z”信號(hào)的波形

圖4-15給定濾波器的頻率響應(yīng)從濾波器的頻率響應(yīng)可以看出,這是一個(gè)低通濾波器。

(4)用lsim求“dot”和“dash”信號(hào)通過濾波器的響應(yīng)的MATLAB程序如下：

figure

ydash＝lsim(bf,af,dash,t(1∶length(dash)；

ydot＝lsim(bf,af,dot,t(1∶length(dot)；

subplot(2,1,1)；plot(t(1∶length(dash,dash,t(1∶length(dash,ydash,′--′)

legend(′dash′,′ydash′)

subplot(2,1,2)；plot(t(1∶length(dot,dot,t(1∶length(dot,ydot,′--′)

legend(′dot′,′ydot′)

運(yùn)行結(jié)果如圖4-16所示。

圖4-16原信號(hào)和輸出信號(hào)的波形圖4-16中,圖(a)為“dash”信號(hào)和通過濾波器的輸出信號(hào)“ydash”；圖(b)為“dot”信號(hào)和輸出信號(hào)“ydot”。從信號(hào)波形上可以看出,原信號(hào)和輸出信號(hào)的波形基本一致,輸出信號(hào)在幅度上有所改變,在時(shí)間上有一些延遲。這個(gè)結(jié)果說明這兩個(gè)信號(hào)都是由低頻分量組成的,它們的傅里葉變換的頻率范圍都大致位于給出的低通濾波器的通帶之內(nèi)。

(5)用lsim求已調(diào)信號(hào)y(t)通過濾波器的響應(yīng)的MATLAB程序?yàn)椋?/p>

figure

y＝dash.*cos(2*pi*f1*t(1∶length(dash)；

yo＝lsim(bf,af,y,t(1∶length(dash)；

plot(t(1∶length(dash,y,t(1∶length(dash,yo,′r*′)；

legend(′y＝dash*cos(2pi*f1*t)′,′yo-systemoutputofy′)

運(yùn)行結(jié)果如圖4-17所示。

圖4-17已調(diào)信號(hào)y(t)和輸出信號(hào)yo(t)的波形圖4-17中輸出信號(hào)yo(t)幾乎為零,這說明“dash”信號(hào)被cos(2πf1t)調(diào)制后,信號(hào)的頻譜基本被搬移到低通濾波器的通帶之外,所以已調(diào)信號(hào)y(t)通過濾波器后,幾乎所有的頻率分量都被衰減,輸出信號(hào)yo(t)幾乎為零。

(6)用解析法求信號(hào)的傅里葉變換分別為：(7)把(6)中的三個(gè)表達(dá)式相加,得到的是信號(hào)x(t)被cos(2πf1t)調(diào)制的結(jié)果,表達(dá)式為

該已調(diào)信號(hào)的頻譜為

已知f1＝200,f2＝400,那么,2π(f2－f1)＝400π,2π(f2＋f1)＝1200π,4πf1＝800π。從已調(diào)信號(hào)的頻譜可知,如果要從信號(hào)x(t)中提取信號(hào)m1(t),可以先用cos(2πf1t)調(diào)制信號(hào)x(t),然后讓已調(diào)信號(hào)通過一個(gè)低通濾波器,得到M1(ω),即可得到m1(t)。

圖4-18信號(hào)m1(t)的波形考察題目給出的低通濾波器,其截止頻率滿足要求,所以,只要讓已調(diào)信號(hào)通過這個(gè)低通濾波器,得到的輸出信號(hào)就是與m1(t)波形相同的信號(hào)。

MATLAB程序如下：

figure

m1＝lsim(bf,af,x.*cos(2*pi*f1*t),t)；

plot(t,m1)

輸出波形如圖4-18所示。比對(duì)輸出波形和字母的莫爾斯碼信號(hào)的波形,可以判斷出這是字母D(-··)。

(8)根據(jù)(7)的結(jié)果,要從信號(hào)x(t)中提取信號(hào)m2(t),可以先用sin(2πf2t)調(diào)制信號(hào)x(t),然后讓已調(diào)信號(hào)通過給定低通濾波器,得到M2(ω),