主程序如下：
Fd=1;
Fs=3*Fd;
Fd=1;
Fs=3*Fd;
M=1;
M=16;
SNR_db=0:12;
SNR1_db=0:0.1:12;
%誤碼情況仿真
for n=1:length(SNR_db)
Eb_N0=10^(SNR_db(n)/10);
sgma=sqrt(1/(8*Eb_N0));
x=randint(10000,1,M);
y=dmodce(x,Fd,Fs,'psk',M);
ynoise=y+sqrt(Fs/Fd)*sgma*(randn(length(y),1)+j*randn(length(y),1));%加噪
z=ddemodce(ynoise,Fd,Fs,'psk',M);
[number,pm(n)]=symerr(x,z);
end;
for m=1:length(SNR1_db)
Eb_N0=10^(SNR1_db(m)/10);
pm_theory(m)=2*Qfunct(sqrt(2*log2(M)*Eb_N0)*sin(pi/M));%理論誤符號率
pe_theory(m)=pm_theory(m)/log2(M);
end;
semilogy(SNR_db,pm,'*',SNR1_db,pm_theory);%作圖
xlabel('信噪比(dB)');
ylabel('誤符號率');

Qfunct函數的.M文件：
function [y]=Qfunct(x)
%[y]=Qfunct(x)
%Q函數
%y=1/sqrt(2*pi)*integral from x to inf of exp(-t^2/2)dt.
y=(1/2)*erfc(x/sqrt(2));