近在做MSK的正交调制解调仿真，有一些心得体会，也踩了一些坑，再次记录一下。
我的目标是，自己写一个正交调制解调，其性能和matlab自带的mskmod和mskdemod函数性能相近，功能相同，输出相同。

一、 mskmod函数初探

发现了一些mskmod函数的一些秘密，在此记录一下

mskmod - Minimum shift keying modulation
This MATLAB function outputs the complex envelope y of the modulation of the
message signal x using differentially encoded minimum shift keying (MSK)
modulation.

y = mskmod(x,nsamp)
y = mskmod(x,nsamp,dataenc)
y = mskmod(x,nsamp,dataenc,ini_phase)
[y,phaseout] = mskmod(...)

从help mskmod函数我们可以得到mskmod函数的使用方法

1. 基础用法

最基础的是直接输入需要调制的信息序列和输入上采样倍数

信息序列：输入需要是01序列，不能是+1/-1这种
上采样：即一个波形采样多少个点，比如一个波形时间为1s，采样率为100Hz，那么上采样倍数就是100，如果一个波形时间为2s，采样率为100Hz，上采样倍数就是50

在实际的仿真代码中，我是这样设置的：

R  = 1e1;   % The bit speed
Fc = 1e1;   % The carrier frequency
Fs = 1e2;   % the sample frequency
Ts = 1/Fs;
Upsample = Fs/R;

2.编码

mskmod函数第三个参数是是否做数据差分编码，默认情况下，mskmod是’diff’方式进行调制，也就是默认情况下，mskmod函数认为你的输入是已经经过差分编码的

’nondiff’ 表示输入的数据是未经差分编码
’diff’ 表示输入的数据是经过差分编码的

但是要注意，这里有一个大坑！！！！
起初我认为，我在mskmod函数里输入diff，是让mskmod函数帮我做一个差分编码，但是仿真结果总是出错，后来我扣了扣mskmod代码注解的字眼

dataenc can be either ‘diff’ for differentially encoded MSK or ‘nondiff’ for nondifferentially encoded MSK.

这里他用的encoded，意思是已经编码过的，非常坑！！
言归正传，由于msk有相位模糊的问题，因此一般是加上差分编码，这里附上差分编码的代码：

% msk差分编码
b0=1;
for i=1:N
    encode_output(i)=b0*bitstream(i);%对应bk
    b0=bitstream(i);
end

这里也有一个误区，有两种常见的差分编码形式，需要好好区分，用对才可以

msk需要的差分编码是如下图这种方式

而另一种差分编码（查分曼彻斯特码）的编码方式通常用在bpsk里：

b0=1;
for i=1:N
    encode_output(i)=-b0*bitstream(i);
    b0=encode_output(i);
end  %差分编码

3.初相

初相这个我没有太多研究，一般都是设置为0，默认也是0

二、 正交调制

经过对mskmod函数的探索，我将mskmod配置到我需要的形式，然后进行正交调制，我的目标是，正交调制的输出和mskmod函数输出相同，参考的是书上这段表格

mskmod函数代码为：

tx_data = mskmod(encode_output1,Upsample,'diff',0*pi);

正交调制的代码为：

    I=[];Q=[];%奇数进I路、偶数进Q路
    bitstream3 = [bitstream1(1),bitstream1];
    I = bitstream3(1:2:length(bitstream3));
    Q = bitstream3(2:2:length(bitstream3));

    %调制
    I_data=[];Q_data=[];
    base_wave=-1/R:1/Fs:1/R-1/Fs;

    % I_data = kron(I,cos(pi*base_wave/(2*T)));
    n = -Upsample:length(I)*2*Upsample-Upsample-1;
    cos_I = cos(pi*n*Ts*R./2);
    I_data_temp = kron(I,ones(1,2*Upsample)).* cos_I;
    I_data = I_data_temp(Upsample+1:end+Upsample*(mod(N,2)-1));%取出零时刻的值
    %此处考虑了总共有奇数个bit需要传输的情况，实际上不需要考虑

    % Q_data = kron(Q,sin(pi*base_wave/(2*R)));
    n = 0:length(Q)*2*Upsample-1;
    sin_Q = sin(pi*n*Ts*R./2);
    Q_data_temp = kron(Q,ones(1,2*Upsample)) .* sin_Q;
    Q_data = Q_data_temp(1:end-Upsample*mod(N,2));%%取出零时刻的值
    %此处考虑了总共有奇数个bit需要传输的情况，实际上不需要考虑

    % Q路延迟半个符号
    number_delay=length(base_wave)/2;
    % Q_data1=[zeros(1,number_delay),Q_data(1:length(Q_data)-number_delay)];

    % 载波信息
    bit_t=0:1/Fs:N/R-1/Fs;
    I_carrier=cos(2*pi*Fc*bit_t);
    Q_carrier=sin(2*pi*Fc*bit_t);

    % 调制
    I_band = I_data.*I_carrier;
    Q_band = Q_data.*Q_carrier;

    % 两路合成一路
    MSK_signal = I_band + 1i*Q_band;

其正交调制框图参考书上的结构：

最终得到的波形图为

从上到下依次是，i路数据，q路数据，两路合并的数据。
通过对比可以发现，与书上是一样的，证明我们做的结果是正确的

3. 解调

解调部分就比较简单，就是正常的正交解调

    I_output = MSK_receive.*I_carrier;
    I_filter_ouput_delay=[I_output(number_delay+1:end) zeros(1,number_delay)];
    Q_output = MSK_receive.*Q_carrier;

这里需要对i路进行一个补偿和延迟，这是因为：
仔细观察书上的那个表格，其实要传输的序列应该是从0时刻开始，而-2和-1时刻是人为补偿的数据，为了能够进行差分编码，因此i路传输的数据其实是6个数据，前面包含了一个自己设置的1，而这个1与我们的数据没有关系，因此在解调是应该去除半个周期的头部，但是这样操作会使得i路数据变短了，无法解调出最后一个数据，因为缺了半个周期，因此我们这里进行补零，补零不影响数据的解调，只是补足了缺少的数据。

解调时，先进行一次本地相关，然后积分判决就可以得到iq数据，代码如下：

%第二次本地相关解调
   I_dected = real(I_filter_ouput_delay).*cos_I(2*number_delay+1:end);
   Q_dected = imag(Q_output).*sin_Q;
   
  %积分判决
    for i=1:round(N/2)
       I_recover(i) = sign(sum(I_dected((i-1)*2*number_delay+1:i*2*number_delay)));
       Q_recover(i) = sign(sum(Q_dected((i-1)*2*number_delay+1:i*2*number_delay)));
    end

% 合并IQ两路
    bit_recover = zeros(1,N);
    bit_recover(1:2:end) = Q_recover;
    bit_recover(2:2:end) = I_recover;
    bit_recover = (bit_recover+1) / 2;

至此，我们完成了所有的代码，可以比对数据得到误码率了
仿真得到的误码率曲线为：

看到，仿真曲线与理论曲线基本一致，因此实现了我们的目标

完整的代码我放在这里：

%%%%%%%%%%%%%%%%%%%%%  MSK调制解调 %%%%%%%%%
%%%%% data:2021年8月10日  author:卡布基诺加勺糖 %%%%%%

%%%%%程序说明
%MSK调制解调器的仿真
%调制方式：MSK 编码：差分编码（正交调制等效于差分编码？？）
%噪声类型：线性高斯白噪声
%解调解调方式：采用正交调制、相干解调

%%%%%仿真环境
% 软件版本：MATLAB2019a

%% %%%%%%%%%%%%%%%%% 参数设定 %%%%%%%%%%%%%%%%%%%%%%
clc;clear all;
tic
R  = 1e1;   % The bit speed
Fc = 1e1; % The carrier frequency
Fs = 1e2;   % the sample frequency
Ts = 1/Fs;
Upsample = Fs/R;

%% %%%%%%%%%%%%%%%%% 主程序 %%%%%%%%%%%%%%%%%%%%%%
EbN0_loop   = [0   2   4   6   8   ]   ;   % signal to noise rate 
N_loop      = [1e4 5e4 1e5 5e5 1e6 ];          %循环点数
for loop=1:length(EbN0_loop)
    %%%%%%%%%%%%%%%%%%% 信源 %%%%%%%%%%%%%%%%%%%%%%
%     loop=1; % use for debug
    % 生成信息序列
    N = N_loop(loop);
    bitstream=randi([0,1],1,N);  % randomly generate the bitstream
    bitstream = [1 1 1 0 0 1 0 1 1 0];
    bitstream1=2*bitstream-1;  % 双极性

    %%%%%%%%%%%%%%%%%%% 编码 %%%%%%%%%%%%%%%%%%%%%%
    % b0=1;
    % for i=1:N
    %     encode_output(i)=-b0*bitstream1(i);
    %     b0=encode_output(i);
    % end  %dpsk差分编码

    b0=1;
    for i=1:N
        encode_output(i)=b0*bitstream1(i);%对应bk
        b0=bitstream1(i);
    end  % msk差分编码

    %%%%%%%%%%%%%%%%%%% 调制 %%%%%%%%%%%%%%%%%%%%%%
    % % 用matlab已有函数mskmod进行调制
    %     encode_output1 = (encode_output + 1)/2;
    %     tx_data = mskmod(encode_output1,Upsample,'diff',0*pi);
    % % mskmod: the elements of x must be 0 or 1

    % 自己写的正交调制
    I=[];Q=[];%奇数进I路、偶数进Q路
    bitstream3 = [bitstream1(1),bitstream1];
    I = bitstream3(1:2:length(bitstream3));
    Q = bitstream3(2:2:length(bitstream3));

    %调制
    I_data=[];Q_data=[];
    base_wave=-1/R:1/Fs:1/R-1/Fs;

    n = -Upsample:length(I)*2*Upsample-Upsample-1;
    cos_I = cos(pi*n*Ts*R./2);
    I_data_temp = kron(I,ones(1,2*Upsample)).* cos_I;
    I_data = I_data_temp(Upsample+1:end+Upsample*(mod(N,2)-1));%取出零时刻的值
    %此处考虑了总共有奇数个bit需要传输的情况，实际上不需要考虑


    n = 0:length(Q)*2*Upsample-1;
    sin_Q = sin(pi*n*Ts*R./2);
    Q_data_temp = kron(Q,ones(1,2*Upsample)) .* sin_Q;
    Q_data = Q_data_temp(1:end-Upsample*mod(N,2));%%取出零时刻的值
    %此处考虑了总共有奇数个bit需要传输的情况，实际上不需要考虑

    % Q路延迟半个符号
    number_delay=length(base_wave)/2;

    % 载波信息
    bit_t=0:1/Fs:N/R-1/Fs;
    I_carrier=cos(2*pi*Fc*bit_t);
    Q_carrier=sin(2*pi*Fc*bit_t);

    % 调制
    I_band = I_data.*I_carrier;
    Q_band = Q_data.*Q_carrier;

    % 两路合成一路
    MSK_signal = I_band + 1i*Q_band;

    %%%%%%%%%%%%%%%%%%%%%%% 信道 %%%%%%%%%%%%%%%%%%%%%%%
    % 加入噪声
    EbN0 = EbN0_loop(loop);
    snr  = EbN0 - 10*log10(Upsample);
    MSK_receive = awgn(MSK_signal,snr,'measured');

    %%%%%%%%%%%%%%%%%%% 解调 %%%%%%%%%%%%%%%%%%%%%%
    I_output = MSK_receive.*I_carrier;
    I_filter_ouput_delay=[I_output(number_delay+1:end) zeros(1,number_delay)];
    I_dected = real(I_filter_ouput_delay).*cos_I(2*number_delay+1:end);
    
    Q_output = MSK_receive.*Q_carrier;
    Q_dected = imag(Q_output).*sin_Q;

    for i=1:round(N/2)
        I_recover(i) = sign(sum(I_dected((i-1)*2*number_delay+1:i*2*number_delay)));
        Q_recover(i) = sign(sum(Q_dected((i-1)*2*number_delay+1:i*2*number_delay)));
    end
    
    bit_recover = zeros(1,N);
    bit_recover(1:2:end) = Q_recover;
    bit_recover(2:2:end) = I_recover;
    bit_recover = (bit_recover+1) / 2;

    error_rate(loop) = sum(abs(bit_recover - bitstream)) / N ;

end

%%%%%%%%%%%%%%%%%画图%%%%%%%%%%%%%%%%%

%save('error_rate.mat','error_rate');
toc
disp('save success')
    
figure()
ber_standard = berawgn(EbN0_loop,'msk','on','coherent');
semilogy(EbN0_loop,ber_standard);
hold on
semilogy(EbN0_loop,error_rate);
legend('参考误码率','仿真误码率')
grid on

写在后面

都看到这里了，我就告诉你一个小秘密，不知道你注意到了没有，代码中没有低通滤波！！

是的，我故意的，没想到吧，但是别急着骂我，上面代码仿真出来的结果肯定和理论值相同的，代码没问题，仿真没问题，结果也没问题，那是我有问题 ？

从理论上来说，与本地载波自相关后肯定要过一个低通滤波来滤除高频项，而且我与本地载波进行了两次自相关，应该要用两个低通滤波才行，这就涉及到下面我说的两个问题了
.
.

这里有我本人也没搞清楚的两个点，希望大家一起讨论：

1. 先说第二次自相关和第二个理论上的低通滤波：
   这里其实我是用累加来代替低通滤波了，高频项在累加过程中，sin或cos函数转了一圈，积分为0，因此变相的完成了低通滤波
2. 第一次与本地载波自相关后我没有加入低通滤波
   其实这里我本来在代码中写了低通滤波器的，并且仿真得到的结果与理论值也是一样的，但是我突然想，如果我去掉这个低通滤波，我只用一个低通滤波可不可以？
   这里以I路的公式为例，如果第一次相关后不做低通滤波，那么高频项会一直保留，在第二次相关时产生高频。一开始我直观的想，我进行累加，第二次相关的2到5项都会积分得到零，所以仿真中我去掉了第一个低通滤波，进行仿真，结果和不进行第一次低通滤波的误码率相同。
   但是实际上，只有第二项可以在累加求和中得到零，其他项由于存在Wc频率的sin或cos，导致积分不为零，因此，从理论上讲这会影响我的误码率，但是仿真结果是正确的，而且我也考虑了Wc和Ts的倍数关系，我将两者设置成不是整倍数的数值，仿真结果也是正确的，这让我很费解
   

有没有大佬解惑，求大佬解惑！！