今天这个是真正的过了信道的信号。

单载波QPSK。QPSK算是通过了。

简直就是个速成班啊……惨…………

关于CMA的介绍，可以参考Matlab的文档。

Reviews：Matlab中关于CMA的介绍

%%

% 单载波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)

部分代码以后公开。

100MBuad，QPSK。RRC，alpha = 0.3。