通信原理：QPSK频带调制理论及matlab仿真实验（包含多普勒频偏）

最新推荐文章于 2025-04-02 09:06:44 发布

黎明所至

最新推荐文章于 2025-04-02 09:06:44 发布

阅读量632

点赞数 15

文章标签： matlab 开发语言

本文链接：https://blog.youkuaiyun.com/liming_dawn_h/article/details/145372705

版权

QPSK概念

QPSK基本原理

对应相位值：

QPSK为两路正交的2PSK

调制：

解调：

科斯塔斯环

科斯塔斯环(Costas)法又称同相正交环法或边环法。它仍然利用锁相环提取载频，但是不需要对接收信号作平方运算就能得到载频输出。误差信号是由两路相乘及低通滤波器提供的。

代码模块（不含多普勒频偏）

% 不包含多普勒频偏的仿真
clear;
clc;
close all
format long

%% 参数设置
N = 300;       %比特数
T = 10^-3;        %比特周期
rolloff = 0.25; % 滚降系数
span = 8; % 滤波器的跨度
sps = 20; % 每个符号的采样点数
fc = 2000;       %载波频率
Rb=1/T;         %比特速率
Rs=Rb/2;        %符号速率
Fs = Rb*sps;     %抽样频率
t=0:(1/Fs):(N*T-1/Fs); %时间
SNR = 5; %信噪比

%% 信号发生
signal = randi([0,1], 1, N); %随机产生的比特数0、1
signal = 2*signal-1;      %单极性变为双极性-1、1

%% 串/并转换
Isignal=[]; Qsignal=[];
for i=1:N
    if mod(i, 2)~=0
        Isignal = [Isignal, signal(i)];
    else
        Qsignal = [Qsignal, signal(i)];
    end
end

%% 绘图比较I、Q比特流
bit_data=[];
for i=1:N
    bit_data = [bit_data, signal(i)*ones(1, T*Fs)]; 
end
I_data=[];
Q_data=[];
for i=1:N/2
    I_data = [I_data, Isignal(i)*ones(1, T*Fs)];
    Q_data = [Q_data, Qsignal(i)*ones(1, T*Fs)];
end

figure
subplot(3, 1, 1)
plot(0:(1/Fs):(N*T/2-1/Fs),I_data);  
legend('Isignal')  
xlabel('时间t');
ylabel('幅值');
subplot(3, 1, 2)
plot(0:(1/Fs):(N*T/2-1/Fs), Q_data);  
legend('Qsignal')   
xlabel('时间t');
ylabel('幅值');
subplot(3, 1, 3)
plot(0:(1/Fs):(N*T-1/Fs), bit_data); 
legend(' QPSK signal')  
xlabel('时间t');
ylabel('幅值');

%% 脉冲成型
% 平方根升余弦FIR滤波器
h = rcosdesign(rolloff, span, sps);

% 上采样
Isignal_upsampled = upsample(Isignal , sps);
Qsignal_upsampled = upsample(Qsignal , sps);

% 脉冲整形
Isignal_filtered = conv(Isignal_upsampled, h,'same');
Qsignal_filtered = conv(Qsignal_upsampled, h,'same');

%% 正交调制
% 载波信号
t2=0:(1/Fs):(N*T/2-1/Fs);
I_carrier = 2*Isignal_filtered.*cos(2*pi*fc*t2);       %I路载波信号
Q_carrier = -2*Qsignal_filtered.*sin(2*pi*fc*t2);  %Q路载波信号
QPSK_signal = I_carrier + Q_carrier;

%% 绘图比较载波后图像 
figure
subplot(3, 1, 1)
plot(0:(1/Fs):(N*T/2-1/Fs), I_carrier);  
legend('载波I signal')   
xlabel('时间t');
ylabel('幅值');
subplot(3, 1, 2)
plot(0:(1/Fs):(N*T/2-1/Fs), Q_carrier);  
legend('载波Q signal')     
xlabel('时间t');
ylabel('幅值');
subplot(3, 1, 3)
plot(0:(1/Fs):(N*T/2-1/Fs), QPSK_signal); 
legend('载波QPSK signal')     
xlabel('时间t');
ylabel('幅值');

%% 通过信道传输
QPSK_recovered = awgn(QPSK_signal, SNR);

%% 相干解调
t2=0:1/Fs:N*T/2-1/Fs;
I_recovered=2*QPSK_recovered.*cos(2*pi*fc*t2);
Q_recovered=-2*QPSK_recovered.*sin(2*pi*fc*t2);

%% 匹配滤波
I_matched = conv(I_recovered, fliplr(h),'same');
Q_matched = conv(Q_recovered, fliplr(h),'same');
QPSK_matched=I_matched+Q_matched;

%% 绘图匹配滤波后图像 
figure();
subplot(3, 1, 1)
plot(0:(1/Fs):(N*T/2-1/Fs), I_matched );  
legend('匹配滤波后I signal') 
xlabel('时间t');
ylabel('幅值');
subplot(3, 1, 2)
plot(0:(1/Fs):(N*T/2-1/Fs), I_matched );  
legend('匹配滤波后Q signal')     
xlabel('时间t');
ylabel('幅值');
subplot(3, 1, 3)
plot(0:(1/Fs):(N*T/2-1/Fs), QPSK_matched);  
legend('匹配滤波后QPSK signal') 
xlabel('时间t');
ylabel('幅值');

%% 抽样判决
% 提取采样点
I_recovered=I_matched(1:sps:end);
Q_recovered=Q_matched(1:sps:end);

%星座图（建议N取10^4以上再看）
QPSK_star=I_recovered+j*Q_recovered;
scatterplot(QPSK_star);
title('QPSK-awgn星座图 SNR=5')
xlim([-4.5,4.5]);
ylim([-4.5,4.5]);

%% 并/串变换
bit_recovered = [];
for i=1:N
    if mod(i,2)~=0
        bit_recovered = [bit_recovered, 2*(I_recovered((i-1)/2+1)>0)-1];  %奇数取I路信息
    else
        bit_recovered = [bit_recovered, 2*(Q_recovered(i/2)>0)-1];        %偶数取Q路信息
    end
end

%% 生成波形绘图
recover_data = [];
for i = 1:N
    recover_data = [recover_data, bit_recovered(i)*ones(1, T*Fs)];
end
I_recover_data = [];
Q_recover_data = [];
for i = 1:N/2
    I_recover_data = [I_recover_data, (2*(I_recovered(i)>0)-1)*ones(1, T*Fs)];
    Q_recover_data = [Q_recover_data, (2*(I_recovered(i)>0)-1)*ones(1, T*Fs)];
end
figure
subplot(3, 1, 1)
plot(t2, I_recover_data);
legend('接收I路信息')   
xlabel('时间t');
ylabel('幅值');
subplot(3, 1, 2)
plot(t2, Q_recover_data); 
legend('接收Q路信息') 
xlabel('时间t');
ylabel('幅值');
subplot(3, 1, 3)
plot(t, recover_data);  
legend('接收比特信息')   
xlabel('时间t');
ylabel('幅值');

%% 计算误码率(误比特率)
BER=bit_error_rate(signal,bit_recovered)

%计算误比特率的函数
function ber = bit_error_rate(sent_bits, received_bits)
% 计算比特之间的差异
errors = sum(sent_bits ~= received_bits);
% 计算误比特率
ber = errors / length(sent_bits);
end

效果：

代码模块（含多普勒频偏）

%% 包含多普勒频偏的仿真
clear;
clc;
close all;
format long

%% 参数设置
N_test = 5000;       % 检测比特数
N_ifm = 10000;          % 消息比特数
N= N_test+N_ifm;
T = 10^-4;     % 比特周期
rolloff = 0.35; % 滚降系数
span = 8; % 滤波器的跨度
sps = 20; % 每个符号的采样点数
fc = 6000;       % 载波频率
fd = 200;  % 多普勒频偏
Rb = 1/T;         % 比特速率
Rs = Rb/2;        % 符号速率
Fs = Rb*sps;      % 抽样频率
t_test = 0:(1/Fs):(N_test*T-1/Fs); % 检测时间
t = 0:(1/Fs):(N*T-1/Fs); % 时间
SNR = 10; % 信噪比/过小影响科斯塔斯环求载频
Ts=1/Fs;

%% 信号发生
signal_test = ones(1, N_test);
signal = [signal_test,randi([0,1], 1, N_ifm)]; % 随机产生的比特数0、1
signal = 2*signal-1;      % 单极性变为双极性-1、1

%% 串/并转换
Isignal=[]; Qsignal=[];
for i=1:N
    if mod(i, 2)~=0
        Isignal = [Isignal, signal(i)];
    else
        Qsignal = [Qsignal, signal(i)];
    end
end

%% 脉冲成型
% 平方根升余弦FIR滤波器
h = rcosdesign(rolloff, span, sps);

% 上采样
Isignal_upsampled = upsample(Isignal , sps);
Qsignal_upsampled = upsample(Qsignal , sps);

% 脉冲整形
Isignal_filtered = conv(Isignal_upsampled, h,'same');
Qsignal_filtered = conv(Qsignal_upsampled, h,'same');

%% 正交调制+fd多普勒频偏
% 载波信号
t2=0:(1/Fs):(N*T/2-1/Fs);
I_carrier = 2*Isignal_filtered.*cos(2*pi*(fc+fd)*t2);       % I路载波信号
Q_carrier = -2*Qsignal_filtered.*sin(2*pi*(fc+fd)*t2);    % Q路载波信号
QPSK_signal = I_carrier + 1j*Q_carrier;

%% 通过信道传输
QPSK_recovered = awgn(QPSK_signal, SNR);

%% 带通滤波器
% 设计巴特沃斯带通滤波器
n = 5; % 滤波器的阶数
Wn = [fc-fd/2, fc+fd/2]/(Fs/2); % 归一化截止频率
b = butter(n, Wn, 'bandpass'); 

% 应用滤波器
I_filtered = filter(b, 1, real(QPSK_recovered));
Q_filtered = filter(b, 1, imag(QPSK_recovered));
QPSK_filtered = I_filtered + 1j*Q_filtered;

%% 求fc＋fd
Sample_Length = length(t_test);
Signal_Channel= 0.5*QPSK_recovered(1:Sample_Length);
symbol_rate = Rb;
peak_freq=COSTAS(length(Signal_Channel),Signal_Channel,symbol_rate,Ts)

%% 相干解调
% 使用fc＋fd相干解调，N_ifm越大越容易放大误差
t2=0:1/Fs:N*T/2-1/Fs;
I_recovered=2*real(QPSK_recovered).*cos(2*pi*peak_freq*t2);
Q_recovered=-2*imag(QPSK_recovered).*sin(2*pi*peak_freq*t2);

%% 匹配滤波
I_matched = conv(I_recovered, fliplr(h), 'same');
Q_matched = conv(Q_recovered, fliplr(h), 'same');
QPSK_matched = I_matched + Q_matched;

%% 抽样判决
I_samples = I_matched(1:sps:end);
Q_samples = Q_matched(1:sps:end);
recovered_bits = [];
for k = 1:N
    if mod(k, 2) == 1
        recovered_bits(k) = 2 * (I_samples((k+1)/2) > 0) - 1;
    else
        recovered_bits(k) = 2 * (Q_samples(k/2) > 0) - 1;
    end
end

%% 计算误码率(误比特率)
BER = bit_error_rate(signal, recovered_bits)

%% 绘制星座图
scatterplot(I_samples + 1j * Q_samples);
title(['QPSK星座图 SNR=', num2str(SNR), 'dB']);
xlabel('In-phase');
ylabel('Quadrature');
grid on;

%% 计算误比特率的函数
function ber = bit_error_rate(sent_bits, received_bits)
    errors = sum(sent_bits ~= received_bits);
    ber = errors / length(sent_bits);
end

%% COSTAS环载波恢复处理过程函数
function peak_freq=COSTAS(Sample_Length,Signal_Channel,symbol_rate,Ts)
decimator = 1;
fs=1/Ts;
% 锁相环参数清零和初始化
Signal_PLL = zeros(1,Sample_Length);
I_PLL = zeros(1,Sample_Length);
Q_PLL = zeros(1,Sample_Length);
NCO_Phase = zeros(1,Sample_Length);
Discriminator_Out = zeros(1,Sample_Length);
Freq_Control = zeros(1,Sample_Length);
PLL_Phase_Part = zeros(1,Sample_Length);
PLL_Freq_Part = zeros(1,Sample_Length);
Ko = 1;
% 压控振荡器增益
Kd = 1;
% 鉴相器增益
K = Ko*Kd;
% 总增益
sigma = 0.707;
% 环路阻尼系数,一般设为0.707
BL = 0.98*symbol_rate;
% 环路等效噪声带宽
Wn = 8*sigma*BL/(1+4*sigma^2);
% 环路自由震荡角频率
T_nco = Ts*decimator;
% 压控振荡器NCO频率字更新周期
K1 = (2*sigma*Wn*T_nco)/(K);
% 环路滤波器系数K1
K2 = ((T_nco*Wn)^2)/(K);
% 环路滤波器系数K2

for i = 2:Sample_Length
    Signal_PLL(i) = Signal_Channel(i)*exp(-j*mod(NCO_Phase(i-1),2*pi));
    % 得到环路滤波器前的相乘器的输入
    % 鉴相器的I路输入信息数据
    I_PLL(i) = real(Signal_PLL(i));
    % 鉴相器的Q路输入信息数据
    Q_PLL(i) = imag(Signal_PLL(i));
    Discriminator_Out(i) = sign(I_PLL(i))*Q_PLL(i)/abs(Signal_PLL(i));
    % 鉴相器的输出误差电压信号
    % 环路滤波器对鉴相器输出的误差电压信号处理后得到锁相环相位响应函数
    PLL_Phase_Part(i) = Discriminator_Out(i)*K1;
    Freq_Control(i) = PLL_Phase_Part(i) + PLL_Freq_Part(i-1);
    % 控制压控振荡器的输出信号频率
    PLL_Freq_Part(i) = Discriminator_Out(i)*K2 + PLL_Freq_Part(i-1);
    % 环路滤波器对鉴相器输出的误差电压信号处理后得到锁相环频率响应函数
    NCO_Phase(i) = NCO_Phase(i-1) + Freq_Control(i);
    % 压控振荡器进行相位调整
end
% COSTAS环载波恢复处理过程到此结束，下面计算提取的载波频率
m=cos(NCO_Phase(end/100:end));
% 计算信号的FFT
N = length(m); % 信号长度
M = fft(m); % 计算FFT
% 计算双边频谱的幅度
P2 = abs(M/N);
% 计算单边频谱的幅度（由于FFT结果是对称的，我们只需要一半的频谱）
P1 = P2(1:N/2+1);
P1(2:end-1) = 2*P1(2:end-1);
% 定义频率向量
f = fs*(0:(N/2))/N;
% 找到频谱的峰值所在频率
[~, idx] = max(P1); % 找到峰值的索引
peak_freq = f(idx); % 对应的频率
end

结果展示：

检测的频率以及误码率如下：