为什么求出来的高阶累积量和图中的不一样啊，差的很远，psk和fsk都正常的

Question

clear
M = 16;
nSymb = 1000; % 符号数量
data = randi([0 M-1],nSymb,1);
s = qammod(data,M);
%scatterplot(s);
ss1 = awgn(s,10,'measured');
m20 = mean(ss1.*ss1);
m21 = mean(ss1.*conj(ss1));
m40 = mean(ss1.*ss1.*ss1.*ss1);
m41 = mean(ss1.*ss1.*ss1.*conj(ss1));
m42 = mean(ss1.*ss1.*conj(ss1).*conj(ss1));
m60 = mean(ss1.*ss1.*ss1.*ss1.*ss1.*ss1);
m63 = mean(ss1.*ss1.*ss1.*conj(ss1).*conj(ss1).*conj(ss1));

c20 = m20;
c21 = m21;
c40 = m40-3.*m20.^2;
c41 = m41-3.*m21.*m20;
c42 = m42-abs(m20).^2-2.*m21.^2;
c60 = m60-15.*m40.*m20+30.*m20.^3;
c63 = m63-6.*m41.*m20-9.*m42.*m21+18.*m21.*m20.^2+12.*m21.^3;
        
disp([abs(c20),abs(c21),abs(c40),abs(c41),abs(c42),abs(c60),abs(c63)]);

结果是 1.0e+03 *

0.0003    0.0111    0.0570    0.0051    0.0764    0.1358    2.5277

Answer 1

引自免费微信小程序：皆我百晓生

在分析您提供的MATLAB代码和期望的结果之前，让我们先理解一下高阶累积量（Cumulants）的概念。对于一个随机变量X，其高阶累积量定义为：

[ C_n(X) = \frac{1}{\mu^n} E[(X - \mu)^n] ]

其中 ( \mu ) 是随机变量的均值，( E[\cdot] ) 表示数学期望。这里我们计算的是复数信号的高阶累积量，用于描述信号的非高斯特性。

您的代码中，您已经计算了二阶、四阶和六阶累积量。然而，您提到的结果与图表中的数据不匹配，这可能是由于几个原因：

1. 图表中的数据可能对应不同的调制类型，而您的代码针对的是16-QAM调制。
2. 图表中的数据可能是在不同的信噪比（SNR）条件下得到的，而您的代码使用了一个固定的SNR（10dB）。
3. 图表中的累积量可能被归一化，以便于比较不同调制方案的性能。
4. 图表中的累积量可能指的是信号功率谱密度（PSD）上的累积量，而不是实际信号的累积量。

为了进一步分析，我们需要知道图表的具体来源以及它所展示的数据是如何获取的。然而，我们可以尝试调整代码以适应不同的条件并查看结果是否更接近图表中的数据。

首先，我们可以尝试在不同的SNR下运行代码，看看是否能得到类似的结果。例如，您可以创建一个循环，改变AWGN添加器的噪声水平，并记录每个SNR下的累积量。同时，确保您使用正确的调制方式（比如图表中显示的2ASK、BPSK、QPSK等）。然后，您可以绘制累积量随SNR变化的曲线，与图表进行比较。

这是一个修改后的版本，用于计算不同SNR下的累积量：

clear
M = 16;
nSymb = 1000; % 符号数量
SNRs = [0:5:20]; % 不同的SNR范围
c20_all = zeros(1, length(SNRs));
c21_all = zeros(1, length(SNRs));
c40_all = zeros(1, length(SNRs));
c41_all = zeros(1, length(SNRs));
c42_all = zeros(1, length(SNRs));
c60_all = zeros(1, length(SNRs));
c63_all = zeros(1, length(SNRs));

for iSNR = 1:length(SNRs)
    SNR = SNRs(iSNR);
    data = randi([0 M-1], nSymb, 1);
    s = qammod(data, M);
    ss1 = awgn(s, SNR, 'measured');
    
    % 计算累积量
    m20 = mean(ss1.*ss1);
    m21 = mean(ss1.*conj(ss1));
    m40 = mean(ss1.*ss1.*ss1.*ss1);
    m41 = mean(ss1.*ss1.*ss1.*conj(ss1));
    m42 = mean(ss1.*ss1.*conj(ss1).*conj(ss1));
    m60 = mean(ss1.*ss1.*ss1.*ss1.*ss1.*ss1);
    m63 = mean(ss1.*ss1.*ss1.*conj(ss1).*conj(ss1).*conj(ss1));
    
    c20 = m20;
    c21 = m21;
    c40 = m40 - 3*m20.^2;
    c41 = m41 - 3*m21.*m20;
    c42 = m42 - abs(m20).^2 - 2*m21.^2;
    c60 = m60 - 15*m40.*m20 + 30*m20.^3;
    c63 = m63 - 6*m41.*m20 - 9*m42.*m21 + 18*m21.*m20.^2 + 12*m21.^3;

    c20_all(iSNR) = abs(c20);
    c21_all(iSNR) = abs(c21);
    c40_all(iSNR) = abs(c40);
    c41_all(iSNR) = abs(c41);
    c42_all(iSNR) = abs(c42);
    c60_all(iSNR) = abs(c60);
    c63_all(iSNR) = abs(c63);
end

figure;
hold on;
plot(SNRs, c20_all, '-o', 'DisplayName', 'C20');
plot(SNRs, c21_all, '-o', 'DisplayName', 'C21');
plot(SNRs, c40_all, '-o', 'DisplayName', 'C40');
plot(SNRs, c41_all, '-o', 'DisplayName', 'C41');
plot(SNRs, c42_all, '-o', 'DisplayName', 'C42');
plot(SNRs, c60_all, '-o', 'DisplayName', 'C60');
plot(SNRs, c63_all, '-o', 'DisplayName', 'C63');
xlabel('SNR (dB)');
ylabel('Cumulant Magnitude');
legend('show');
grid on;

这个代码段会计算不同SNR下各个累积量的绝对值，并将它们绘制成曲线。您可以根据需要调整SNR范围和步长。

如果您能提供图表的具体来源和解释，我可以为您提供更精确的帮助。此外，为了获得准确的比较，您可能还需要考虑其他因素，如信号的带宽、符号速率等。