调制识别HOM去噪方法

本文章复现的是Sofiane Kharbech等，于2020年在IEEE Wireless Communications Letters发表的论文《Denoising Higher-Order Moments for Blind Digital Modulation Identification in Multiple-Antenna Systems》（多天线系统中数字调制盲识别的高阶矩去噪方法）

论文复现

本文章复现的是Sofiane Kharbech等，于2020年在IEEE Wireless Communications Letters发表的论文《Denoising Higher-Order Moments for Blind Digital Modulation Identification in Multiple-Antenna Systems》（多天线系统中数字调制盲识别的高阶矩去噪方法）

首先编码三个函数文件，分别实现：

  1. f_HOS_Extraction.m
    从接收信号中提取高阶统计量（HOS）

    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
    function HOS=f_HOS_Extraction(M,SNR,K,Nt,Nr,pbee,cfo,phznoise)
    switch M
       case 1
          x = randi([0,1],1,Nt*K);
          x_mod = pskmod(x,2);
       case 2
          x = randi([0,3],1,Nt*K);
          x_mod = pskmod(x,4);
       case 3
          x = randi([0,7],1,Nt*K);
          x_mod = pskmod(x,8);
       case 4
          x = randi([0,3],1,Nt*K);
          x_mod = pammod(x,4);
       case 5
          x = randi([0,7],1,Nt*K);
          x_mod = pammod(x,8);
       case 6
          x = randi([0,15],1,Nt*K);
          x_mod = qammod(x,16);
       otherwise
          disp('Unknown modulation type');
    end
    scale = modnorm(x_mod,'avpow',1);   % scaling factor
    x_mod = x_mod*scale;                % scaling avg. power to 1W
    %------------------------------------------------> Channel (MIMO + AWGN)
    H = randn(Nr,Nt)+1i*randn(Nr,Nt);
    y_mimo = H*reshape(x_mod,Nt,[]);
    cfo_mat = repmat(exp(2i*pi*(1:1:K)*cfo),Nr,1);
    y_mimo = y_mimo.*cfo_mat;
    if phznoise ~=0
       for cptNr=1:Nr
           y_mimo(cptNr,:) = phznoise(y_mimo(cptNr,:).').';
       end
    end
    y_mimo_SNR = awgn(y_mimo,SNR,'measured');
    Pb = mean(abs(reshape(y_mimo_SNR-y_mimo,1,[]).^2)) + pbee; % noise power

    R = mean(real(x_mod).^4) / mean(real(x_mod).^2);
    [W,yeg] = f_SCMA(Nt,Nr,y_mimo_SNR,R);
    Pbf = Pb*W.'*conj(W); % power of the filtered noise

    % bf = W.'*(y_mimo_SNR-y_mimo)  % fchk
    % mean(abs(bf(1,:).^2))         % fchk

    %----> HOS extraction

    for cpt=1:Nt
        HOS(cpt,:) = f_CalcHOS(yeg(cpt,:),Pbf(cpt,cpt));
    end

  2. f_CalcHOS.m
    该函数计算输入信号的各种高阶矩和累积量。它计算矩，如m20， m40, m60，以及累积量，如c40， c60等。

    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
    % x  : the signal
    % Pb : the noise power
    function hos=f_CalcHOS(x, Pb)
    x=x-mean(x);

    m20=real(mean(x.^2));
    m22=real(mean(conj(x).^2));
    m40=real(mean(x.^4));
    m60=real(mean(x.^6));

    m21y=real(mean(x.*(conj(x))));
    m41y=real(mean((x.^3).*(conj(x))));
    m42y=real(mean((x.^2).*(conj(x).^2)));
    m43y=real(mean(x.*(conj(x).^3)));
    m61y=real(mean((x.^5).*(conj(x))));
    m62y=real(mean((x.^4).*(conj(x).^2)));
    m63y=real(mean((x.^3).*(conj(x).^3)));
    m84y=real(mean((x.^4).*(conj(x).^4)));

    % denoising
    m21=m21y-Pb;
    m41=m41y-3*m20*Pb;
    m42=m42y-2*Pb^2-4*Pb*m21;
    m43=m43y-3*m22*Pb;
    m61=m61y-5*m40*Pb;
    m62=m62y-12*m20*Pb^2-8*m41*Pb;
    m63=m63y-18*m21*Pb^2-9*m42*Pb-6*Pb^3;
    m84=m84y-16*m63*Pb-72*m42*Pb^2-96*m21*Pb^3-24*Pb^4;

    %--------------------------> Cumulants <--------------------------%

    c40=m40-3*(m20.^2);
    c60=m60-15*m20*m40+30*(m20.^3);

    c41=m41-3*m20*m21;
    c42=m42-(abs(m20).^2)-2*(m21.^2);
    c61=m61-5*m21*m40-10*m20*m41+30*m21*(m20.^2);
    c62=m62-6*m42*m20-8*m21*m41-m22*m40+6*m22*m20^2+24*m20*m21^2;
    c63=m63-9*m42*m21+12*(m21.^3)-3*m20*m43-3*m22*m41+18*m20*m21*m22;

    %--------------------------> Normalization <--------------------------%

    m40_n=m40/m21^2;
    m41_n=m41/m21^2;
    m42_n=m42/m21^2;
    m60_n=m60/m21^3;
    m61_n=m61/m21^3;
    m62_n=m62/m21^3;
    m63_n=m63/m21^3;
    m84_n=m84/m21^4;

    c40_n=c40/m21^2;
    c41_n=c41/m21^2;
    c42_n=c42/m21^2;
    c60_n=c60/m21^3;
    c61_n=c61/m21^3;
    c62_n=c62/m21^3;
    c63_n=c63/m21^3;

    m41y_n=m41y/m21y^2;
    m42y_n=m42y/m21y^2;
    m61y_n=m61y/m21y^3;
    m62y_n=m62y/m21y^3;
    m63y_n=m63y/m21y^3;
    m84y_n=m84y/m21y^4;

    hos=real([m40_n m41_n m41y_n m42_n m42y_n m60_n m61_n m61y_n m62_n m62y_n m63_n m63y_n m84_n m84y_n ...
              c40_n c41_n c42_n c60_n c61_n c62_n c63_n]);

  3. f_SCMA.m
    该函数为MIMO系统执行信号盲源分离（BSS）。利用提取的信号特征和噪声特征，对接收信号进行白化变换，然后利用随机梯度下降（SGD）估计BSS的分离矩阵W。

    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
    function [G,yeg]=f_SCMA(Nt,Nr,y_mimo_SNR,R)
    %----------------------------------> Pre-processing: batch whitening
    Ry = (y_mimo_SNR*y_mimo_SNR')/size(y_mimo_SNR,2);
    [U,SIG] = eig(Ry);                          % Eigendecomposition
    SIG = SIG(Nr:-1:(Nr-Nt+1),Nr:-1:(Nr-Nt+1)); % Required condition: Nr>=Nt
    U = U(:,Nr:-1:(Nr-Nt+1));
    B = (SIG^-0.5)*U'; % Whitening matrix
    YW = B*y_mimo_SNR; % Whited signals

    %----------------------------------> Estimating the separator W    
    W = eye(Nt);                       % Initialization
    % W = toeplitz(1:Nt);
    % W = toeplitz(1:Nr,1:Nt);
    mu=1e-2;                           % Step size
    % mu=5e-3;
    for ns=1:size(YW,2)                % Stochastic gradient descent (SGD)
        YW_ns = YW(:,ns);
        z_ns = W.'*YW_ns;
        for n=1:Nt
            e = (real(z_ns(n)).^2-R).*real(z_ns(n));
            W(:,n) = W(:,n)-mu.*e.*conj(YW_ns);
        end
        W = f_GramSchmidt(W); % Post-processing: orthonormalization
    end

    %----------------------------------> Applying BSS
    yeg = W.'*YW;
    G = B.'*W;
    function Q = f_GramSchmidt(W)
    Q = zeros(size(W));
    Q(:,1) = W(:,1)/norm(W(:,1));
    for i=2:size(W,2)
        s_proj = 0;
        for j=1:i-1
            s_proj = s_proj+((Q(:,j)'*W(:,i))/(Q(:,j)'*Q(:,j)))*Q(:,j);
        end
        e = W(:,i)-s_proj;
        Q(:,i) = e/norm(e);
    end

随后，我们需要生成高阶统计（HOS）特征的数据集。模拟各种调制方案下的MIMO传输，计算每个传输信号的HOS特征，并将数据组织成矩阵。生成的数据集是一个三维矩阵（mat_HOS），形状为==特征向量 × 样本指标 × MIMO系统的天线指标==。

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
% 1) 本文件用于生成指定数字调制方式和传输配置的高阶统计量(HOS)数据集
% 2) 生成的数据集名为'mat_HOS'，为三维矩阵：
%    - 第一维：HOS向量（特征维度）
%    - 第二维：样本索引（HOS向量编号）
%    - 第三维：天线编号（MIMO系统接收天线）
% 3) 若'rng shuffle'不可用，使用以下替代方案：
% clk=clock;
% RandStream.setGlobalStream(RandStream('mt19937ar','Seed',str2num(strrep(num2str(clk(4:6)),' ',''))*1e3));
% 4) 若未自动启用并行处理，执行以下操作：
% parpool('local',3)        % 启动3个本地工作进程
% Linux系统可能需要先运行：
% ! sudo chmod -R 777 /home/user_name/.matlab/local_cluster_jobs/R2017a
% 注：若无并行计算工具箱，需将最内层循环的parfor改为for
% 5) 运行结束后保存工作区：
% 需保留所有变量，供DetectionResults.m生成结果
rng shuffle
clc
clear
%--------------------------------------------> Parameters
Nt = 2;    % 发射天线数量
Nr = 6;    % 接收天线数量
K = 4000;  % 传输的MIMO符号数量
M = [1 2 3 4 5 6]; % 调制类型索引（对应不同调制方式）, cf. 'f_HOS_Extraction.m' for the names
% M = [1 2 3];
SNR = -2:1:15;     % 信噪比范围（-2dB到15dB）
lMC = 2000;        % 蒙特卡洛模拟次数
eestd = 0;         % 噪声功率估计误差标准差
cfo = 0;           % 载波频率偏移
% cfo = 1e-4;
phznoise = 0;      % 相位噪声模型
% phznoise = comm.PhaseNoise('Level',-3,'FrequencyOffset',2e-3,'SampleRate',1);

%--------------------------------------------> Main processing
mat_HOS = []; % HOS特征（特征×样本×天线）
for cptSNR=1:length(SNR)
    for cptM=1:length(M)
        Pb_ee = eestd*randn(1,lMC);
        parfor cptMC=1:lMC
            % 特征提取函数
            HOS = f_HOS_Extraction(M(cptM),SNR(cptSNR),K,Nt,Nr,Pb_ee(cptMC),cfo,phznoise);
            % 拼接特征矩阵（第三维为天线，第二维为样本）
            mat_HOS = [mat_HOS; permute(HOS,[3 2 1])];    
        end
        disp(['M = ' num2str(M(cptM)) ' completed.'])       
    end
end
mat_HOS = permute(mat_HOS,[2 1 3]);% 调整维度为[特征数×样本数×天线数]
clearvars -except M SNR Nt Nr lMC mat_HOS

最后，我们使用HOS数据集（mat_HOS）根据计算的特征对调制类型进行分类。它通过绘制正确识别概率（PCI）作为信噪比的函数来评估调制识别的性能。

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
% 1) 本文件用于绘制特定场景下正确识别概率（PCI）随信噪比（SNR）变化曲线
% 2) 支持三种特征组合场景：
%    - '仅使用高阶累积量（HOC）'
%    - 'HOC + 未去噪高阶矩（HOM）'
%    - 'HOC + 去噪后HOM'
% 3) 运行前需加载'GenDataset.m'生成的所有变量
HOSth = { % Theoretical HOS值，对应论文表2
'   '  'B-PSK' 'Q-PSK' '8-PSK' '4-ASK' '8-ASK' '16-QAM'    ;
'M40'    1      -1       0       1.64    1.77    -0.67     ;
'M41'    1       0       0       1.64    1.77     0        ;
'M42'    1       1       1       1.64    1.77     1.32     ;
'M60'    1       0       0       2.92    3.62     0        ;
'M61'    1      -1       0       2.92    3.62    -1.32     ;
'M62'    1       0       0       2.92    3.62     0        ;
'M63'    1       1       1       2.92    3.62     1.96     ;
'M84'    1       1       1       5.25    7.92     3.12     ;
'C40'   -2      -1       0      -1.36   -1.42    -0.68     ;
'C41'   -2       0       0      -1.36   -1.42     0        ;
'C42'   -2      -1      -1      -1.36   -1.42    -0.68     ;
'C60'    16      0       0       8.32    7.19     0        ;
'C61'    16      4       0       8.32    7.19     2.08     ;
'C62'    16      0       0       8.32    7.19     0        ;
'C63'    16      4       4       8.32    7.19     2.08     ;
};
HOSth = cell2mat(HOSth(2:end,2:end));
lSNR = length(SNR); % SNR点数
lM = length(M); % 调制类型数
lMMC = lM*lMC; % 总样本数 = 调制类型数*蒙特卡洛次数
stat = zeros(4+Nt,lSNR*lMMC); % 统计矩阵初始化
count = 1;
for cptSNR=1:lSNR
    for cptM=1:lM
        for cptMC=1:lMC           
            stat(1,count) = SNR(cptSNR);% 记录SNR
            stat(2,count) = M(cptM); % 记录真实调制索引
            count = count+1;      
        end
    end
end
%% 下面的代码使用那种方式就取消注释哪种方式即可
% HOSth = HOSth(9:end,:);        % 'HOC only'时取消注释
for cptNt=1:Nt
    % HOC only
    %mat_HOS_cptNt = mat_HOS(15:end,:,cptNt);
    % HOC + 未去噪HOM
    mat_HOS_cptNt = mat_HOS([1 3 5 6 8 10 12 14 15 16 17 18 19 20 21],:,cptNt);
    % HOC + 去噪HOM
    % mat_HOS_cptNt = mat_HOS([1 2 4 6 7 9 11 13 15 16 17 18 19 20 21],:,cptNt);
    [~,resUz] = min(pdist2(HOSth',mat_HOS_cptNt','euclidean')); % MD classifier
    stat(2+cptNt,:)=resUz;
end
% 多天线结果融合
if Nt==1
   stat(4,:) = stat(3,:); % 单天线直接复制
else
   stat(3+Nt,:) = mode(stat(3:2+Nt,:)); % frequency (number of occurrences)
end
stat(4+Nt,:) = double(stat(2,:)==stat(3+Nt,:)); % 生成二进制正确性标记向量（1=正确，0=错误）
pci = zeros(1,lSNR);
end_ = 0;
% 按SNR分段统计正确识别率
for cptSNR=1:lSNR
    beg_ = end_+1;
    end_ = beg_+lMMC-1; % 当前SNR对应的数据段
    pci(cptSNR) = sum(stat(4+Nt,beg_:end_))/lMMC; % 正确率
end
figure
plot(SNR,pci,'-o','LineWidth',1.5)   % 添加样式：实线+圆圈标记
xlabel('SNR', 'FontSize', 12)
ylabel('正确识别概率', 'FontSize', 12)
grid on
set(gca, 'FontWeight', 'bold')       % 坐标轴字体加粗

到这里，对论文的复现就完成了。

One More Thing

如果需要实现对某一信号分类的效果，也就是输入信号x，输出信号的调制方式1-6的效果，可以对代码做如下改进。

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
% 主函数，用于调制方式识别
function modulation_type = modulate_and_identify(input_signal, SNR, Nt, Nr, K, pbee, cfo, phznoise)
    % 1. 生成高阶统计量（HOS）
    HOS = f_HOS_Extraction(input_signal, SNR, K, Nt, Nr, pbee, cfo, phznoise);
    % 2. 理论HOS值，用于后续匹配
    HOSth = { 
        '   '  'B-PSK' 'Q-PSK' '8-PSK' '4-ASK' '8-ASK' '16-QAM';
        'M40'   1     -1       0       1.64   1.77   -0.67;
        'M41'   1      0       0       1.64   1.77    0;
        'M42'   1      1       1       1.64   1.77    1.32;
        'M60'   1      0       0       2.92   3.62    0;
        'M61'   1     -1       0       2.92   3.62   -1.32;
        'M62'   1      0       0       2.92   3.62    0;
        'M63'   1      1       1       2.92   3.62    1.96;
        'M84'   1      1       1       5.25   7.92    3.12;
        'C40'  -2     -1       0      -1.36  -1.42   -0.68;
        'C41'  -2      0       0      -1.36  -1.42    0;
        'C42'  -2     -1      -1      -1.36  -1.42   -0.68;
        'C60'  16      0       0       8.32   7.19    0;
        'C61'  16      4       0       8.32   7.19    2.08;
        'C62'  16      0       0       8.32   7.19    0;
        'C63'  16      4       4       8.32   7.19    2.08;
    };
    HOSth = cell2mat(HOSth(2:end,2:end));
    % 3. 使用高阶矩进行调制识别
    [~, resUz] = min(pdist2(HOSth', HOS', 'euclidean')); % 最小距离分类法
    % 4. 输出识别结果
    modulation_type = get_modulation_type(resUz);
    disp(['识别的调制方式为: ', modulation_type]);
end

% 调制方式的映射函数
function modulation_type = get_modulation_type(index)
    modulation_types = {'B-PSK', 'Q-PSK', '8-PSK', '4-ASK', '8-ASK', '16-QAM'};
    modulation_type = modulation_types{index}; % 获取对应的调制方式
end

% HOS提取函数
function HOS = f_HOS_Extraction(input_signal, SNR, K, Nt, Nr, pbee, cfo, phznoise)
    % 信号预处理
    x = input_signal;  % 输入信号
    scale = modnorm(x, 'avpow', 1);   % 归一化
    x_mod = x * scale;  % 归一化信号
    % MIMO信道模拟
    H = randn(Nr, Nt) + 1i * randn(Nr, Nt);
    y_mimo = H * reshape(x_mod, Nt, []);
    % 载波频率偏移（CFO）和相位噪声模型
    cfo_mat = repmat(exp(2i * pi * (1:1:K) * cfo), Nr, 1);
    y_mimo = y_mimo .* cfo_mat;
    if phznoise ~= 0
        for cptNr = 1:Nr
            y_mimo(cptNr, :) = phznoise(y_mimo(cptNr, :).');
        end
    end
    % 加噪声
    y_mimo_SNR = awgn(y_mimo, SNR, 'measured');
    Pb = mean(abs(reshape(y_mimo_SNR - y_mimo, 1, []).^2)) + pbee; % 噪声功率估计
    % BSS（盲源分离）
    R = mean(real(x_mod).^4) / mean(real(x_mod).^2);
    [W, yeg] = f_SCMA(Nt, Nr, y_mimo_SNR, R);
    Pbf = Pb * W.' * conj(W); % 滤波后的噪声功率
    % HOS特征提取
    for cpt = 1:Nt
        HOS(cpt, :) = f_CalcHOS(yeg(cpt, :), Pbf(cpt, cpt));
    end
end

% SCMA函数：用于盲源分离
function [G, yeg] = f_SCMA(Nt, Nr, y_mimo_SNR, R)
    Ry = (y_mimo_SNR * y_mimo_SNR') / size(y_mimo_SNR, 2);
    [U, SIG] = eig(Ry);  % 特征分解
    SIG = SIG(Nr:-1:(Nr-Nt+1), Nr:-1:(Nr-Nt+1)); % 确保Nr>=Nt
    U = U(:, Nr:-1:(Nr-Nt+1));
    B = (SIG^-0.5) * U';  % 白化矩阵
    YW = B * y_mimo_SNR;  % 白化信号
    % 初始化分离矩阵
    W = eye(Nt);
    mu = 1e-2;  % 步长
    % 随机梯度下降（SGD）算法
    for ns = 1:size(YW, 2)
        YW_ns = YW(:, ns);
        z_ns = W.' * YW_ns;
        for n = 1:Nt
            e = (real(z_ns(n)).^2 - R) .* real(z_ns(n));
            W(:, n) = W(:, n) - mu * e .* conj(YW_ns);
        end
        W = f_GramSchmidt(W); % 正交化处理
    end
    yeg = W.' * YW;
    G = B.' * W;
end

% Gram-Schmidt正交化处理
function Q = f_GramSchmidt(W)
    Q = zeros(size(W));
    Q(:, 1) = W(:, 1) / norm(W(:, 1));

    for i = 2:size(W, 2)
        s_proj = 0;
        for j = 1:i-1
            s_proj = s_proj + ((Q(:, j)' * W(:, i)) / (Q(:, j)' * Q(:, j))) * Q(:, j);
        end

        e = W(:, i) - s_proj;
        Q(:, i) = e / norm(e);
    end
end

% HOS计算函数
function hos = f_CalcHOS(x, Pb)
    x = x - mean(x);

    m20 = real(mean(x.^2));
    m22 = real(mean(conj(x).^2));
    m40 = real(mean(x.^4));
    m60 = real(mean(x.^6));

    m21y = real(mean(x .* (conj(x))));
    m41y = real(mean((x.^3) .* (conj(x))));
    m42y = real(mean((x.^2) .* (conj(x).^2)));
    m43y = real(mean(x .* (conj(x).^3)));
    m61y = real(mean((x.^5) .* (conj(x))));
    m62y = real(mean((x.^4) .* (conj(x).^2)));
    m63y = real(mean((x.^3) .* (conj(x).^3)));
    m84y = real(mean((x.^4) .* (conj(x).^4)));

    m21 = m21y - Pb;
    m41 = m41y - 3 * m20 * Pb;
    m42 = m42y - 2 * Pb^2 - 4 * Pb * m21;
    m43 = m43y - 3 * m22 * Pb;
    m61 = m61y - 5 * m40 * Pb;
    m62 = m62y - 12 * m20 * Pb^2 - 8 * m41 * Pb;
    m63 = m63y - 18 * m21 * Pb^2 - 9 * m42 * Pb - 6 * Pb^3;
    m84 = m84y - 16 * m63 * Pb - 72 * m42 * Pb^2 - 96 * m21 * Pb^3 - 24 * Pb^4;

    c40 = m40 - 3 * (m20.^2);
    c60 = m60 - 15 * m20 * m40 + 30 * (m20.^3);

    c41 = m41 - 3 * m20 * m21;
    c42 = m42 - (abs(m20).^2) - 2 * (m21.^2);
    c61 = m61 - 5 * m21 * m40 - 10 * m20 * m41 + 30 * m21 * (m20.^2);
    c62 = m62 - 6 * m42 * m20 - 8 * m21 * m41 - m22 * m40 + 6 * m22 * m20^2 + 24 * m20 * m21^2;
    c63 = m63 - 9 * m42 * m21 + 12 * (m21.^3) - 3 * m20 * m43 - 3 * m22 * m41 + 18 * m20 * m21 * m22;

    hos = real([m40 m41 m42 m60 m61 m62 m63 m84 c40 c41 c42 c60 c61 c62 c63]);
end

Nt = 2;       % 发射天线数量
Nr = 6;       % 接收天线数量
K = 4000;     % 传输的MIMO符号数量
SNR = 10;     % 信噪比（dB）
pbee = 0.01;  % 噪声功率估计误差
cfo = 1e-4;   % 载波频率偏移
phznoise = 0; % 相位噪声模型（0代表无相位噪声）
% 生成一个模拟的输入信号（例如，16-QAM调制）
M = 6; % 假设已知的调制类型为16-QAM（M = 6）
x = randi([0, 15], 1, Nt * K); % 生成16-QAM信号的符号
x_mod = qammod(x, 16);          % 调制为16-QAM信号
input_signal = x_mod;           % 输入信号
% 调用调制识别函数
modulation_type = modulate_and_identify(input_signal, SNR, Nt, Nr, K, pbee, cfo, phznoise);
% 显示结果
disp(['识别的调制方式为: ', modulation_type]);
论文复现 > 信号处理
#论文复现 #信号处理 #调制识别
调制识别HOM去噪方法
http://example.com/2025/05/24/调制识别HOM去噪方法/
作者
Yifan Xie
发布于
2025年5月24日
许可协议