matlab编写2psk、qpsk、8psk的调制解调及其蒙特卡罗仿真

以下是MATLAB代码示例，实现了2PSK、QPSK、8PSK的调制解调及其蒙特卡罗仿真。

2PSK调制解调：

调制：

clc; clear all; close all;
M = 2; %调制阶数
N = 1000; %符号数
fc = 10^6; %载波频率
fs = 10*fc; %采样频率
T = 1/fc; %符号周期
t = [0:1/fs:N*T-1/fs]; %时间序列
f = 2*fc; %信号频率
phi = 0; %相位
a = rand(1,N)>0.5; %随机产生二进制数字
s = sqrt(2/T)*a-1; %基带信号
ss = s.*cos(2*pi*f*t+phi); %调制信号
figure(1);
subplot(2,1,1);plot(t,s);title('基带信号');xlabel('时间');
subplot(2,1,2);plot(t,ss);title('调制信号');xlabel('时间');

解调：

ss = ss.*cos(2*pi*f*t+phi); %解调信号
[b,a] = butter(5,fc/(fs/2)); %设计滤波器
s_hat = filter(b,a,ss); %滤波
figure(2);
subplot(2,1,1);plot(t,ss);title('解调信号');xlabel('时间');
subplot(2,1,2);plot(t,s_hat);title('解调后信号');xlabel('时间');

2PSK蒙特卡罗仿真：

clc; clear all; close all;
M = 2; %调制阶数
N = 1000; %符号数
fc = 10^6; %载波频率
fs = 10*fc; %采样频率
T = 1/fc; %符号周期
f = 2*fc; %信号频率
phi = 0; %相位
EbN0_dB = [0:10]; %信噪比范围
for ii = 1:length(EbN0_dB)
    EbN0 = 10^(EbN0_dB(ii)/10); %线性信噪比
    a = rand(1,N)>0.5; %随机产生二进制数字
    s = sqrt(2/T)*a-1; %基带信号
    ss = s.*cos(2*pi*f*t+phi); %调制信号
    %加入高斯白噪声
    P = mean(abs(ss).^2); %信号功率
    N0 = P/EbN0; %噪声功率
    noise = sqrt(N0/2)*(randn(size(t))+j*randn(size(t))); %高斯白噪声
    r = ss + noise; %接收信号
    %解调
    r = r.*cos(2*pi*f*t+phi); %解调信号
    [b,a] = butter(5,fc/(fs/2)); %设计滤波器
    s_hat = filter(b,a,r); %滤波
    %判决
    a_hat = s_hat>0;
    %计算误码率
    err(ii) = sum(a~=a_hat);
end
%理论误码率
theory_err = 0.5*erfc(sqrt(EbN0));
%画图
semilogy(EbN0_dB,err/N,'ro-',EbN0_dB,theory_err,'b*-');
xlabel('Eb/N0(dB)');ylabel('误码率');grid;
legend('实际误码率','理论误码率');

QPSK调制解调：

调制：

clc; clear all; close all;
M = 4; %调制阶数
N = 1000; %符号数
fc = 10^6; %载波频率
fs = 10*fc; %采样频率
T = 1/fc; %符号周期
t = [0:1/fs:N*T-1/fs]; %时间序列
f = 2*fc; %信号频率
phi = 0; %相位
a = randi([0,3],1,N); %随机产生4进制数字
s = (1-2*a(1:2:end)) + j*(1-2*a(2:2:end)); %基带信号
ss = s.*exp(j*2*pi*f*t); %调制信号
figure(1);
subplot(2,1,1);plot(t,real(s));title('基带实部');xlabel('时间');
subplot(2,1,2);plot(t,imag(s));title('基带虚部');xlabel('时间');
figure(2);
subplot(2,1,1);plot(t,real(ss));title('调制实部');xlabel('时间');
subplot(2,1,2);plot(t,imag(ss));title('调制虚部');xlabel('时间');

解调：

ss = ss.*exp(-j*2*pi*f*t); %解调信号
[b,a] = butter(5,fc/(fs/2)); %设计滤波器
s_hat = filter(b,a,ss); %滤波
figure(3);
subplot(2,1,1);plot(t,real(ss));title('解调实部');xlabel('时间');
subplot(2,1,2);plot(t,imag(ss));title('解调虚部');xlabel('时间');
figure(4);
subplot(2,1,1);plot(t,real(s_hat));title('解调后实部');xlabel('时间');
subplot(2,1,2);plot(t,imag(s_hat));title('解调后虚部');xlabel('时间');

QPSK蒙特卡罗仿真：

clc; clear all; close all;
M = 4; %调制阶数
N = 1000; %符号数
fc = 10^6; %载波频率
fs = 10*fc; %采样频率
T = 1/fc; %符号周期
f = 2*fc; %信号频率
phi = 0; %相位
EbN0_dB = [0:10]; %信噪比范围
for ii = 1:length(EbN0_dB)
    EbN0 = 10^(EbN0_dB(ii)/10); %线性信噪比
    a = randi([0,3],1,N); %随机产生4进制数字
    s = (1-2*a(1:2:end)) + j*(1-2*a(2:2:end)); %基带信号
    ss = s.*exp(j*2*pi*f*t); %调制信号
    %加入高斯白噪声
    P = mean(abs(ss).^2); %信号功率
    N0 = P/EbN0; %噪声功率
    noise = sqrt(N0/2)*(randn(size(t))+j*randn(size(t))); %高斯白噪声
    r = ss + noise; %接收信号
    %解调
    r = r.*exp(-j*2*pi*f*t); %解调信号
    [b,a] = butter(5,fc/(fs/2)); %设计滤波器
    s_hat = filter(b,a,r); %滤波
    %判决
    a_hat = [real(s_hat)>0; imag(s_hat)>0];
    a_hat = a_hat(:)';
    a = (dec2bin(a,2)-'0')';
    a = a(:)';
    %计算误码率
    err(ii) = sum(a~=a_hat);
end
%理论误码率
theory_err = 0.5*erfc(sqrt(EbN0));
%画图
semilogy(EbN0_dB,err/N,'ro-',EbN0_dB,theory_err,'b*-');
xlabel('Eb/N0(dB)');ylabel('误码率');grid;
legend('实际误码率','理论误码率');

8PSK调制解调：

调制：

clc; clear all; close all;
M = 8; %调制阶数
N = 1000; %符号数
fc = 10^6; %载波频率
fs = 10*fc; %采样频率
T = 1/fc; %符号周期
t = [0:1/fs:N*T-1/fs]; %时间序列
f = 2*fc; %信号频率
phi = 0; %相位
a = randi([0,7],1,N); %随机产生8进制数字
s = exp(j*2*pi/8*(a)); %基带信号
ss = s.*exp(j*2*pi*f*t); %调制信号
figure(1);
subplot(2,1,1);plot(t,real(s));title('基带实部');xlabel('时间');
subplot(2,1,2);plot(t,imag(s));title('基带虚部');xlabel('时间');
figure(2);
subplot(2,1,1);plot(t,real(ss));title('调制实部');xlabel('时间');
subplot(2,1,2);plot(t,imag(ss));title('调制虚部');xlabel('时间');

解调：

ss = ss.*exp(-j*2*pi*f*t); %解调信号
[b,a] = butter(5,fc/(fs/2)); %设计滤波器
s_hat = filter(b,a,ss); %滤波
figure(3);
subplot(2,1,1);plot(t,real(ss));title('解调实部');xlabel('时间');
subplot(2,1,2);plot(t,imag(ss));title('解调虚部');xlabel('时间');
figure(4);
subplot(2,1,1);plot(t,real(s_hat));title('解调后实部');xlabel('时间');
subplot(2,1,2);plot(t,imag(s_hat));title('解调后虚部');xlabel('时间');

8PSK蒙特卡罗仿真：

clc; clear all; close all;
M = 8; %调制阶数
N = 1000; %符号数
fc = 10^6; %载波频率
fs = 10*fc; %采样频率
T = 1/fc; %符号周期
f = 2*fc; %信号频率
phi = 0; %相位
EbN0_dB = [0:10]; %信噪比范围
for ii = 1:length(EbN0_dB)
    EbN0 = 10^(EbN0_dB(ii)/10); %线性信噪比
    a = randi([0,7],1,N); %随机产生8进制数字
    s = exp(j*2*pi/8*(a)); %基带信号
    ss = s.*exp(j*2*pi*f*t); %调制信号
    %加入高斯白噪声
    P = mean(abs(ss).^2); %信号功率
    N0 = P/EbN0; %噪声功率
    noise = sqrt(N0/2)*(randn(size(t))+j*randn(size(t))); %高斯白噪声
    r = ss + noise; %接收信号
    %解调
    r = r.*exp(-j*2*pi*f*t); %解调信号
    [b,a] = butter(5,fc/(fs/2)); %设计滤波器
    s_hat = filter(b,a,r); %滤波
    %判决
    a_hat = angle(s_hat);
    a_hat = round(a_hat/(2*pi/8));
    a_hat = dec2bin(a_hat,3)-'0';
    a_hat = a_hat(:)';
    a = (dec2bin(a,3)-'0')';
    a = a(:)';
    %计算误码率
    err(ii) = sum(a~=a_hat);
end
%理论误码率
theory_err = 2/3*erfc(sqrt(3/2*EbN0)*sin(pi/8));
%画图
semilogy(EbN0_dB,err/N,'ro-',EbN0_dB,theory_err,'b*-');
xlabel('Eb/N0(dB)');ylabel('误码率');grid;
legend('实际误码率','理论误码率');