通信系统仿真速成第4天：频分复用FDM

别人整天学机器学习、深度学习、数据挖掘……我还（“被”）搞这些…呵呵………

频分复用(FDM，Frequency Division Multiplexing)就是将用于传输信道的总带宽划分成若干个子频带（或称子信道），每一个子信道传输1路信号。频分复用要求总频率宽度大于各个子信道频率之和，同时为了保证各子信道中所传输的信号互不干扰，应在各子信道之间设立隔离带，这样就保证了各路信号互不干扰（条件之一）。频分复用技术的特点是所有子信道传输的信号以并行的方式工作，每一路信号传输时可不考虑传输时延，因而频分复用技术取得了非常广泛的应用。频分复用技术除传统意义上的频分复用(FDM)外，还有一种是正交频分复用(OFDM)。

这里要做的是，传统意义上的FDM。

频分复用的基本思想是：要传送的信号带宽是有限的，而线路可使用的带宽则远远大于要传送的信号带宽，通过对多路信号采用不同频率进行调制的方法，使调制后的各路信号在频率位置上错开，以达到多路信号同时在一个信道内传输的目的。因此，频分复用的各路信号是在时间上重叠而在频谱上不重叠的信号。

以上只是一小点说明。

可参考百度百科：频分复用

%%
% Frequency Division Multiplexing 发送端
% 2017年5月24日18:17:40

clear;
close all;
clc

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

%%

% Set up parameters and signals.
baud_rate = 200; % Baud rate

M1 = 16; % Alphabet size for modulation
M2 = 4; % Alphabet size for modulation
Nsym = 10000;

msg1 = randi([0 M1-1],Nsym,1); % Random message
hMod1 = comm.RectangularQAMModulator(M1);
modmsg1 = step(hMod1,msg1); % Modulate using QAM. % 映射后的基带信号
trainlen = 1000; % Length of training sequence
modmsg1 = modmsg1 ./ mean(abs(modmsg1)); % 所谓的幅度“归一化”

msg2 = randi([0 M2-1],Nsym,1); % Random message
hMod2 = comm.RectangularQAMModulator(M2);
modmsg2 = step(hMod2,msg2); % Modulate using QAM. % 映射后的基带信号
modmsg2 = modmsg2 ./ mean(abs(modmsg2));

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;

f_guard1 = 5;
bw = baud_rate*(1+rolloff);
f_carrier1 = f_guard1 + bw/2;
f_b1 = f_carrier1 - bw/2;
f_e1 = f_carrier1 + bw/2;
f_guard = 10;
f_b2 = f_e1 + f_guard;
f_e2 = f_b2 + bw;
f_carrier2 = (f_b2 + f_e2)/2;

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

% == band1
modmsg_upsample1=upsample(modmsg1,sps);
T_pulsecos1=conv(modmsg_upsample1,rrcFilter);
Tlength=length(T_pulsecos1);
txSig1=T_pulsecos1(CP/2:Tlength-CP/2);

% == band2
modmsg_upsample2=upsample(modmsg2,sps);
T_pulsecos2=conv(modmsg_upsample2,rrcFilter);
Tlength=length(T_pulsecos2);
txSig2=T_pulsecos2(CP/2:Tlength-CP/2);

t = (0:1/fs:((length(txSig1)-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);
cos2 = cos(2*pi*f_carrier2 * t);
sin2 = sin(2*pi*f_carrier2 * t);
x_upconv1 = real(txSig1).* cos1 + imag(txSig1) .* sin1;
x_upconv2 = real(txSig2).* cos2 + imag(txSig2) .* sin2;

x_upconv = x_upconv1 + x_upconv2;

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)');

同前面的文章一样。自己造一个信道出来。假装信号受到了畸变。

整个流程的代码如下。RLS和DFE滤波器均衡。

Review：自适应滤波器RLS

%%
% Frequency Division Multiplexing 接收端
% 2017年5月24日18:31:22

clear;
close all;
clc

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

% Set up parameters and signals.
baud_rate = 200; % Baud rate

M1 = 16; % Alphabet size for modulation
M2 = 4; % Alphabet size for modulation
Nsym = 10000;

msg1 = randi([0 M1-1],Nsym,1); % Random message
hMod1 = comm.RectangularQAMModulator(M1);
modmsg1 = step(hMod1,msg1); % Modulate using QAM. % 映射后的基带信号
modmsg1 = modmsg1 ./ mean(abs(modmsg1)); % 所谓的幅度“归一化”

chan1 = [.986; .845; .237; .123+.31i]; % Channel coefficients
filtmsg1 = filter(chan1,1,modmsg1); % Introduce channel distortion.（已经经过信道的畸变的基带复信号，星座点） 
trainlen = 1000; % Length of training sequence


msg2 = randi([0 M2-1],Nsym,1); % Random message
hMod2 = comm.RectangularQAMModulator(M2);
modmsg2 = step(hMod2,msg2); % Modulate using QAM. % 映射后的基带信号
modmsg2 = modmsg2 ./ mean(abs(modmsg2)); % 所谓的幅度“归一化”
chan2 = [.986; .745; .237; .223+.31i]; % Channel coefficients
filtmsg2 = filter(chan2,1,modmsg2); % Introduce channel distortion.（已经经过信道的畸变的基带复信号，星座点） 


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

f_guard1 = 5;
bw = baud_rate*(1+rolloff);
f_carrier1 = f_guard1 + bw/2;
f_b1 = f_carrier1 - bw/2;
f_e1 = f_carrier1 + bw/2;
f_guard = 10;
f_b2 = f_e1 + f_guard;
f_e2 = f_b2 + bw;
f_carrier2 = (f_b2 + f_e2)/2;
% 2. 脉冲成型
% txSig = upfirdn(modmsg, rrcFilter, sps);  % 发送端的基带复波形信号
% txSig = txSig(CP/2:end-CP/2);

% == band1
% modmsg_upsample1=upsample(modmsg1,sps);
modmsg_upsample1=upsample(filtmsg1,sps);
T_pulsecos1=conv(modmsg_upsample1,rrcFilter);
Tlength=length(T_pulsecos1);
txSig1=T_pulsecos1(rrclen/2:Tlength-rrclen/2);

% == band2
% modmsg_upsample2=upsample(modmsg2,sps);
modmsg_upsample2=upsample(filtmsg2,sps);
T_pulsecos2=conv(modmsg_upsample2,rrcFilter);
Tlength=length(T_pulsecos2);
txSig2=T_pulsecos2(rrclen/2:Tlength-rrclen/2);

t = (0:1/fs:((length(txSig1)-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);
cos2 = cos(2*pi*f_carrier2 * t);
sin2 = sin(2*pi*f_carrier2 * t);
x_upconv1 = real(txSig1).* cos1 + imag(txSig1) .* sin1;
x_upconv2 = real(txSig2).* cos2 + imag(txSig2) .* sin2;

x_upconv = x_upconv1 + x_upconv2;

%%
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)');

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

%%

% 1. 同步
x_received = [x_upconv(300:end);x_upconv;x_upconv;x_upconv;x_upconv;x_upconv];
% x_resampled = resample(x_received,1,1);
x_resampled= x_received;
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,(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);
cos2 = cos(2*pi*f_carrier2 * t);
sin2 = sin(2*pi*f_carrier2 * t);


xi_dnconv1 = x_sync .* cos1;
xq_dnconv1 = x_sync .* sin1;
x_dnconv1 = xi_dnconv1 + 1j * xq_dnconv1;

xi_dnconv2 = x_sync .* cos2;
xq_dnconv2 = x_sync .* sin2;
x_dnconv2 = xi_dnconv2 + 1j * xq_dnconv2;

dn_offset = 0;
rxMatchFilt1=conv(x_dnconv1,rrcFilter);
msg_upsample_len1=length(rxMatchFilt1);
rxFilt1=rxMatchFilt1(rrclen/2:msg_upsample_len1-rrclen/2);
rxFilt1=downsample(rxFilt1,sps,dn_offset);

rxMatchFilt2=conv(x_dnconv2,rrcFilter);
msg_upsample_len2=length(rxMatchFilt2);
rxFilt2=rxMatchFilt2(rrclen/2:msg_upsample_len2-rrclen/2);
rxFilt2=downsample(rxFilt2,sps,dn_offset);

close all;
scatterplot(rxFilt1(trainlen:end));
scatterplot(rxFilt2(trainlen:end));

%% LMS or RLS
close all;

% band1
eq1 = lineareq(40, 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(20, lms(0.001)); % LMS
eq1.SigConst = step(hMod1,(0:M1-1)')'; % Set signal constellation. % 标准星座图
[symbolest1,~] = equalize(eq1,rxLmsInput1,modmsg1(1:trainlen)); % Equalize. % 均衡器obj，需要均衡的信号，训练序列

symbolest1 = symbolest1 ./ mean(abs(symbolest1)) .* mean(abs(eq1.SigConst));
rxFilt_disp1 = rxFilt1 ./ mean(abs(rxFilt1)) .* mean(abs(eq1.SigConst));

symbol_eq_lms1 = symbolest1(trainlen*10+1:end);% 仅仅是个示意，所以取了后面的一些点……前面应该还很烂，后面稍微会好一些。前面的点可能还处于在训练阶段……

% band2
% eq2 = lineareq(100, lms(0.001)); % LMS
eq2 = lineareq(25, rls(0.999,0.01));
eq2.SigConst = step(hMod2,(0:M2-1)')'; % Set signal constellation. % 标准星座图
[symbolest2,~] = equalize(eq2,rxLmsInput2,modmsg2(1:trainlen)); % Equalize. % 均衡器obj，需要均衡的信号，训练序列
symbolest2 = symbolest2 ./ mean(abs(symbolest2)) .* mean(abs(eq2.SigConst));
rxFilt_disp2 = rxFilt2 ./ mean(abs(rxFilt2)) .* mean(abs(eq2.SigConst));
symbol_eq_lms2 = symbolest2(trainlen*10+1:end);

scatterplot(symbol_eq_lms1(10000:end)); % 取后面一些点…………
scatterplot(symbol_eq_lms2(10000:end));

%% DFE
% band1
N1 = 3;
N2 = 5;
dfeObj = dfe(N1,N2,rls(0.999,0.001));
% dfeObj = dfe(N1,N2,lms(0.0001)); 
% Set the signal constellation
dfeObj.SigConst = step(hMod1,(0:M1-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;

rxDfe1 = equalize(dfeObj,symbol_eq_lms1);
rxDfe1 = rxDfe1(2:end); %??

symbolest1 = rxDfe1; % with dfe.

% band2
N1 = 3;
N2 = 5;
dfeObj = dfe(N1,N2,rls(0.999,0.001));
% dfeObj = dfe(N1,N2,lms(0.0001)); 
% Set the signal constellation
dfeObj.SigConst = step(hMod2,(0:M2-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;

rxDfe2 = equalize(dfeObj,symbol_eq_lms2);
rxDfe2 = rxDfe2(2:end); %??

symbolest2 = rxDfe2; % with dfe.

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

h = scatterplot(rxFilt_disp2(trainlen+1:end),1,trainlen,'bx'); hold on;
scatterplot(symbolest2(trainlen+1:end),1,trainlen,'r.',h);
scatterplot(eq2.SigConst,1,0,'y*',h);
legend('Filtered signal','Equalized signal',...
   'Ideal signal constellation');
hold off;

function x_sync = sync_two_signals( x_resampled,x_training_wave,idx,N )
% sync_two_signals( x_resampled,x_training_wave,idx )
% x_resampled：收到的信号
% x_training_wave：用发送的信号
% idx：要找同步上的第几段。idx = 1,2,3,...
% N：要取几段出来

[hicorrI,lagsiI]=xcorr(x_resampled,x_training_wave);
[~,offsetindex]=max(abs(hicorrI));
% clc
close all;
figure;plot(lagsiI,abs(hicorrI));
% offsetindex
for n = 2 : idx
    % figure;plot(lagsiI,abs(hicorrI));
    hicorrI(offsetindex) = 0;
    [~,offsetindex]=max(abs(hicorrI));
end
% offsetindex
h=1;

b = lagsiI(offsetindex)+h;
e = lagsiI(offsetindex)+length(x_training_wave)+h-1;
e = e + length(x_training_wave) * (N-1);
% [b,e]
index = b:e;
% [index(1),index(end)]
x_sync=x_resampled(index);

end