通信系统仿真速成第3天：16QAM调制与解调（实验）

本文链接：https://blog.csdn.net/qcyfred/article/details/72529013

无记忆信源的线性调制，本质上有什么区别吗？

调制就是 z = x + iy;解调就是 x = dot(z,(1,0)); y = dot(z,(0,1))。

正交基，内积……。

发送端。

%%
% 单载波16QAM 发送端
% 2017年5月17日17:55:37

clear;
close all;
clc

rand_seed = 0;
rand('seed',rand_seed);
randn('seed',rand_seed);

%%

% Set up parameters and signals.
baud_rate = 50; % Baud rate
f_carrier1 = 45;

M = 16; % Alphabet size for modulation
Nsym = 10000;

msg = randi([0 M-1],Nsym,1); % Random message
hMod = comm.RectangularQAMModulator(M);
modmsg = step(hMod,msg); % Modulate using QAM. % 映射后的基带信号
trainlen = 1000; % Length of training sequence

rolloff = .3; %滚降因子
span = 20 ; %截断长度
sps = 10;  % Samples per symbol
rrcFilter=rcosdesign(rolloff,span,sps,'sqrt'); %根升余弦滚降滤波器，‘sqrt’均方根升余弦；‘normal’升余弦
CP=length(rrcFilter);
fs = baud_rate*sps; % 时间采样率,时间采样间隔为 1/fs 秒
Tsymbol=1/baud_rate;


% 2. 脉冲成型
% txSig = upfirdn(modmsg, rrcFilter, sps);  % 发送端的基带复波形信号
% txSig = txSig(CP/2:end-CP/2);

T_upsample=upsample(modmsg,sps);
T_pulsecos=conv(T_upsample,rrcFilter);

Tlength=length(T_pulsecos);
txSig=T_pulsecos(CP/2:Tlength-CP/2);

t = (0:1/fs:((length(txSig)-1)/fs)).';
T = t(end)+1/fs;
df = 1/T;
freq = -fs/2:df:fs/2-df;
cos1 = cos(2*pi*f_carrier1 * t);
sin1 = sin(2*pi*f_carrier1 * t);
x_upconv = real(txSig).* cos1 + imag(txSig) .* sin1;

figure(1);
plot(freq,20*log10(abs(fftshift(fft(x_upconv))/max(abs(fftshift(fft(x_upconv)))))));
ylim([-100,10])
xlim([0,freq(end)])
grid on;
title('发送信号');
xlabel('频率(Hz)');

save x_upconv_16qam_50MBd.txt x_upconv -ascii

接收端。

%%
% 单载波16QAM 接收端
% CMA均衡，LMS训练，BER循环计算
% 2017年5月18日10:37:53

clear;
close all;
clc

%%
% x_received = Rx_oscilloscope('osc');% 从示波器读数据
%%

rand_seed = 0;
rand('seed',rand_seed);
randn('seed',rand_seed);

% Set up parameters and signals.

M = 16; % Alphabet size for modulation
baud_rate = 50; % Baud rate
f_carrier1 = 45; % Carrier frequency
Nsym = 10000; % Number of symbols

msg = randi([0 M-1],Nsym,1); % Random message
hMod = comm.RectangularQAMModulator(M);
modmsg = step(hMod,msg); % Modulate using QAM. % 映射后的基带信号
trainlen = 2000; % Length of training sequence

rolloff = .3; % 滚降因子
span = 20 ; % 截断长度
sps = 10;  % Samples per symbol
rrcFilter=rcosdesign(rolloff,span,sps,'sqrt'); %根升余弦滚降滤波器，‘sqrt’均方根升余弦；‘normal’升余弦

fs = baud_rate*sps; % 时间采样率,时间采样间隔为 1/fs 秒
Tsymbol=1/baud_rate;

% 2. 脉冲成型
% txSig = upfirdn(modmsg, rrcFilter, sps);  % 发送端的基带复波形信号
rrcLen=length(rrcFilter);
msg_upsample=upsample(modmsg,sps);
msg_pulse_rrc=conv(msg_upsample,rrcFilter);
msg_upsample_len=length(msg_pulse_rrc);
txSig=msg_pulse_rrc(rrcLen/2:msg_upsample_len-rrcLen/2);

t = (0:1/fs:((length(txSig)-1)/fs)).';
T = t(end)+1/fs;
df = 1/T;
freq = -fs/2:df:fs/2-df;
cos1 = cos(2*pi*f_carrier1 * t);
sin1 = sin(2*pi*f_carrier1 * t);
x_upconv = real(txSig).* cos1 + imag(txSig) .* sin1;

%% === 接收端
x_training_wave = x_upconv;
x_training_msg = modmsg;

%%

% 1. 同步
x_resampled = resample(x_received,1,2);
for k = 1 : 5
    try
        x_sync = sync_two_signals( x_resampled,x_training_wave,k,5);
        is_sync = 1;
    catch
        is_sync = 0;
        fprintf('未同步\n');
    end
    if is_sync == 1
        fprintf('同步\n');
        break;
    end
end

x_sync = x_sync/max(abs(x_sync));
x_training_wave = x_training_wave/max(abs(x_training_wave));
%%
% close all;
% plot(x_sync(1e3:2000),'r');hold on;
% plot(x_training_wave(1e3:2000),'b');

%%

% 2. 下变频 + 匹配滤波
% x_sync = x_training_wave;
%
%%
t = (0:1/fs:((length(x_sync)-1)/fs)).';
T = t(end)+1/fs;
df = 1/T;
freq = -fs/2:df:fs/2-df;

figure(1);
plot(freq,20*log10(abs(fftshift(fft(x_sync))/max(abs(fftshift(fft(x_sync)))))));
ylim([-100,10])
xlim([0,freq(end)])
grid on;
xlabel('频率(Hz)');
title('接收信号');

%%

cos1 = cos(2*pi*f_carrier1 * t);
sin1 = sin(2*pi*f_carrier1 * t);

xi_dnconv = x_sync .* cos1;
xq_dnconv = x_sync .* sin1;
x_dnconv= xi_dnconv + 1j * xq_dnconv;

dn_offset = 0;
rxMatchFilt=conv(x_dnconv,rrcFilter);
msg_upsample_len=length(rxMatchFilt);
rxFilt=rxMatchFilt(rrcLen/2:msg_upsample_len-rrcLen/2);
rxFilt=downsample(rxFilt,sps,dn_offset);

close all;
scatterplot(rxFilt);

%% CMA（qcy）
% CMA_taps = 80;
% nWeights = CMA_taps; 
% mius = [1e-5,1e-6];
% symbolcma = CMA_single_pol_mex(rxFilt,CMA_taps,mius);
% rxCma = symbolcma;
% rxCma = rxCma./mean(abs(rxCma));
% scatterplot(rxCma);

%
%% CMA（Matlab）
close all;
nWeights = 8;
stepSize = 0.0001;
alg = cma(stepSize);
eqCMA = lineareq(nWeights,alg);
eqCMA.SigConst = step(hMod,(0:M-1)')';
eqCMA.leakagefactor = 1;
% eqCMA.ResetBeforeFiltering = 0; % Maintain continuity between iterations.
% eqCMA.Weights = [ones(1,length(eqCMA.Weights)-1) 1];
eqCMA.Weights = [zeros(1,length(eqCMA.Weights)-1),1];
[symbolcma,~] = equalize(eqCMA,rxFilt);

%=========================================================
%                                                  CMA引发相位旋转
% The constant modulus algorithm is useful when no training signal is available, 
% and works best for constant modulus modulations such as PSK.
% However, if CMA has no additional side information, it can introduce phase ambiguity.
% For example, CMA might find weights that produce a perfect QPSK constellation but 
% might introduce a phase rotation of 90, 180, or 270 degrees.
% Alternatively, differential modulation can be used to avoid phase ambiguity.
%=========================================================

rxCma = symbolcma./mean(abs(symbolcma));

scatterplot(rxCma);

%% CMMA -- 原理不清的某某算法。代码暂不公开。
close all;
CMMA_taps = 120;
mius = [1e-6,1e-6];
symbolcmma = CMMA_16QAM_single_pol(rxCma,CMMA_taps,mius);
% scatterplot(symbolcma(nWeights:end))

rxCmma = symbolcmma./mean(abs(symbolcmma));

scatterplot(rxCmma);

%% LMS or RLS
close all;

for k = 1 : 5
    try
        rxLmsInput = sync_two_signals( rxCmma,modmsg,k,4);
        is_sync = 1;
    catch
        is_sync = 0;
        fprintf('未同步\n');
    end
    if is_sync == 1
        fprintf('同步\n');
        break;
    end
end

% rxCma = rxFilt; % 没有CMA
% eq1 = lineareq(20, rls(0.999,0.001));
% eq1 = lineareq(6, rls(0.99,0.001)); % Create an equalizer object. % 40 taps，RLS算法，步长0.99，自相关矩阵逆矩阵的初值InvCorrInit
eq1 = lineareq(9, lms(0.001)); % LMS
eq1.SigConst = step(hMod,(0:M-1)')'; % Set signal constellation. % 标准星座图
[symbolest,~] = equalize(eq1,rxLmsInput,modmsg(1:trainlen)); % Equalize. % 均衡器obj，需要均衡的信号，训练序列

symbolest = symbolest ./ mean(abs(symbolest)) .* mean(abs(eq1.SigConst));
rxFilt_disp = rxFilt ./ mean(abs(rxFilt)) .* mean(abs(eq1.SigConst));

symbol_eq_lms = symbolest(trainlen*2+1:end);

scatterplot(symbol_eq_lms);


%%
% DFE
N1 = 3;
N2 = 5;
dfeObj = dfe(N1,N2,rls(0.999,0.0001));
% dfeObj = dfe(N1,N2,lms(0.0001)); 
% Set the signal constellation
dfeObj.SigConst = step(hMod,(0:M-1)')';
% Maintain continuity between calls to equalize
dfeObj.ResetBeforeFiltering = 0;
% Define initial coefficients to help convergence
dfeObj.Weights = zeros(1,N1+N2);
dfeObj.Weights(2) = 1;

rxDfe = equalize(dfeObj,symbol_eq_lms);
rxDfe = rxDfe(2:end); %??

symbolest = rxDfe; % with dfe.


%%

close all;
h = scatterplot(rxFilt_disp(trainlen+1:end),1,trainlen,'bx'); hold on;
scatterplot(symbolest(trainlen+1:end),1,trainlen,'r.',h);
scatterplot(eq1.SigConst,1,0,'y*',h);
legend('Filtered signal','Equalized signal',...
   'Ideal signal constellation');
hold off;

%%
% 计算BER

一个没有数学表达式来源的算法，可以用？

CMA针对恒模调制。针对16QAM这种三个圈的调制，可能会不稳定。