通信系统仿真速成第2天：QPSK调制与解调（实验）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
%%
% 单载波QPSK 接收端
% CMA均衡，LMS训练，BER循环计算
% 2017年5月18日10:37:53

clear;
close all;
clc

%%

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

% Set up parameters and signals.

M = 4; % Alphabet size for modulation
baud_rate = 100; % Baud rate
f_carrier1 = 75; % 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 = 3000; % 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;

% x_received = Rx_oscilloscope('osc');
load ./rxdata/oscCAP_7osc.txt;
x_received=oscCAP_7osc;
%%

% 1. 同步
x_resampled = resample(x_received,1,1);
for k = 1 : 5
    try
        x_sync = sync_two_signals( x_resampled,x_training_wave,k,3);
        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（Matlab）
close all;
nWeights = 20;
stepSize = 0.001;
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(nWeights:end); % 
rxCma = rxCma./mean(abs(rxCma));
scatterplot(rxCma(nWeights+1:end));

%% LMS or RLS
close all;

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

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

%%
% Plot signals.
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,'k*',h);
legend('Filtered signal','Equalized signal',...
   'Ideal signal constellation');
hold off;

%%
% a1 = symbolest(trainlen+1:end);
% a2 = modmsg;
% 之后如何计算BER就不用写了。方法非常magic。

%% 没有相位旋转的时候计算BER
% % Compute error rates with equalization.
% hDemod = comm.RectangularQAMDemodulator(M);
% demodmsg = step(hDemod,symbolest); % Demodulate detected signal from equalizer.
% 
% % demodmsg
% % msg
% 
% % [hicorrI,lagsiI]=xcorr(demodmsg,msg);
% % [~,offsetindex]=max((hicorrI));
% % figure;plot(lagsiI,abs(hicorrI));
% 
% % Create ErrorRate Calculator System object
% serVec = step(comm.ErrorRate,msg(trainlen+1:end),demodmsg(trainlen+1:end));
% srate = serVec(1)
% snum = serVec(2)
% % Convert integers to bits
% hIntToBit = comm.IntegerToBit(log2(M));
% Tx_bit = step(hIntToBit, msg(trainlen+1:end));
% Rx_bit = step(hIntToBit, demodmsg(trainlen+1:end));
% % Calculate BER
% berVec = step(comm.ErrorRate,Rx_bit,Tx_bit);
% brate = berVec(1)
% bnum = berVec(2)