多徑傳輸的源程序： 
(一)，用到的函數

1,f2t函數

   function x=f2t(X)

   global dt df t f T N

   %x=f2t(X)

   %x為時域的取樣值矢量

   %X為x的傅氏變換

   %X與x長度相同并為2的整冪

   %本函數需要一個全局變量dt(時域取樣間隔)

   X=[X(N/2+1:N),X(1:N/2)];

   x=ifft(X)/dt;

   end

 

2,t2f函數。

   function X=t2f(x)

   global dt df N t f T

   %X=t2f(x)

   %x為時域的取樣值矢量

   %X為x的傅氏變換

   %X與x長度相同,并為2的整冪。

   %本函數需要一個全局變量dt(時域取樣間隔)

   H=fft(x);

   X=[H(N/2+1:N),H(1:N/2)]*dt;

   end

 

(二),主程序

   global dt t f df N T

   close all

   clear Eb_N0 Pe

   

   k=input('取樣點數＝2^k, k取13左右');

   if isempty(k), k=13; end

   z=input('每個信號取樣點數＝2^z, z<k');

   if isempty(z), z=3; end

   aa=input('滾降系數=[0.5]');

   if aa==[],aa=0.5;end

   N=2^k;

   L=2^z;M=N/L;

   Rb=2;                     %碼速率是2Mb/s

   Ts=1/Rb;                   %碼元間隔

   dt=Ts/L;                   %時域采樣間隔

   df=1/(N*dt);               %頻域采樣間隔

   T=N*dt;                    %截短時間

   Bs=N*df/2;                 %系統帶寬

   f=[-Bs+df/2:df:Bs];        %頻域橫坐標

   t=[-T/2+dt/2:dt:T/2];      %時域橫坐標

   g=sin(pi*t/Ts).*cos(pi*t*aa/Ts)./[pi*t/Ts.*(1-4*t.^2*aa^2/Ts^2)]; %升余弦脈沖波形

   GG=t2f(g);GG=abs(GG);      %升余弦脈沖的傅式變換

   GT=sqrt(GG);GR=GT;         %最佳系統的發送接收濾波器的傅式變換

   

   for l1=1:20;

   Eb_N0(l1)=(l1-1);          %Eb/N0 in dB

   eb_n0(l1)=10^(Eb_N0(l1)/10);

   Eb=1;  

   n0=Eb/eb_n0(l1);           %信道的噪聲譜密度  

   sita=n0*Bs;               %信道中噪聲功率

   

   n_err=0;                  %誤碼計數

   for l2=1:5; 

   dm=Ts/2;

   C=abs(1-0.5*exp(-j*2*pi*f*dm));%多徑信道

   b=sign(randn(1,M));

   s=zeros(1,N);            %產生沖激序列

   s(L/2:L:N)=b/dt;

   SS=t2f(s);

   S=SS.*GG.*C;             %信道的傅式變化

   a=f2t(S);a=real(a);      %不加噪聲的輸出

   n_ch=sqrt(sita)*randn(size(t)); %信道噪聲

   N_CH=t2f(n_ch);

   nr=real(f2t(N_CH.*GR));  %輸出噪聲

   sr=a+nr;                 %接收信號

   y=sr(L/2:L:N);           %取樣

   bb=sign(y);              %判決

   n_err=n_err+length(find(bb~=b));%錯誤累計 

   end

   Pe(l1)=n_err/(M*l2);

   semilogy(Eb_N0,Pe,'black');   %Pe~Eb/N0曲線畫圖

   xlabel('Eb/N0');ylabel('Pe');title('Pe~Eb/N0曲線');

   eb_n0=10.^(Eb_N0/10);

   hold on

   semilogy(Eb_N0,0.5*erfc(sqrt(eb_n0)));

   axis([0,15,1e-3,1])

   xlabel('Eb/N0')

   ylabel('Pe')

   legend('實際的','理論的')

   end