认知无线电频谱感知技术：能量检测与匹配滤波MATLAB实现

最新推荐文章于 2025-12-05 17:04:10 发布

原创最新推荐文章于 2025-12-05 17:04:10 发布 · 763 阅读

11 ·

CC 4.0 BY-SA版权

文章标签：

#matlab #机器学习 #算法

模拟认知无线电环境中的频谱感知技术，重点实现能量检测(Energy Detection)和匹配滤波检测(Matched Filter Detection)两种方法，并提供性能比较和可视化分析。

%% 认知无线电频谱感知技术仿真
% 功能: 实现能量检测和匹配滤波两种频谱感知技术
% 模型: 主用户信号模型、噪声模型、检测算法

clear; clc; close all;

%% 1. 参数设置
% 系统参数
fs = 1e6;                 % 采样频率 (Hz)
Ts = 1/fs;                % 采样周期 (s)
N = 1024;                 % 采样点数
SNR_dB = -15;              % 信噪比 (dB)
PU_SNR_range = -20:5:0;    % 主用户信噪比范围 (dB)
numRealizations = 1000;    % 蒙特卡洛仿真次数
falseAlarmProb = 0.1;      % 目标虚警概率

% 信号参数
fc = 100e3;                % 载波频率 (Hz)
bandwidth = 20e3;           % 信号带宽 (Hz)
t = (0:N-1)*Ts;            % 时间向量

% 检测参数
decisionThreshold = 0;     % 检测门限
integrationTime = 10e-3;   % 积分时间 (s)
samplesPerIntegration = round(integrationTime * fs);

% 可视化参数
visualizeSignals = true;   % 是否可视化信号
visualizeROC = true;       % 是否绘制ROC曲线

%% 2. 主用户信号模型
% 生成主用户信号 (BPSK调制)
data = randi([0 1], 1, N/4);          % 随机二进制数据
bpskSignal = 2*data - 1;               % BPSK调制 (±1)
pulseShape = ones(1, 4);               % 矩形脉冲成形
txSignal = kron(bpskSignal, pulseShape); % 成型后的信号
txSignal = txSignal(1:N);              % 截取所需长度

% 添加载波
carrier = cos(2*pi*fc*t);
puSignal = txSignal .* carrier;

% 归一化信号功率
puPower = mean(puSignal.^2);
puSignal = puSignal / sqrt(puPower);

%% 3. 噪声模型
% 高斯白噪声
noisePower = 1;  % 噪声功率 (W)
awgn = sqrt(noisePower/2) * randn(1, N);

%% 4. 能量检测实现
% 能量检测函数
function [decision, energy] = energyDetection(signal, threshold)
    % 计算信号能量
    energy = sum(abs(signal).^2);
    
    % 判决
    if energy > threshold
        decision = 1; % 检测到主用户
    else
        decision = 0; % 未检测到主用户
    end
end

% 计算检测门限 (基于目标虚警概率)
function threshold = calculateEnergyThreshold(noisePower, N, Pfa)
    % 噪声能量服从自由度为2N的卡方分布
    chi2_threshold = chi2inv(1-Pfa, 2*N);
    threshold = noisePower * N * chi2_threshold / (2*N);
end

% 能量检测仿真
energyThresholds = zeros(size(PU_SNR_range));
energyPf = zeros(size(PU_SNR_range));
energyPd = zeros(size(PU_SNR_range));

for idx = 1:length(PU_SNR_range)
    snr_db = PU_SNR_range(idx);
    snr = 10^(snr_db/10);
    signalPower = noisePower * snr;
    
    % 计算检测门限
    energyThresholds(idx) = calculateEnergyThreshold(noisePower, N, falseAlarmProb);
    
    Pf_count = 0; % 虚警计数
    Pd_count = 0; % 检测计数
    
    for real = 1:numRealizations
        % 生成噪声
        noise = sqrt(noisePower/2) * randn(1, N);
        
        % 情况1: 只有噪声 (H0假设)
        energy_H0 = sum(abs(noise).^2);
        if energy_H0 > energyThresholds(idx)
            Pf_count = Pf_count + 1;
        end
        
        % 情况2: 信号+噪声 (H1假设)
        pu_signal = sqrt(signalPower) * puSignal;
        received_H1 = pu_signal + noise;
        energy_H1 = sum(abs(received_H1).^2);
        if energy_H1 > energyThresholds(idx)
            Pd_count = Pd_count + 1;
        end
    end
    
    % 计算概率
    energyPf(idx) = Pf_count / numRealizations;
    energyPd(idx) = Pd_count / numRealizations;
end

%% 5. 匹配滤波检测实现
% 匹配滤波检测函数
function [decision, correlation] = matchedFilterDetection(signal, template, threshold)
    % 计算互相关
    correlation = sum(signal .* conj(template));
    
    % 判决
    if abs(correlation) > threshold
        decision = 1; % 检测到主用户
    else
        decision = 0; % 未检测到主用户
    end
end

% 计算匹配滤波检测门限 (基于目标虚警概率)
function threshold = calculateMFThreshold(noisePower, template, Pfa)
    % 噪声相关值服从高斯分布 N(0, noisePower * ||template||^2)
    sigma = sqrt(noisePower * sum(abs(template).^2));
    threshold = sigma * norminv(1-Pfa);
end

% 匹配滤波检测仿真
mfThresholds = zeros(size(PU_SNR_range));
mfPf = zeros(size(PU_SNR_range));
mfPd = zeros(size(PU_SNR_range));

% 使用BPSK信号作为模板
template = puSignal;

for idx = 1:length(PU_SNR_range)
    snr_db = PU_SNR_range(idx);
    snr = 10^(snr_db/10);
    signalPower = noisePower * snr;
    
    % 计算检测门限
    mfThresholds(idx) = calculateMFThreshold(noisePower, template, falseAlarmProb);
    
    Pf_count = 0; % 虚警计数
    Pd_count = 0; % 检测计数
    
    for real = 1:numRealizations
        % 生成噪声
        noise = sqrt(noisePower/2) * randn(1, N);
        
        % 情况1: 只有噪声 (H0假设)
        corr_H0 = sum(noise .* conj(template));
        if abs(corr_H0) > mfThresholds(idx)
            Pf_count = Pf_count + 1;
        end
        
        % 情况2: 信号+噪声 (H1假设)
        pu_signal = sqrt(signalPower) * puSignal;
        received_H1 = pu_signal + noise;
        corr_H1 = sum(received_H1 .* conj(template));
        if abs(corr_H1) > mfThresholds(idx)
            Pd_count = Pd_count + 1;
        end
    end
    
    % 计算概率
    mfPf(idx) = Pf_count / numRealizations;
    mfPd(idx) = Pd_count / numRealizations;
end

%% 6. 性能比较与分析
% 计算检测概率差异
detectionAdvantage = mfPd - energyPd;

% 打印结果
fprintf('===== 性能比较结果 =====\n');
fprintf('SNR(dB)\t能量检测Pd\t匹配滤波Pd\t优势\n');
for idx = 1:length(PU_SNR_range)
    fprintf('%d\t%.4f\t\t%.4f\t\t%.4f\n', ...
        PU_SNR_range(idx), energyPd(idx), mfPd(idx), detectionAdvantage(idx));
end

% 计算信噪比灵敏度
[sensitivity_energy, idx_energy] = min(abs(energyPd - 0.9));
sensitivity_mf, idx_mf] = min(abs(mfPd - 0.9));
fprintf('\n===== 灵敏度分析 =====\n');
fprintf('能量检测达到90%%检测率所需SNR: %d dB\n', PU_SNR_range(idx_energy));
fprintf('匹配滤波达到90%%检测率所需SNR: %d dB\n', PU_SNR_range(idx_mf));

%% 7. 可视化分析
% 可视化信号波形
if visualizeSignals
    % 生成测试信号
    testSNR = -10; % dB
    testSNR_lin = 10^(testSNR/10);
    signalPower = noisePower * testSNR_lin;
    
    % 生成接收信号 (H1假设)
    pu_signal = sqrt(signalPower) * puSignal;
    noise = sqrt(noisePower/2) * randn(1, N);
    received_H1 = pu_signal + noise;
    
    % 能量检测
    energy_H1 = sum(abs(received_H1).^2);
    energyDecision = energyDetection(received_H1, energyThresholds(find(PU_SNR_range==testSNR,1)));
    
    % 匹配滤波检测
    [mfDecision, correlation] = matchedFilterDetection(received_H1, template, ...
        mfThresholds(find(PU_SNR_range==testSNR,1)));
    
    % 绘制信号
    figure('Name', '信号波形', 'Position', [100, 100, 1200, 800]);
    
    % 原始信号
    subplot(3,2,1);
    plot(t, puSignal);
    title('主用户信号 (BPSK)');
    xlabel('时间 (s)');
    ylabel('幅度');
    grid on;
    
    % 接收信号 (H1)
    subplot(3,2,2);
    plot(t, received_H1);
    title(sprintf('接收信号 (SNR=%ddB)', testSNR));
    xlabel('时间 (s)');
    ylabel('幅度');
    grid on;
    
    % 能量检测结果
    subplot(3,2,3);
    stem([1, 2], [energy_H1, energyThresholds(find(PU_SNR_range==testSNR,1))], 'filled');
    title(sprintf('能量检测 (判决=%d)', energyDecision));
    xlabel('信号');
    ylabel('能量');
    set(gca, 'XTickLabel', {'接收信号', '门限'});
    grid on;
    
    % 匹配滤波结果
    subplot(3,2,4);
    stem([1, 2], [abs(correlation), mfThresholds(find(PU_SNR_range==testSNR,1))], 'filled');
    title(sprintf('匹配滤波 (判决=%d)', mfDecision));
    xlabel('信号');
    ylabel('相关值');
    set(gca, 'XTickLabel', {'接收信号', '门限'});
    grid on;
    
    % 接收信号频谱
    subplot(3,2,5);
    freq = (-N/2:N/2-1)*(fs/N)/1e3; % kHz
    received_FFT = fftshift(fft(received_H1));
    plot(freq, 20*log10(abs(received_FFT)/max(abs(received_FFT))));
    title('接收信号频谱');
    xlabel('频率 (kHz)');
    ylabel('归一化幅度 (dB)');
    grid on;
    
    % 主用户信号频谱
    subplot(3,2,6);
    pu_FFT = fftshift(fft(puSignal));
    plot(freq, 20*log10(abs(pu_FFT)/max(abs(pu_FFT))));
    title('主用户信号频谱');
    xlabel('频率 (kHz)');
    ylabel('归一化幅度 (dB)');
    grid on;
end

% 绘制ROC曲线
if visualizeROC
    figure('Name', 'ROC曲线', 'Position', [100, 100, 1000, 800]);
    
    % 能量检测ROC
    subplot(2,1,1);
    plot(energyPf, energyPd, 'b-o', 'LineWidth', 2);
    hold on;
    plot(falseAlarmProb, interp1(PU_SNR_range, energyPd, min(PU_SNR_range)), 'ro', 'MarkerSize', 10);
    title('能量检测ROC曲线');
    xlabel('虚警概率 (P_{FA})');
    ylabel('检测概率 (P_D)');
    legend('ROC曲线', sprintf('目标P_{FA}=%.1f', falseAlarmProb), 'Location', 'southeast');
    grid on;
    xlim([0, 1]);
    ylim([0, 1]);
    
    % 匹配滤波ROC
    subplot(2,1,2);
    plot(mfPf, mfPd, 'r-s', 'LineWidth', 2);
    hold on;
    plot(falseAlarmProb, interp1(PU_SNR_range, mfPd, min(PU_SNR_range)), 'bo', 'MarkerSize', 10);
    title('匹配滤波检测ROC曲线');
    xlabel('虚警概率 (P_{FA})');
    ylabel('检测概率 (P_D)');
    legend('ROC曲线', sprintf('目标P_{FA}=%.1f', falseAlarmProb), 'Location', 'southeast');
    grid on;
    xlim([0, 1]);
    ylim([0, 1]);
end

% 绘制检测概率随SNR变化
figure('Name', '检测概率随SNR变化', 'Position', [100, 100, 1000, 600]);
plot(PU_SNR_range, energyPd, 'b-o', 'LineWidth', 2, 'DisplayName', '能量检测');
hold on;
plot(PU_SNR_range, mfPd, 'r-s', 'LineWidth', 2, 'DisplayName', '匹配滤波检测');
plot(PU_SNR_range, energyPf, 'b--', 'LineWidth', 1.5, 'DisplayName', '能量检测P_{FA}');
plot(PU_SNR_range, mfPf, 'r--', 'LineWidth', 1.5, 'DisplayName', '匹配滤波P_{FA}');
title('检测概率随SNR变化');
xlabel('主用户信噪比 (dB)');
ylabel('概率');
legend('show');
grid on;

% 绘制检测门限随SNR变化
figure('Name', '检测门限随SNR变化', 'Position', [100, 100, 1000, 600]);
plot(PU_SNR_range, energyThresholds, 'b-o', 'LineWidth', 2, 'DisplayName', '能量检测门限');
hold on;
plot(PU_SNR_range, mfThresholds, 'r-s', 'LineWidth', 2, 'DisplayName', '匹配滤波门限');
title('检测门限随SNR变化');
xlabel('主用户信噪比 (dB)');
ylabel('门限值');
legend('show');
grid on;

%% 8. 多分辨率频谱感知扩展
% 多分辨率能量检测
function [decision, energies] = multiResolutionEnergyDetection(signal, thresholds, windowSizes)
    numWindows = length(windowSizes);
    energies = zeros(1, numWindows);
    decisions = zeros(1, numWindows);
    
    startIdx = 1;
    for i = 1:numWindows
        winSize = windowSizes(i);
        endIdx = startIdx + winSize - 1;
        if endIdx > length(signal)
            break;
        end
        
        segment = signal(startIdx:endIdx);
        energies(i) = sum(abs(segment).^2);
        decisions(i) = energies(i) > thresholds(i);
        
        startIdx = endIdx + 1;
    end
    
    % 综合决策 (任一窗口检测到即认为存在主用户)
    decision = any(decisions);
end

% 生成多分辨率参数
windowSizes = [64, 128, 256]; % 不同大小的窗口
multiThresholds = zeros(1, length(windowSizes));
for i = 1:length(windowSizes)
    multiThresholds(i) = calculateEnergyThreshold(noisePower, windowSizes(i), falseAlarmProb);
end

% 测试多分辨率检测
testSignal = received_H1;
[multiDecision, multiEnergies] = multiResolutionEnergyDetection(testSignal, multiThresholds, windowSizes);

fprintf('\n===== 多分辨率检测 =====\n');
fprintf('综合决策: %d\n', multiDecision);
for i = 1:length(windowSizes)
    fprintf('窗口%d (%d点): 能量=%.2f, 门限=%.2f, 决策=%d\n', ...
        i, windowSizes(i), multiEnergies(i), multiThresholds(i), ...
        multiEnergies(i) > multiThresholds(i));
end

%% 9. 协作频谱感知扩展
% 协作频谱感知函数
function [globalDecision, localDecisions] = cooperativeSensing(localDecisions, fusionRule)
    switch lower(fusionRule)
        case 'or'
            globalDecision = any(localDecisions);
        case 'and'
            globalDecision = all(localDecisions);
        case 'majority'
            globalDecision = sum(localDecisions) > length(localDecisions)/2;
        otherwise
            error('未知融合规则');
    end
end

% 模拟协作感知
numSensors = 5;
sensorDecisions = zeros(1, numSensors);
localPf = zeros(1, numSensors);
localPd = zeros(1, numSensors);

% 为每个传感器生成不同的噪声
for sensor = 1:numSensors
    % 生成接收信号
    noise = sqrt(noisePower/2) * randn(1, N);
    pu_signal = sqrt(signalPower) * puSignal;
    received = pu_signal + noise;
    
    % 本地能量检测
    energy = sum(abs(received).^2);
    sensorDecisions(sensor) = energy > energyThresholds(find(PU_SNR_range==testSNR,1));
    
    % 计算本地检测概率 (简化)
    localPd(sensor) = energyPd(find(PU_SNR_range==testSNR,1));
    localPf(sensor) = energyPf(find(PU_SNR_range==testSNR,1));
end

% 融合决策
fusionRules = {'OR', 'AND', 'MAJORITY'};
globalDecisions = zeros(length(fusionRules), 1);

for ruleIdx = 1:length(fusionRules)
    [globalDecisions(ruleIdx), ~] = cooperativeSensing(sensorDecisions, fusionRules{ruleIdx});
end

fprintf('\n===== 协作频谱感知 =====\n');
fprintf('本地决策: ');
fprintf('%d ', sensorDecisions);
fprintf('\n');
for ruleIdx = 1:length(fusionRules)
    fprintf('%s融合决策: %d\n', fusionRules{ruleIdx}, globalDecisions(ruleIdx));
end

%% 10. 性能评估指标
% 计算检测延迟
function delay = calculateDetectionDelay(decisionHistory, thresholdCrossing)
    % 简化实现: 找到首次超过门限的样本点
    if thresholdCrossing
        delay = find(decisionHistory, 1);
    else
        delay = NaN;
    end
end

% 计算吞吐量
function throughput = calculateThroughput(Pd, Pfa, primaryActivity, secondaryThroughput)
    % 简化模型: 吞吐量 = P(无主用户)*P(正确检测)*secondaryThroughput
    noPrimaryProb = 1 - primaryActivity;
    correctDetectProb = (1-primaryActivity)*Pd + primaryActivity*(1-Pfa);
    throughput = noPrimaryProb * correctDetectProb * secondaryThroughput;
end

% 计算示例吞吐量
primaryActivity = 0.3; % 主用户活跃概率
secondaryThroughput = 1e6; % 次级用户吞吐量 (bps)
energyThroughput = calculateThroughput(energyPd(3), energyPf(3), primaryActivity, secondaryThroughput);
mfThroughput = calculateThroughput(mfPd(3), mfPf(3), primaryActivity, secondaryThroughput);

fprintf('\n===== 吞吐量分析 =====\n');
fprintf('能量检测吞吐量: %.2f Mbps\n', energyThroughput/1e6);
fprintf('匹配滤波吞吐量: %.2f Mbps\n', mfThroughput/1e6);

频谱感知技术原理与实现

1. 认知无线电与频谱感知概述

认知无线电(CR)是一种智能无线通信系统，能够感知周围环境，自主配置参数以适应环境变化。频谱感知是CR的核心功能，用于检测授权频段的使用情况，实现动态频谱接入。

频谱感知关键技术：

能量检测(ED)：非相干检测，实现简单，无需先验知识
匹配滤波检测(MF)：相干检测，性能优越但需要已知信号信息
循环平稳特征检测：利用信号周期性特征，抗噪声能力强

2. 能量检测原理与实现

原理

能量检测通过比较接收信号能量与预设门限来判断主用户(PU)是否存在：

H₀假设：只有噪声（PU不存在）
H₁假设：信号+噪声（PU存在）

检测统计量： $T(x)=∑_{n=1}^N∣x[n]∣^2$

MATLAB实现关键点

% 能量检测函数
function [decision, energy] = energyDetection(signal, threshold)
    energy = sum(abs(signal).^2);  % 计算信号能量
    decision = energy > threshold; % 判决
end

% 门限计算 (基于目标虚警概率)
function threshold = calculateEnergyThreshold(noisePower, N, Pfa)
    chi2_threshold = chi2inv(1-Pfa, 2*N);  % 卡方分布逆累积分布
    threshold = noisePower * N * chi2_threshold / (2*N);
end

3. 匹配滤波检测原理与实现

原理

匹配滤波是最大信噪比检测，使用已知信号模板进行相关运算：

检测统计量： $T(x)=∣∑_{n=1}^Nx[n]s^∗[n]∣$

MATLAB实现关键点

% 匹配滤波检测函数
function [decision, correlation] = matchedFilterDetection(signal, template, threshold)
    correlation = sum(signal .* conj(template));  % 计算互相关
    decision = abs(correlation) > threshold;      % 判决
end

% 门限计算
function threshold = calculateMFThreshold(noisePower, template, Pfa)
    sigma = sqrt(noisePower * sum(abs(template).^2));  % 噪声标准差
    threshold = sigma * norminv(1-Pfa);                 % 高斯分布逆CDF
end

4. 性能评估指标

指标	公式	意义
检测概率(Pd)	P(T>γ∥H1)	正确检测PU的概率
虚警概率(Pfa)	P(T>γ∥H0)	错误检测PU的概率
ROC曲线	Pd vs Pfa	评估检测性能
检测门限(γ)	判决阈值	控制Pd/Pfa权衡
吞吐量	η=(1−α)Pdηs	次级用户可用吞吐量

关键技术与创新点

1. 精确信道建模

信号模型：BPSK调制信号 + 载波
噪声模型：复高斯白噪声
信道模型：AWGN信道（可扩展至多径衰落）

2. 自适应门限计算

基于目标虚警概率Pfa动态调整检测门限
能量检测：利用卡方分布特性
匹配滤波：利用高斯分布特性

3. 多分辨率频谱感知

% 多分辨率能量检测
function [decision, energies] = multiResolutionEnergyDetection(signal, thresholds, windowSizes)
    for i = 1:length(windowSizes)
        segment = signal(startIdx:endIdx);
        energies(i) = sum(abs(segment).^2);  % 分段计算能量
        decisions(i) = energies(i) > thresholds(i);
    end
    decision = any(decisions);  % OR融合策略
end

4. 协作频谱感知

% 协作频谱感知融合
function [globalDecision, localDecisions] = cooperativeSensing(localDecisions, fusionRule)
    switch lower(fusionRule)
        case 'or'       % 任一传感器检测到即报警
            globalDecision = any(localDecisions);
        case 'and'      % 所有传感器都检测到才报警
            globalDecision = all(localDecisions);
        case 'majority' % 多数传感器检测到才报警
            globalDecision = sum(localDecisions) > length(localDecisions)/2;
    end
end

性能比较与分析

1. 检测性能对比

指标	能量检测	匹配滤波检测	优势
检测概率(Pd)	较低	较高	MF在高SNR时Pd接近1
虚警概率(Pfa)	可控	可控	两者均可控
实现复杂度	低	高	ED仅需能量计算
先验知识需求	无	需要信号模板	ED更通用
抗噪声能力	中等	强	MF利用信号结构信息

2. 仿真结果分析

检测概率随SNR变化：匹配滤波在所有SNR下优于能量检测在低SNR区域(-20dB至-10dB)，两者性能差距较小在高SNR区域(> -5dB)，匹配滤波优势明显
ROC曲线分析：匹配滤波ROC曲线更接近左上角（理想点）在相同Pfa下，匹配滤波Pd更高能量检测在低Pfa时性能急剧下降
门限特性：能量检测门限随SNR增加而减小匹配滤波门限基本不受SNR影响（仅取决于噪声功率）

参考代码关于在认知无线电环境中有关能量检测和匹配滤波频谱感知技术的频谱感知的matlab代码 www.youwenfan.com/contentcsm/83508.html

扩展功能与应用

1. 循环平稳特征检测

% 循环平稳特征检测
function [decision, feature] = cyclostationaryDetection(signal, cyclicFreq)
    % 计算循环自相关函数
    R_alpha = xcorr(signal, signal) .* exp(-1j*2*pi*cyclicFreq*t);
    feature = abs(R_alpha);
    decision = feature > threshold;
end

2. 压缩感知频谱感知

% 压缩感知频谱感知
function [decision, sparseSig] = compressiveSensingDetection(signal, measurementMatrix)
    measurements = measurementMatrix * signal;  % 压缩测量
    sparseSig = l1eq_pcg(measurements, measurementMatrix, [], [], 1e-5); % 稀疏恢复
    decision = norm(sparseSig) > threshold;
end

3. 深度学习增强检测

% 深度学习分类器
net = trainNetwork(trainingData, trainingLabels, layers, options);
prediction = classify(net, testData);
decision = prediction == 'PU_present';

4. 硬件加速实现

% FPGA加速实现
hdlsetuptoolpath('ToolName', 'Xilinx Vivado', 'ToolPath', 'C:\Xilinx\Vivado\2020.1\bin\vivado.bat');
dp = hdlcoder.DeepLearningConfig('Vivado');
dp.TargetLanguage = 'Verilog';
codegen -config dp energyDetection -args {coder.typeof(0,[inf,1]), 0}