（內容源自詳解MATLAB／SIMULINK 通信系統建模與仿真?? 劉學勇編著第三章內容，有興趣的讀者請閱讀原書）

一、離散傅里葉變換

?

?

clear all
n=0:30;%信號的時間范圍
x=sin(0.2*n).*exp(-0.1*n);
k=0:30;%頻率范圍
N=31;
Wnk=exp(-j*2*pi/N).^(n'*k);%DFT公式
X=x*Wnk;
subplot(2,1,1);stem(n,x);title('序列x');
subplot(2,1,2);stem(-15:15,[abs(X(17:end)) abs(X(1:16))]);%將DFT下標重新排列（DFT的默認下標為[0,N-1]）
title('X幅度');

二、循環卷積

clear all
h=[6 3 4 2 1 -2];
x=[3 2 6  7 -1 -3];
h1=fliplr(h)%反轉序列h
H=toeplitz(h,[h(1) h1(1:5)]);%生成循環矩陣
%T = toeplitz(c,r) 返回非對稱托普利茨矩陣，
%其中 c 作為第一列，r 作為第一行。如果 c 和 r 的首個元素不同，toeplitz 將發出警告并使用列元素作為對角線
y=H*x';%計算循環卷積序列H=fft(h);
X=fft(x);%兩個序列的DFT
Y=H.*X;%循環卷積序列的DFT
y1=ifft(Y);%求循環卷積subplot(2,1,1);stem(y);title('直接計算');
subplot(2,1,2);stem(y1);title('DFT計算');

?

這里用兩種方得到的結果相同證明DFT運算內部是采用循環卷積而不是線性卷積

三、利用DFT求線性卷積

?

clear all
n1=0:20;
n2=0:10;
h=sinc(0.2*n1);
x=exp(-0.2*n2);
y=conv(x,h);%1.直接時域卷積h1=[h zeros(1,length(x)-1)];
x1=[x zeros(1,length(h)-1)];%對x(n),h(n)進行擴展
H1=fft(h1);
X1=fft(x1);%計算DFT
Y1=H1.*X1;
y1=ifft(Y1);%2.利用IDFT算出時域卷積結果
subplot(2,1,1);stem(y);title('直接計算');
subplot(2,1,2);stem(y1);title('DFT計算');

?

雖然DFT內部采用的是循環卷積，但是我們可以利用DFT實現線性卷積結果的求解

四、Hilbert變換

clear all
ts=0.01;%這里是連續信號，所以需要對其進行采樣才能處理，因為載波頻率為20HZ，由采樣定理至少需要40HZ，這里采用0.01s作為時間間隔
fs=1/ts;
t=0:ts:10;
df=fs/length(t);%DFT的頻譜分辨率
f=-50:df:50-df;%生成頻率矢量（用于畫圖）
x=exp(-10*abs(t-5)).*cos(2*pi*20*t);
X=fft(x)/fs;%求信號的頻譜，因為原始信號是模擬信號，由抽樣定理，必須求出FFT后再除以fs才是原模擬信號的頻譜xa=hilbert(x);%求解析信號
Xa=fft(xa)/fs;%求解析信號的頻譜
subplot(2,1,1);plot(t,x);%利用時間分辨率畫圖
title('信號x');xlabel('時間t');
subplot(2,1,2);plot(f,fftshift(abs(X)));title('信號x幅度譜');xlabel('頻率f');
%fftshift將0頻分量移到中心，重新排列abs(X),由于X是復數，所以利用abs求包絡
figure
subplot(2,1,1);plot(t,abs(xa));title('信號xa包絡');xlabel('時間t');
subplot(2,1,2);plot(f,fftshift(abs(Xa)));title('信號xa幅度譜');xlabel('頻率f');

?可以看到解析信號只包含正頻率部分，且幅度值是原始信號頻譜幅度值的兩倍

五、帶通信號的低通表示

clear all
ts=0.01;%這里是連續信號，所以需要對其進行采樣才能處理，因為載波頻率為20HZ，由采樣定理至少需要40HZ，這里采用0.01s作為時間間隔
fs=1/ts;
t=0:ts:10;
df=fs/length(t);%DFT的頻譜分辨率
f=-50:df:50-df;%生成頻率矢量（用于畫圖）
x=exp(-10*abs(t-5)).*cos(2*pi*20*t);
X=fft(x)/fs;%求信號的頻譜，因為原始信號是模擬信號，由抽樣定理，必須求出FFT后再除以fs才是原模擬信號的頻譜xa=hilbert(x);%求解析信號fc1=20;
x11=xa.*exp(-j*2*pi*fc1*t);%3-41公式，求20hz的低通信號
X11=fft(x11)/fs;%模擬信號，fft之后要除以fs
subplot(2,1,1);plot(t,real(x11));title('fc=20HZ時的低通信號同相分量');xlabel('時間t');
%同相分量就是解析信號中的實部
subplot(2,1,2);plot(f,fftshift(abs(X11)));title('fc=20HZ時的低通信號幅度譜');xlabel('頻率f');fc2=10;
x12=xa.*exp(-j*2*pi*fc2*t);
X12=fft(x12)/fs;
figure
subplot(2,1,1);plot(t,real(x12));title('fc=10HZ時的低通信號同相分量');xlabel('時間t');
subplot(2,1,2);plot(f,fftshift(abs(X12)));title('fc=10HZ時的低通信號幅度譜');xlabel('頻率f');

總結來說，就是利用解析信號和頻移公式實現了帶通信號的低通表示。

六、平穩隨機過程

clear all
N1=2000;
N2=100;
x=randn(N2,N1);%產生一個100行，2000列的高斯分布隨機數
for ii=1:N2%對每行進行循環[Rx(ii,:),lags]=xcorr(x(ii,:),50,'coeff');%xcorr是用來計算自相關函數，第一個參數是自相關函數序列，%第二個參數50是自相關值的最大時間偏移，lags就是[-50:50],可以理解為x中的每一行有2000個數，經過xcorr后每一行轉化為101個數的%自相關函數Sf(ii,:)=fftshift(abs(fft(Rx(ii,:))));%功率譜密度是自相關函數的傅里葉變換
endRx_av=sum(Rx)/N2;%求100行的平均值
Sf_av=sum(Sf)/N2;subplot(2,1,1);plot(lags,Rx_av);title('自相關函數')%畫出在最大時間偏移上的自相關函數的變化情況
subplot(2,1,2);plot(lags,Sf_av);title('功率譜密度')
axis([-50 50 0 2]);

?

七、帶通隨機過程

clear all
ts=0.002;
tao=-1:ts:1;%時間序列，因為自相關函數是連續函數，所以計算時需要抽樣，因為函數是sinc形式，所以在【-1，1】之外值會很小
%可以忽略不計
B=20;
f0=100;
R=sinc(2*B*tao).*cos(2*pi*f0*tao);fs=1/ts;
df=fs/length(tao);%頻域分辨率，用采樣頻率/時間序列的長度
f=-fs/2:df:fs/2-df;
S=fft(R)/fs;%因為這里自相關函數是連續函數，所以求功率譜密度時還要在除以fs
%類似于之前的例子中模擬信號求fft之后還要再除以fs才是正確地結果
subplot(2,1,1);plot(tao,R);title('自相關函數');xlabel('tao');ylabel('R');
subplot(2,1,2);plot(f,fftshift(abs(S)));title('功率譜密度');xlabel('f');ylabel('S');

八、隨機過程通過線性系統

clear all
N1=2000;
N2=100;
x=randn(N2,N1);%相同于例3.17
for ii=1:N2%對每一行進行循環y(ii,1)=x(ii,1);%由濾波器的響應可以推出來y(n)=0.6*y(n-1)+x(n),y(-1)=0,所以y(0)=x(0);for jj=2:N1y(ii,jj)=0.6*y(ii,jj-1)+x(ii,jj);%遞推公式end[Ry(ii,:),lags]=xcorr( y(ii,:),50,'coeff');Sf(ii,:)=fftshift(abs(fft(Ry(ii,:))));
endRy_av=sum(Ry)/N2;%平均化
Sf_av=sum(Sf)/N2;
subplot(2,1,1);plot(lags,Ry_av);title('自相關函數');
subplot(2,1,2);plot(lags,Sf_av);title('功率譜密度');