[实战] M-QAM信号生成（完整C代码）

M-QAM信号生成

不依赖任何第三方库，用C语言生成各种QAM信号，并使用python进行数据频谱分析，确认C语言实现的正确性。
话不多说，实战开始。

M-QAM调制原理与实现

M-QAM调制原理

M-QAM（正交幅度调制）是一种将数字信号映射到正交载波上的调制技术，通过同时改变载波的幅度和相位来传输信息。

核心原理：

1. 星座图映射：将输入的比特流分组，每组k=log₂(M)比特映射为复平面上的一个点
   $\text{scale}=\frac{A}{\sqrt{\frac{2(M-1)}{3}}}$

2. 调制公式：
   $s(t)=I(t)\cdot \cos (2\pi f_{c}t)-Q(t)\cdot \sin (2\pi f_{c}t)$

3. 调制特性：

   • 频谱效率：$\eta ={\log }_2(M)\text{bps/Hz}$
   • 误码率与信噪比关系：$P_b\approx \frac{4}{{\log }_{2}M}\left(1-\frac{1}{\sqrt{M}}\right)Q\left(\sqrt{\frac{3\cdot \text{SNR}\cdot {\log }_{2}M}{M-1}}\right)$
     其中：

• $Q(x)$ 是高斯误差函数，表征噪声影响
• SNR为信噪比，调制阶数 $M$ 增大时，相同SNR下误码率显著上升
• 公式反映了高阶调制在提升频谱效率的同时需牺牲可靠性

5. 升余弦滤波器：
   $h(t)=\{\begin{array}{ll}1-\alpha +\frac{4\alpha }{\pi }& t=0\\ \frac{\alpha }{\sqrt{2}}\left[(1+\frac{2}{\pi })\sin (\frac{\pi }{4\alpha })+(1-\frac{2}{\pi })\cos (\frac{\pi }{4\alpha })\right]& |t|=\frac{T}{4\alpha }\\ \frac{\sin (\pi t/T(1-\alpha ))+4\alpha t/T\cos (\pi t/T(1+\alpha ))}{\pi t/T[1-(4\alpha t/T{)}^2]}& \text{otherwise}\end{array}$
6. AWGN生成：
   $\text{noise}=\sigma \times \sqrt{-2\ln U_1}\cos (2\pi U_2)$
   其中 $\sigma =\sqrt{\frac{P_{\text{signal}}}{10^{{\text{SNR}}_{\text{dB}}/10}}}$

常见可配置参数

C语言实现 (mqam_generator.c)

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <complex.h>
#include <string.h>

typedef struct {
    int M;               // 调制阶数
    int num_symbols;     // 符号数量
    int samples_per_symbol; 
    double carrier_freq; // 载波频率 (Hz)
    double symbol_rate;  // 符号速率 (Baud)
    double snr_db;       // 信噪比 (dB)
    double carrier_phase;// 载波相位 (度)
    double amplitude;    // 信号幅度
    int pulse_shaping;   // 脉冲成型
    double rolloff;      // 升余弦滚降系数
} MQAMParams;

// 升余弦脉冲成型
double rrc_pulse(double t, double T, double alpha) {
    if (fabs(t) < 1e-8) return 1.0 - alpha + 4*alpha/M_PI;
    if (fabs(fabs(t) - T/(4*alpha)) < 1e-8) {
        double t1 = (1+2/M_PI)*sin(M_PI/(4*alpha));
        double t2 = (1-2/M_PI)*cos(M_PI/(4*alpha));
        return (alpha/sqrt(2))*(t1 + t2);
    }
    double num = sin(M_PI*t/T*(1-alpha)) + 4*alpha*t/T*cos(M_PI*t/T*(1+alpha));
    double den = M_PI*t/T*(1 - pow(4*alpha*t/T, 2));
    return num/den;
}

// 添加高斯白噪声
void add_awgn(float complex *signal, int len, double snr_db) {
    // 计算信号功率
    double signal_power = 0.0;
    for (int i = 0; i < len; i++) {
        signal_power += creal(signal[i])*creal(signal[i]) + cimag(signal[i])*cimag(signal[i]);
    }
    signal_power /= len;
    
    double snr_linear = pow(10.0, snr_db/10.0);
    double noise_power = signal_power / snr_linear;
    double noise_std = sqrt(noise_power/2);
    
    for (int i = 0; i < len; i++) {
        double u1, u2;
        do {
            u1 = rand() / (RAND_MAX + 1.0);
        } while (u1 < 1e-10);
        u2 = rand() / (RAND_MAX + 1.0);
        
        double z0 = sqrt(-2.0*log(u1)) * cos(2*M_PI*u2);
        double z1 = sqrt(-2.0*log(u1)) * sin(2*M_PI*u2);
        signal[i] += noise_std * (z0 + z1 * I);
    }
}

void generate_mqam(const char *filename, MQAMParams params) {
    // 计算参数
    int total_samples = params.num_symbols * params.samples_per_symbol;
    float complex *signal = calloc(total_samples, sizeof(float complex));
    float complex *symbols = malloc(params.num_symbols * sizeof(float complex));
    
    srand(time(NULL));
    double phase_rad = params.carrier_phase * M_PI / 180.0;
    double T = 1.0 / params.symbol_rate;
    double dt = T / params.samples_per_symbol;
    
    // 星座图生成
    int grid_size = sqrt(params.M);
    double d = 2.0 / (grid_size - 1);
    double scale = params.amplitude / sqrt(2.0*(params.M-1)/3.0);
    
    // 生成符号
    for (int i = 0; i < params.num_symbols; i++) {
        int symbol_idx = rand() % params.M;
        int x = symbol_idx % grid_size;
        int y = symbol_idx / grid_size;
        
        double data_I = (x - grid_size/2 + 0.5) * d;
        double data_Q = (y - grid_size/2 + 0.5) * d;
        symbols[i] = scale * (data_I + I * data_Q);
    }
    
    // 计算升余弦滤波器
    int span = 3;
    int filter_length = 2 * span * params.samples_per_symbol + 1;
    double *rrc_values = malloc(filter_length * sizeof(double));
    double t_step = T / params.samples_per_symbol;
    double energy = 0.0;
    
    // 生成滤波器系数
    for (int i = 0; i < filter_length; i++) {
        double t = -span*T + i * t_step;
        rrc_values[i] = rrc_pulse(t, T, params.rolloff);
        energy += rrc_values[i] * rrc_values[i];
    }
    
    // 能量归一化
    double norm_factor = sqrt(energy);
    for (int i = 0; i < filter_length; i++) {
        rrc_values[i] /= norm_factor;
    }
    
    // 脉冲成型和调制
    for (int i = 0; i < total_samples; i++) {
        double t = i * dt;
        float complex baseband = 0 + 0*I;
        
        if (params.pulse_shaping) {
            for (int j = -span; j <= span; j++) {
                int sym_idx = i/params.samples_per_symbol - j;
                if (sym_idx >= 0 && sym_idx < params.num_symbols) {
                    double tau = t - (sym_idx + 0.5) * T;
                    double pulse_val = rrc_pulse(tau, T, params.rolloff) / norm_factor;
                    baseband += symbols[sym_idx] * pulse_val;
                }
            }
        } else {
            int sym_idx = i / params.samples_per_symbol;
            if (sym_idx < params.num_symbols)
                baseband = symbols[sym_idx];
        }
        
        // 载波调制
        double carrier_phase = 2 * M_PI * params.carrier_freq * t + phase_rad;
        signal[i] = baseband * cexp(I * carrier_phase);
    }
    
    // 添加噪声
    if (params.snr_db > 0) {
        add_awgn(signal, total_samples, params.snr_db);
    }
    
    // 写入文件
    FILE *file = fopen(filename, "wb");
    fwrite(signal, sizeof(float complex), total_samples, file);
    fclose(file);
    
    free(signal);
    free(symbols);
    free(rrc_values);
}

int main() {
    MQAMParams params = {
        .M = 16,
        .num_symbols = 10000,
        .samples_per_symbol = 4,
        .carrier_freq = 10e3,
        .symbol_rate = 1e3,
        .snr_db = 20.0,
        .carrier_phase = 45.0,
        .amplitude = 0.8,
        .pulse_shaping = 1,
        .rolloff = 0.35
    };
    
    generate_mqam("16qam.bin", params);
    return 0;
}

Python验证脚本 (mqam_analyzer.py)

# -*- coding: utf-8 -*-
"""
M-QAM信号分析器
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal, fft

def analyze_mqam(filename, params):
    # 读取文件
    data = np.fromfile(filename, dtype=np.complex64)
    total_samples = len(data)
    sps = params['samples_per_symbol']
    num_symbols = params['num_symbols']
    sample_rate = params['symbol_rate'] * sps
    
    # 下变频
    t = np.arange(total_samples) / sample_rate
    carrier_phase = 2 * np.pi * params['carrier_freq'] * t + np.deg2rad(params['carrier_phase'])
    carrier = np.exp(-1j * carrier_phase)
    baseband = data * carrier
    
    # 匹配滤波
    span = 3
    if params['pulse_shaping']:
        t_symbol = np.linspace(-span, span, 2 * span * sps + 1)
        rrc = np.zeros(len(t_symbol))
        T_symbol = 1 / params['symbol_rate']
        for i, tau in enumerate(t_symbol):
            rrc[i] = rrc_pulse(tau, T_symbol, params['rolloff'])
        
        # 能量归一化
        rrc_norm = np.linalg.norm(rrc)
        rrc /= rrc_norm
        baseband = np.convolve(baseband, rrc, mode='same')
    
    # 符号抽取 - 精确在符号中心点抽取
    offset = span * sps  # 补偿滤波器延迟
    # 寻找最佳抽样点 (峰值附近)
    test_index = offset + sps // 2
    test_sample = baseband[test_index:test_index+100*sps:sps]
    avg_power = np.mean(np.abs(test_sample)**2)
    
    # 在符号中心点抽取
    start_index = offset + sps // 2
    valid_symbols = min(num_symbols, (total_samples - start_index) // sps)
    symbols = baseband[start_index:start_index + valid_symbols * sps:sps]
    
    # 计算理论平均功率 (用于SNR和EVM)
    ref_power = params['amplitude'] ** 2
    
    # 计算实际信号功率
    actual_power = np.mean(np.abs(symbols)**2)
    
    # 计算SNR
    ideal_symbols, evm_per_symbol = calculate_evm(symbols, params)
    noise = symbols - ideal_symbols
    noise_power = np.mean(np.abs(noise)**2)
    snr_measured = 10 * np.log10(ref_power / noise_power)  # 使用理论功率作为参考
    
    # 计算EVM (RMS)
    evm_rms = 100 * np.sqrt(noise_power / ref_power)
    
    # 星座图
    plt.figure(figsize=(15, 10))
    plt.subplot(2, 2, 1)
    plt.scatter(np.real(symbols), np.imag(symbols), alpha=0.6)
    plt.title(f'{params["M"]}-QAM Constellation\nMeasured SNR: {snr_measured:.2f} dB')
    plt.xlabel('I')
    plt.ylabel('Q')
    plt.grid(True)
    plt.axis('equal')
    
    # 眼图
    plt.subplot(2, 2, 2)
    eye_length = 2  # 2个符号周期
    num_eyes = min(200, valid_symbols - eye_length)
    for i in range(100, 100 + num_eyes):
        start_idx = offset + i * sps
        end_idx = start_idx + eye_length * sps
        segment = np.real(baseband[start_idx:end_idx])
        plt.plot(np.linspace(0, eye_length, len(segment)), segment, 'b-', alpha=0.1)
    plt.title('Eye Diagram (I Component)')
    plt.grid(True)
    
    # 频谱
    plt.subplot(2, 2, 3)
    nfft = 4096
    f = np.linspace(-sample_rate/2, sample_rate/2, nfft)
    spectrum = np.abs(fft.fftshift(fft.fft(data, nfft)))
    plt.plot(f, 20*np.log10(spectrum + 1e-10))
    plt.title('Power Spectral Density')
    plt.xlabel('Frequency (Hz)')
    plt.ylabel('dB')
    plt.grid(True)
    plt.ylim(-100, 50)
    
    # EVM分布
    plt.subplot(2, 2, 4)
    plt.hist(evm_per_symbol, bins=50, alpha=0.7)
    plt.title(f'EVM Distribution\nMean EVM: {np.mean(evm_per_symbol):.2f}% | RMS EVM: {evm_rms:.2f}%')
    plt.xlabel('EVM (%)')
    plt.ylabel('Count')
    plt.grid(True)
    
    plt.tight_layout()
    plt.savefig(f'{filename[:-4]}_analysis.png')
    plt.show()
    
    return {
        'measured_snr': snr_measured,
        'evm': evm_rms  # 使用RMS EVM作为结果
    }

def rrc_pulse(t, T, alpha):
    if abs(t) < 1e-8: 
        return 1 - alpha + 4*alpha/np.pi
    if abs(abs(t) - T/(4*alpha)) < 1e-8:
        t1 = (1+2/np.pi)*np.sin(np.pi/(4*alpha))
        t2 = (1-2/np.pi)*np.cos(np.pi/(4*alpha))
        return (alpha/np.sqrt(2))*(t1 + t2)
    num = np.sin(np.pi*t/T*(1-alpha)) + 4*alpha*t/T*np.cos(np.pi*t/T*(1+alpha))
    den = np.pi*t/T*(1 - (4*alpha*t/T)**2)
    return num/den

def calculate_evm(symbols, params):
    M = params['M']
    grid_size = int(np.sqrt(M))
    d = 2.0 / (grid_size - 1)
    scale = params['amplitude'] / np.sqrt(2*(M-1)/3)
    
    # 生成理想星座
    ideal_constellation = np.zeros(M, dtype=complex)
    for i in range(M):
        x = i % grid_size
        y = i // grid_size
        I = (x - grid_size/2 + 0.5) * d
        Q = (y - grid_size/2 + 0.5) * d
        ideal_constellation[i] = scale * (I + 1j*Q)
    
    # 计算每个符号的EVM (用于直方图)
    evm_per_symbol = []
    ideal_symbols = np.zeros_like(symbols)
    for i, sym in enumerate(symbols):
        distances = np.abs(ideal_constellation - sym)
        closest = np.argmin(distances)
        ideal_sym = ideal_constellation[closest]
        ideal_symbols[i] = ideal_sym
        error = sym - ideal_sym
        # 使用理论参考功率归一化
        evm_per_symbol.append(100 * np.abs(error) / np.abs(ideal_sym))
    
    return ideal_symbols, np.array(evm_per_symbol)

# 参数配置
params = {
    'M': 16,
    'num_symbols': 10000,
    'samples_per_symbol': 4,
    'carrier_freq': 10e3,
    'symbol_rate': 1e3,
    'carrier_phase': 45.0,
    'snr_db': 20.0,
    'amplitude': 0.8,
    'pulse_shaping': True,
    'rolloff': 0.35
}

# 计算采样率
params['sample_rate'] = params['symbol_rate'] * params['samples_per_symbol']

# 分析信号
results = analyze_mqam("16qam.bin", params)

# 结果验证
print(f"Measured SNR: {results['measured_snr']:.2f} dB (Expected: {params['snr_db']} dB)")
print(f"RMS EVM: {results['evm']:.2f}% (Expected: {100 * np.sqrt(10**(-params['snr_db']/10)):.2f}%)")

验证指标

1. EVM计算：
   $\text{EVM}=\frac{1}{N}{\sum }_{i=1}^N\frac{|S_i-{\hat{S}}_i|}{|{\hat{S}}_i|}\times 100\%$
2. 信噪比测量：
   ${\text{SNR}}_{\text{meas}}=10{\log }_{10}\left(\frac{P_{\text{signal}}}{P_{\text{noise}}}\right)$
   $P_{\text{noise}}=\frac{1}{N}\sum |S_i-{\hat{S}}_i{|}^2$

验证流程

1. 编译运行C程序：

gcc mqam_generator.c -o mqam_generator -lm
./mqam_generator

2. 生成文件：

• 16qam.bin：16-QAM调制信号（复数格式）

3. 运行Python分析：

python mqam_analyzer.py

4. 验证输出：

• 自动生成分析图表：

  • 星座图：显示信号在I-Q平面的分布
  • 眼图：评估信号质量及时序特性
  • 功率谱密度：显示信号频谱特性
  • EVM分布：误差向量幅度直方图
    

• 控制台输出关键指标：

    Measured SNR: 21.17 dB (Expected: 20.0 dB)
    RMS EVM: 8.74% (Expected: 10.00%)
  

5. 分析图表说明：

• 星座图：清晰显示16-QAM的4×4网格结构，点聚集程度反映噪声水平
• 眼图：张开的"眼睛"表示良好的信号完整性
• 频谱图：显示信号带宽和频谱效率
• EVM分布：评估调制质量，值越小越好

技术说明

1. 关键算法：

   • 星座映射：将比特流映射到复平面网格
   • 升余弦滤波：实现脉冲成型，减少码间干扰
   • AWGN生成：使用Box-Muller变换生成高斯噪声
   • EVM计算：误差向量幅度评估调制质量

2. 验证指标：

   • SNR测量：比较理论值与实际测量值
   • EVM分析：评估信号失真程度
   • 星座图完整性：检查符号聚类情况
   • 频谱特性：验证信号带宽符合预期

该实现完整覆盖了M-QAM信号生成、参数配置、数据存储和验证分析的全流程，可通过修改参数结构体生成不同阶数的QAM信号（4-QAM、64-QAM等）。

研究学习不易，点赞易。
工作生活不易，收藏易，点收藏不迷茫 ：）