直扩信号参数估计：载频、码速率和扩频增益

MATLAB程序用于估计接收到的直扩信号的载频、码速率和扩频增益：

主程序文件

%% 直扩信号参数估计：载频、码速率和扩频增益
clear; close all; clc;

fprintf('=== 直扩信号参数估计 ===\n');

%% 生成测试直扩信号（如果无真实信号）
[receivedSignal, params] = generateDSSSTestSignal();

%% 信号预处理
[preprocessedSignal, preprocessingInfo] = preprocessSignal(receivedSignal, params);

%% 载频估计
[carrierFreq, carrierEstimationInfo] = estimateCarrierFrequency(preprocessedSignal, params);

%% 码速率估计
[codeRate, codeEstimationInfo] = estimateCodeRate(preprocessedSignal, params);

%% 扩频增益估计
[spreadingGain, gainEstimationInfo] = estimateSpreadingGain(preprocessedSignal, codeRate, params);

%% 参数验证和精度分析
[validationResults, confidence] = validateEstimations(carrierFreq, codeRate, spreadingGain, params);

%% 结果显示和可视化
displayResults(carrierFreq, codeRate, spreadingGain, validationResults, confidence, params);

%% 性能分析
analyzePerformance(carrierFreq, codeRate, spreadingGain, params);

fprintf('\n直扩信号参数估计完成！\n');

测试信号生成

function [receivedSignal, params] = generateDSSSTestSignal()
    % 生成直扩测试信号
    
    fprintf('生成直扩测试信号...\n');
    
    % 信号参数
    params = struct();
    params.true.carrierFreq = 2.4e9;          % 真实载频 2.4 GHz
    params.true.codeRate = 10e6;              % 真实码速率 10 MHz
    params.true.spreadingGain = 127;          % 真实扩频增益 127
    params.true.samplingRate = 100e6;         % 采样率 100 MHz
    params.true.signalDuration = 100e-6;      % 信号持续时间 100 μs
    
    % 信道参数
    params.channel.SNR = 10;                  % 信噪比 dB
    params.channel.frequencyOffset = 50e3;    % 频率偏移 50 kHz
    params.channel.phaseNoise = 0.01;         % 相位噪声
    
    % 计算信号参数
    numSamples = round(params.true.signalDuration * params.true.samplingRate);
    t = (0:numSamples-1) / params.true.samplingRate;
    
    % 生成扩频码 (Gold码或m序列)
    spreadingCode = generateSpreadingCode(params.true.spreadingGain);
    
    % 生成信息序列
    infoBits = randi([0 1], 1, ceil(numSamples / params.true.spreadingGain));
    
    % 扩频调制
    basebandSignal = spreadSignal(infoBits, spreadingCode, params);
    
    % 载波调制
    carrier = exp(1j * 2 * pi * params.true.carrierFreq * t);
    modulatedSignal = basebandSignal .* carrier;
    
    % 添加频率偏移和相位噪声
    frequencyOffset = exp(1j * 2 * pi * params.channel.frequencyOffset * t);
    phaseNoise = exp(1j * params.channel.phaseNoise * randn(size(t)));
    
    modulatedSignal = modulatedSignal .* frequencyOffset .* phaseNoise;
    
    % 添加噪声
    receivedSignal = addNoise(modulatedSignal, params.channel.SNR);
    
    fprintf('测试信号生成完成:\n');
    fprintf('  真实载频: %.3f GHz\n', params.true.carrierFreq/1e9);
    fprintf('  真实码速率: %.3f MHz\n', params.true.codeRate/1e6);
    fprintf('  真实扩频增益: %d\n', params.true.spreadingGain);
    fprintf('  信号长度: %d 采样点\n', numSamples);
end

function spreadingCode = generateSpreadingCode(gain)
    % 生成扩频码 (m序列)
    
    % 使用线性反馈移位寄存器生成m序列
    registerLength = ceil(log2(gain + 1));
    
    % 选择本原多项式 (示例)
    switch registerLength
        case 7  % 127位m序列
            taps = [7, 1];  % x^7 + x + 1
        case 6  % 63位m序列
            taps = [6, 1];  % x^6 + x + 1
        case 5  % 31位m序列
            taps = [5, 2];  % x^5 + x^2 + 1
        otherwise
            taps = [registerLength, 1];  % 默认
    end
    
    % 初始化寄存器
    register = ones(1, registerLength);
    spreadingCode = zeros(1, 2^registerLength - 1);
    
    for i = 1:length(spreadingCode)
        % 输出位
        spreadingCode(i) = register(end);
        
        % 计算反馈
        feedback = mod(sum(register(taps)), 2);
        
        % 移位
        register = [feedback, register(1:end-1)];
    end
    
    % 转换为BPSK调制 (-1, +1)
    spreadingCode = 2 * spreadingCode - 1;
    
    % 截断到指定增益
    spreadingCode = spreadingCode(1:gain);
end

function spreadSignal = spreadSignal(infoBits, spreadingCode, params)
    % 扩频调制
    
    numChipsPerBit = length(spreadingCode);
    numBits = length(infoBits);
    
    % 将信息比特映射到BPSK
    mappedBits = 2 * infoBits - 1;
    
    % 扩频
    spreadSignal = zeros(1, numBits * numChipsPerBit);
    for i = 1:numBits
        startIdx = (i-1) * numChipsPerBit + 1;
        endIdx = i * numChipsPerBit;
        spreadSignal(startIdx:endIdx) = mappedBits(i) * spreadingCode;
    end
    
    % 调整采样率
    samplesPerChip = round(params.true.samplingRate / params.true.codeRate);
    if samplesPerChip > 1
        % 上采样
        spreadSignal = repelem(spreadSignal, samplesPerChip);
    end
end

function noisySignal = addNoise(signal, SNR_dB)
    % 添加高斯白噪声
    
    signalPower = mean(abs(signal).^2);
    noisePower = signalPower / (10^(SNR_dB/10));
    noise = sqrt(noisePower/2) * (randn(size(signal)) + 1j*randn(size(signal)));
    
    noisySignal = signal + noise;
end

信号预处理

function [preprocessedSignal, info] = preprocessSignal(receivedSignal, params)
    % 信号预处理
    
    fprintf('信号预处理...\n');
    
    info = struct();
    
    % 1. 自动增益控制
    preprocessedSignal = automaticGainControl(receivedSignal);
    
    % 2. 带通滤波（如果知道大致频带）
    if isfield(params, 'estimatedCarrierRange')
        preprocessedSignal = bandpassFilter(preprocessedSignal, params);
    end
    
    % 3. 下变频到基带（使用估计的载频或已知范围）
    if isfield(params, 'estimatedCarrierFreq')
        preprocessedSignal = downconvertToBaseband(preprocessedSignal, params);
    end
    
    % 4. 计算信号统计特性
    info.signalPower = mean(abs(preprocessedSignal).^2);
    info.noisePower = estimateNoisePower(preprocessedSignal);
    info.estimatedSNR = 10 * log10(info.signalPower / info.noisePower);
    
    fprintf('  信号功率: %.2f dB\n', 10*log10(info.signalPower));
    fprintf('  噪声功率: %.2f dB\n', 10*log10(info.noisePower));
    fprintf('  估计SNR: %.2f dB\n', info.estimatedSNR);
    
    fprintf('信号预处理完成\n');
end

function processedSignal = automaticGainControl(signal)
    % 自动增益控制
    
    targetPower = 1.0;
    currentPower = mean(abs(signal).^2);
    gain = sqrt(targetPower / currentPower);
    
    processedSignal = signal * gain;
end

function filteredSignal = bandpassFilter(signal, params)
    % 带通滤波
    
    samplingRate = params.true.samplingRate;
    
    if isfield(params, 'estimatedCarrierRange')
        % 使用估计的载频范围
        f_low = params.estimatedCarrierRange(1);
        f_high = params.estimatedCarrierRange(2);
    else
        % 默认频带
        f_low = 2.3e9;
        f_high = 2.5e9;
    end
    
    % 设计带通滤波器
    nyquist = samplingRate / 2;
    Wn = [f_low, f_high] / nyquist;
    
    if Wn(2) >= 1
        Wn(2) = 0.99;
    end
    
    [b, a] = butter(6, Wn, 'bandpass');
    filteredSignal = filter(b, a, signal);
end

function basebandSignal = downconvertToBaseband(signal, params)
    % 下变频到基带
    
    samplingRate = params.true.samplingRate;
    t = (0:length(signal)-1) / samplingRate;
    
    if isfield(params, 'estimatedCarrierFreq')
        carrierFreq = params.estimatedCarrierFreq;
    else
        carrierFreq = params.true.carrierFreq;  % 使用真实值（实际中未知）
    end
    
    % 生成本地载波
    localCarrier = exp(-1j * 2 * pi * carrierFreq * t);
    
    % 下变频
    basebandSignal = signal .* localCarrier;
    
    % 低通滤波
    cutoffFreq = min(10e6, samplingRate/4);  % 假设信号带宽小于10MHz
    Wn = cutoffFreq / (samplingRate/2);
    [b, a] = butter(6, Wn, 'low');
    basebandSignal = filter(b, a, basebandSignal);
end

function noisePower = estimateNoisePower(signal)
    % 估计噪声功率
    
    % 使用小波变换估计噪声
    [c, l] = wavedec(real(signal), 5, 'db4');
    noisePower = median(abs(c)) / 0.6745;
    noisePower = noisePower^2;
end

载频估计

function [carrierFreq, info] = estimateCarrierFrequency(signal, params)
    % 估计载波频率
    
    fprintf('估计载波频率...\n');
    
    info = struct();
    samplingRate = params.true.samplingRate;
    
    % 方法1: 功率谱密度峰值检测
    [carrierFreq1, confidence1] = estimateByPSD(signal, samplingRate);
    
    % 方法2: 循环平稳特征检测
    [carrierFreq2, confidence2] = estimateByCyclostationarity(signal, samplingRate);
    
    % 方法3: 高阶统计量
    [carrierFreq3, confidence3] = estimateByHigherOrderStats(signal, samplingRate);
    
    % 融合估计结果
    carrierFreq = fuseEstimations(carrierFreq1, carrierFreq2, carrierFreq3, ...
                                 confidence1, confidence2, confidence3);
    
    info.methods.psd = carrierFreq1;
    info.methods.cyclo = carrierFreq2;
    info.methods.hos = carrierFreq3;
    info.confidence = [confidence1, confidence2, confidence3];
    info.finalEstimation = carrierFreq;
    
    % 误差分析
    if isfield(params.true, 'carrierFreq')
        trueFreq = params.true.carrierFreq;
        error = abs(carrierFreq - trueFreq);
        relativeError = error / trueFreq * 100;
        
        info.trueValue = trueFreq;
        info.absoluteError = error;
        info.relativeError = relativeError;
        
        fprintf('  真实载频: %.6f GHz\n', trueFreq/1e9);
        fprintf('  估计载频: %.6f GHz\n', carrierFreq/1e9);
        fprintf('  绝对误差: %.3f kHz\n', error/1e3);
        fprintf('  相对误差: %.3f%%\n', relativeError);
    else
        fprintf('  估计载频: %.6f GHz\n', carrierFreq/1e9);
    end
    
    fprintf('载频估计完成\n');
end

function [carrierFreq, confidence] = estimateByPSD(signal, samplingRate)
    % 通过功率谱密度估计载频
    
    % 计算功率谱密度
    [psd, f] = pwelch(signal, hamming(1024), 512, 1024, samplingRate);
    
    % 寻找峰值
    [peaks, peakLocs] = findpeaks(psd, 'SortStr', 'descend', 'NPeaks', 5);
    
    if isempty(peaks)
        carrierFreq = 0;
        confidence = 0;
        return;
    end
    
    % 选择最高峰值对应的频率
    mainPeakLoc = peakLocs(1);
    carrierFreq = f(mainPeakLoc);
    
    % 计算置信度（基于峰值突出程度）
    if length(peaks) > 1
        prominence = (peaks(1) - peaks(2)) / peaks(1);
        confidence = min(1, prominence * 2);
    else
        confidence = 0.7;
    end
    
    % 可视化
    if false  % 设置为true来显示PSD
        figure;
        plot(f/1e6, 10*log10(psd));
        xlabel('频率 (MHz)');
        ylabel('功率谱密度 (dB/Hz)');
        title('功率谱密度分析');
        grid on;
        hold on;
        plot(carrierFreq/1e6, 10*log10(peaks(1)), 'ro', 'MarkerSize', 10);
        legend('PSD', '估计载频');
    end
end

function [carrierFreq, confidence] = estimateByCyclostationarity(signal, samplingRate)
    % 通过循环平稳特征估计载频
    
    % 计算循环谱
    alpha = 0:0.01:0.5;  % 循环频率范围
    cyclicSpectrum = zeros(length(alpha), 1024);
    
    for i = 1:length(alpha)
        % 简化的循环谱计算
        shiftedSignal = signal .* exp(-1j * 2 * pi * alpha(i) * (0:length(signal)-1));
        [psd, ~] = pwelch(shiftedSignal, hamming(256), 128, 256, samplingRate);
        cyclicSpectrum(i, :) = psd(1:1024);
    end
    
    % 寻找循环谱中的特征
    featureMap = sum(cyclicSpectrum, 2);
    [~, maxIdx] = max(featureMap);
    
    % 估计载频（基于循环频率）
    carrierFreq = alpha(maxIdx) * samplingRate;
    
    % 置信度基于特征突出程度
    sortedFeatures = sort(featureMap, 'descend');
    if length(sortedFeatures) > 1
        confidence = (sortedFeatures(1) - sortedFeatures(2)) / sortedFeatures(1);
    else
        confidence = 0.6;
    end
    
    confidence = min(0.9, confidence);
end

function [carrierFreq, confidence] = estimateByHigherOrderStats(signal, samplingRate)
    % 通过高阶统计量估计载频
    
    % 计算四阶累积量（对载频敏感）
    signalReal = real(signal);
    signalImag = imag(signal);
    
    % 简化的高阶统计量分析
    kurtosisReal = kurtosis(signalReal);
    kurtosisImag = kurtosis(signalImag);
    
    % 基于峰度的频率估计（启发式方法）
    freqShift = (kurtosisReal - 3) * 1e6;  % 缩放因子
    
    % 初始频率估计
    [psd, f] = pwelch(signal, hamming(512), 256, 512, samplingRate);
    [~, initialFreqIdx] = max(psd);
    initialFreq = f(initialFreqIdx);
    
    carrierFreq = initialFreq + freqShift;
    
    % 置信度基于峰度值
    confidence = min(0.8, abs(kurtosisReal - 3) / 10);
end

function fusedFreq = fuseEstimations(freq1, freq2, freq3, conf1, conf2, conf3)
    % 融合多个估计结果
    
    estimations = [freq1, freq2, freq3];
    confidences = [conf1, conf2, conf3];
    
    % 归一化置信度权重
    totalConfidence = sum(confidences);
    if totalConfidence == 0
        weights = [1/3, 1/3, 1/3];
    else
        weights = confidences / totalConfidence;
    end
    
    % 加权平均
    fusedFreq = sum(estimations .* weights);
    
    % 如果估计值差异太大，选择最可信的
    stdEstimations = std(estimations);
    if stdEstimations > 0.1 * mean(estimations)  % 差异超过10%
        [~, bestIdx] = max(confidences);
        fusedFreq = estimations(bestIdx);
    end
end

码速率估计

function [codeRate, info] = estimateCodeRate(signal, params)
    % 估计码速率
    
    fprintf('估计码速率...\n');
    
    info = struct();
    samplingRate = params.true.samplingRate;
    
    % 方法1: 自相关函数分析
    [codeRate1, confidence1] = estimateByAutocorrelation(signal, samplingRate);
    
    % 方法2: 小波变换分析
    [codeRate2, confidence2] = estimateByWavelet(signal, samplingRate);
    
    % 方法3: 循环自相关
    [codeRate3, confidence3] = estimateByCyclicAutocorrelation(signal, samplingRate);
    
    % 融合估计结果
    codeRate = fuseCodeRateEstimations(codeRate1, codeRate2, codeRate3, ...
                                      confidence1, confidence2, confidence3);
    
    info.methods.autocorr = codeRate1;
    info.methods.wavelet = codeRate2;
    info.methods.cyclic = codeRate3;
    info.confidence = [confidence1, confidence2, confidence3];
    info.finalEstimation = codeRate;
    
    % 误差分析
    if isfield(params.true, 'codeRate')
        trueRate = params.true.codeRate;
        error = abs(codeRate - trueRate);
        relativeError = error / trueRate * 100;
        
        info.trueValue = trueRate;
        info.absoluteError = error;
        info.relativeError = relativeError;
        
        fprintf('  真实码速率: %.3f MHz\n', trueRate/1e6);
        fprintf('  估计码速率: %.3f MHz\n', codeRate/1e6);
        fprintf('  绝对误差: %.3f kHz\n', error/1e3);
        fprintf('  相对误差: %.3f%%\n', relativeError);
    else
        fprintf('  估计码速率: %.3f MHz\n', codeRate/1e6);
    end
    
    fprintf('码速率估计完成\n');
end

function [codeRate, confidence] = estimateByAutocorrelation(signal, samplingRate)
    % 通过自相关函数估计码速率
    
    % 计算自相关函数
    autocorr = xcorr(real(signal), 'coeff');
    autocorr = autocorr(length(signal):end);  % 取正时延部分
    
    % 寻找周期性峰值
    [peaks, peakLocs] = findpeaks(autocorr, 'MinPeakHeight', 0.1, 'MinPeakDistance', 10);
    
    if length(peakLocs) < 2
        codeRate = samplingRate / 100;  % 默认值
        confidence = 0.3;
        return;
    end
    
    % 计算峰值间隔（码片周期）
    peakIntervals = diff(peakLocs);
    chipPeriod = median(peakIntervals);
    
    codeRate = samplingRate / chipPeriod;
    
    % 置信度基于峰值规律性
    intervalStd = std(peakIntervals);
    confidence = max(0, 1 - intervalStd / chipPeriod);
    
    % 可视化
    if false
        figure;
        plot((0:length(autocorr)-1)/samplingRate*1e6, autocorr);
        xlabel('时延 (μs)');
        ylabel('自相关系数');
        title('自相关函数分析');
        grid on;
        hold on;
        plot(peakLocs/samplingRate*1e6, peaks, 'ro');
        legend('自相关函数', '峰值');
    end
end

function [codeRate, confidence] = estimateByWavelet(signal, samplingRate)
    % 通过小波变换估计码速率
    
    % 使用小波变换检测码片边界
    signalReal = real(signal);
    
    % 选择合适的小波基
    [c, l] = wavedec(signalReal, 5, 'db4');
    
    % 分析细节系数
    detailCoeffs = detcoef(c, l, 3);  % 选择适当的分解层数
    
    % 检测细节系数中的跳变点（码片边界）
    threshold = 3 * median(abs(detailCoeffs)) / 0.6745;
    significantPoints = find(abs(detailCoeffs) > threshold);
    
    if length(significantPoints) < 2
        codeRate = samplingRate / 100;
        confidence = 0.3;
        return;
    end
    
    % 计算跳变点间隔
    intervals = diff(significantPoints);
    chipSamples = median(intervals);
    
    % 考虑小波分解的尺度
    scaleFactor = 2^3;  % 第3层分解
    actualChipSamples = chipSamples * scaleFactor;
    
    codeRate = samplingRate / actualChipSamples;
    
    % 置信度基于跳变点的规律性
    confidence = max(0, 1 - std(intervals) / median(intervals));
end

function [codeRate, confidence] = estimateByCyclicAutocorrelation(signal, samplingRate)
    % 通过循环自相关估计码速率
    
    signalReal = real(signal);
    N = length(signalReal);
    
    % 计算循环自相关
    maxLag = min(1000, round(N/10));
    cyclicCorr = zeros(1, maxLag);
    
    for tau = 1:maxLag
        if tau < N
            product = signalReal(1:end-tau) .* signalReal(tau+1:end);
            cyclicCorr(tau) = mean(product);
        end
    end
    
    % 寻找周期性
    [peaks, peakLocs] = findpeaks(cyclicCorr, 'MinPeakHeight', 0.05*max(cyclicCorr));
    
    if length(peakLocs) < 2
        codeRate = samplingRate / 100;
        confidence = 0.3;
        return;
    end
    
    % 估计码片周期
    fundamentalPeriod = findFundamentalPeriod(peakLocs);
    codeRate = samplingRate / fundamentalPeriod;
    
    % 置信度
    confidence = min(0.8, length(peaks) / 10);
end

function period = findFundamentalPeriod(peakLocs)
    % 寻找基本周期
    
    if length(peakLocs) < 2
        period = peakLocs(1);
        return;
    end
    
    % 计算所有峰值间隔
    intervals = diff(peakLocs);
    
    % 使用直方图找到最常见的间隔
    [counts, edges] = histcounts(intervals, 'BinMethod', 'integers');
    [~, maxIdx] = max(counts);
    period = round(mean([edges(maxIdx), edges(maxIdx+1)]));
end

function fusedRate = fuseCodeRateEstimations(rate1, rate2, rate3, conf1, conf2, conf3)
    % 融合码速率估计结果
    
    estimations = [rate1, rate2, rate3];
    confidences = [conf1, conf2, conf3];
    
    % 去除异常值
    validMask = confidences > 0.2;
    if sum(validMask) == 0
        fusedRate = mean(estimations);
        return;
    end
    
    validEstimations = estimations(validMask);
    validConfidences = confidences(validMask);
    
    % 加权平均
    totalConfidence = sum(validConfidences);
    weights = validConfidences / totalConfidence;
    
    fusedRate = sum(validEstimations .* weights);
end

扩频增益估计

function [spreadingGain, info] = estimateSpreadingGain(signal, codeRate, params)
    % 估计扩频增益
    
    fprintf('估计扩频增益...\n');
    
    info = struct();
    samplingRate = params.true.samplingRate;
    
    % 方法1: 自相关函数分析
    [gain1, confidence1] = estimateGainByAutocorrelation(signal, codeRate, samplingRate);
    
    % 方法2: 功率谱分析
    [gain2, confidence2] = estimateGainByPowerSpectrum(signal, codeRate, samplingRate);
    
    % 方法3: 统计特性分析
    [gain3, confidence3] = estimateGainByStatistics(signal, codeRate, samplingRate);
    
    % 融合估计结果
    spreadingGain = fuseGainEstimations(gain1, gain2, gain3, confidence1, confidence2, confidence3);
    
    info.methods.autocorr = gain1;
    info.methods.spectrum = gain2;
    info.methods.stats = gain3;
    info.confidence = [confidence1, confidence2, confidence3];
    info.finalEstimation = spreadingGain;
    
    % 误差分析
    if isfield(params.true, 'spreadingGain')
        trueGain = params.true.spreadingGain;
        error = abs(spreadingGain - trueGain);
        relativeError = error / trueGain * 100;
        
        info.trueValue = trueGain;
        info.absoluteError = error;
        info.relativeError = relativeError;
        
        fprintf('  真实扩频增益: %d\n', trueGain);
        fprintf('  估计扩频增益: %d\n', spreadingGain);
        fprintf('  绝对误差: %d\n', error);
        fprintf('  相对误差: %.3f%%\n', relativeError);
    else
        fprintf('  估计扩频增益: %d\n', spreadingGain);
    end
    
    fprintf('扩频增益估计完成\n');
end

function [gain, confidence] = estimateGainByAutocorrelation(signal, codeRate, samplingRate)
    % 通过自相关函数估计扩频增益
    
    % 计算自相关函数
    autocorr = xcorr(real(signal), 'coeff');
    autocorr = autocorr(length(signal):end);
    
    % 寻找周期性模式
    [peaks, peakLocs] = findpeaks(autocorr, 'MinPeakHeight', 0.05);
    
    if length(peakLocs) < 3
        gain = 63;  % 默认值
        confidence = 0.3;
        return;
    end
    
    % 分析峰值间隔模式
    intervals = diff(peakLocs);
    chipSamples = round(median(intervals));
    
    % 估计信息符号周期（寻找更大的周期）
    symbolIntervals = findSymbolPeriod(autocorr, chipSamples);
    
    if isempty(symbolIntervals)
        gain = 63;
        confidence = 0.3;
        return;
    end
    
    symbolSamples = median(symbolIntervals);
    gain = round(symbolSamples / chipSamples);
    
    % 置信度
    confidence = min(0.8, length(peaks) / 20);
end

function symbolIntervals = findSymbolPeriod(autocorr, chipSamples)
    % 寻找符号周期
    
    % 在自相关函数中寻找比码片周期更长的周期
    searchRange = round(chipSamples * 2):round(chipSamples * 100);
    symbolCorr = zeros(size(searchRange));
    
    for i = 1:length(searchRange)
        tau = searchRange(i);
        if tau < length(autocorr)
            symbolCorr(i) = autocorr(tau);
        end
    end
    
    % 寻找局部最大值
    [peaks, locs] = findpeaks(symbolCorr, 'MinPeakHeight', 0.1);
    
    if isempty(peaks)
        symbolIntervals = [];
        return;
    end
    
    symbolIntervals = searchRange(locs);
end

function [gain, confidence] = estimateGainByPowerSpectrum(signal, codeRate, samplingRate)
    % 通过功率谱估计扩频增益
    
    % 计算功率谱
    [psd, f] = pwelch(real(signal), hamming(512), 256, 512, samplingRate);
    
    % 分析频谱特征
    bandwidth = estimateBandwidth(psd, f);
    
    % 估计处理增益
    estimatedGain = round(bandwidth / codeRate);
    
    % 常见的扩频增益值
    commonGains = [7, 15, 31, 63, 127, 255, 511, 1023];
    
    % 选择最接近的常见值
    [~, closestIdx] = min(abs(estimatedGain - commonGains));
    gain = commonGains(closestIdx);
    
    % 置信度
    confidence = 0.6;
end

function bandwidth = estimateBandwidth(psd, f)
    % 估计信号带宽
    
    maxPSD = max(psd);
    threshold = maxPSD / 2;  % -3dB带宽
    
    aboveThreshold = psd > threshold;
    
    if ~any(aboveThreshold)
        bandwidth = f(end) - f(1);
        return;
    end
    
    startIdx = find(aboveThreshold, 1, 'first');
    endIdx = find(aboveThreshold, 1, 'last');
    
    bandwidth = f(endIdx) - f(startIdx);
end

function [gain, confidence] = estimateGainByStatistics(signal, codeRate, samplingRate)
    % 通过统计特性估计扩频增益
    
    signalReal = real(signal);
    
    % 计算信号的峰度（对扩频信号有特征性）
    kurt = kurtosis(signalReal);
    
    % 估计码片数（启发式方法）
    estimatedChips = round(1000 / (kurt - 1));  % 经验公式
    
    % 映射到常见的扩频增益
    commonGains = [7, 15, 31, 63, 127, 255, 511, 1023];
    [~, closestIdx] = min(abs(estimatedChips - commonGains));
    gain = commonGains(closestIdx);
    
    % 置信度基于峰度值
    confidence = min(0.7, abs(kurt - 1.5) / 5);
end

function fusedGain = fuseGainEstimations(gain1, gain2, gain3, conf1, conf2, conf3)
    % 融合扩频增益估计结果
    
    gains = [gain1, gain2, gain3];
    confidences = [conf1, conf2, conf3];
    
    % 选择最可信的估计
    [maxConf, bestIdx] = max(confidences);
    
    if maxConf > 0.5
        fusedGain = gains(bestIdx);
    else
        % 使用中位数
        fusedGain = median(gains);
    end
    
    % 取整到最接近的2^n-1值
    possibleGains = 2.^(3:10) - 1;
    [~, closestIdx] = min(abs(fusedGain - possibleGains));
    fusedGain = possibleGains(closestIdx);
end

参数验证和结果显示

function [validationResults, confidence] = validateEstimations(carrierFreq, codeRate, spreadingGain, params)
    % 参数验证和精度分析
    
    fprintf('参数验证和精度分析...\n');
    
    validationResults = struct();
    confidenceScores = zeros(1, 3);
    
    % 载频验证
    if isfield(params.true, 'carrierFreq')
        trueFreq = params.true.carrierFreq;
        freqError = abs(carrierFreq - trueFreq);
        freqAccuracy = max(0, 1 - freqError / trueFreq);
        confidenceScores(1) = freqAccuracy;
        
        validationResults.carrier.trueValue = trueFreq;
        validationResults.carrier.estimatedValue = carrierFreq;
        validationResults.carrier.error = freqError;
        validationResults.carrier.accuracy = freqAccuracy;
    else
        confidenceScores(1) = 0.7;  % 默认置信度
    end
    
    % 码速率验证
    if isfield(params.true, 'codeRate')
        trueRate = params.true.codeRate;
        rateError = abs(codeRate - trueRate);
        rateAccuracy = max(0, 1 - rateError / trueRate);
        confidenceScores(2) = rateAccuracy;
        
        validationResults.codeRate.trueValue = trueRate;
        validationResults.codeRate.estimatedValue = codeRate;
        validationResults.codeRate.error = rateError;
        validationResults.codeRate.accuracy = rateAccuracy;
    else
        confidenceScores(2) = 0.7;
    end
    
    % 扩频增益验证
    if isfield(params.true, 'spreadingGain')
        trueGain = params.true.spreadingGain;
        gainError = abs(spreadingGain - trueGain);
        gainAccuracy = max(0, 1 - gainError / trueGain);
        confidenceScores(3) = gainAccuracy;
        
        validationResults.spreadingGain.trueValue = trueGain;
        validationResults.spreadingGain.estimatedValue = spreadingGain;
        validationResults.spreadingGain.error = gainError;
        validationResults.spreadingGain.accuracy = gainAccuracy;
    else
        confidenceScores(3) = 0.7;
    end
    
    % 总体置信度
    confidence = mean(confidenceScores);
    
    fprintf('参数验证完成 - 总体置信度: %.2f\n', confidence);
end

function displayResults(carrierFreq, codeRate, spreadingGain, validationResults, confidence, params)
    % 显示估计结果
    
    fprintf('\n=== 直扩信号参数估计结果 ===\n');
    fprintf('载波频率: %.6f GHz\n', carrierFreq/1e9);
    fprintf('码速率: %.3f MHz\n', codeRate/1e6);
    fprintf('扩频增益: %d\n', spreadingGain);
    fprintf('总体估计置信度: %.2f\n', confidence);
    
    if isfield(validationResults, 'carrier')
        fprintf('\n--- 精度分析 ---\n');
        fprintf('载频误差: %.3f kHz\n', validationResults.carrier.error/1e3);
        fprintf('码速率误差: %.3f kHz\n', validationResults.codeRate.error/1e3);
        fprintf('扩频增益误差: %d\n', validationResults.spreadingGain.error);
    end
    
    % 可视化结果
    visualizeResults(carrierFreq, codeRate, spreadingGain, validationResults, confidence, params);
end

function visualizeResults(carrierFreq, codeRate, spreadingGain, validationResults, confidence, params)
    % 可视化估计结果
    
    figure('Position', [100, 100, 1400, 900]);
    
    % 1. 参数估计结果对比
    subplot(2,3,1);
    if isfield(validationResults, 'carrier')
        trueValues = [validationResults.carrier.trueValue/1e9, ...
                     validationResults.codeRate.trueValue/1e6, ...
                     validationResults.spreadingGain.trueValue];
        estimatedValues = [carrierFreq/1e9, codeRate/1e6, spreadingGain];
        
        bar([trueValues; estimatedValues]');
        set(gca, 'XTickLabel', {'载频 (GHz)', '码速率 (MHz)', '扩频增益'});
        legend('真实值', '估计值', 'Location', 'best');
        title('参数估计对比');
        grid on;
    else
        parameters = [carrierFreq/1e9, codeRate/1e6, spreadingGain];
        bar(parameters);
        set(gca, 'XTickLabel', {'载频 (GHz)', '码速率 (MHz)', '扩频增益'});
        title('参数估计结果');
        grid on;
    end
    
    % 2. 置信度显示
    subplot(2,3,2);
    confidenceValues = [confidence, 1-confidence];
    pie(confidenceValues, {'可信', '不确定'});
    title(sprintf('总体置信度: %.1f%%', confidence*100));
    
    % 3. 误差分析（如果有真实值）
    if isfield(validationResults, 'carrier')
        subplot(2,3,3);
        errors = [validationResults.carrier.error/validationResults.carrier.trueValue*100, ...
                 validationResults.codeRate.error/validationResults.codeRate.trueValue*100, ...
                 validationResults.spreadingGain.error/validationResults.spreadingGain.trueValue*100];
        bar(errors);
        ylabel('相对误差 (%)');
        set(gca, 'XTickLabel', {'载频', '码速率', '扩频增益'});
        title('参数估计相对误差');
        grid on;
    end
    
    % 4. 信号时域波形
    subplot(2,3,4);
    if isfield(params, 'testSignal')
        signal = params.testSignal;
        t = (0:length(signal)-1) / params.true.samplingRate * 1e6;
        plot(t, real(signal(1:min(10000, end))));
        xlabel('时间 (μs)');
        ylabel('幅度');
        title('接收信号时域波形');
        grid on;
    end
    
    % 5. 功率谱密度
    subplot(2,3,5);
    if isfield(params, 'testSignal')
        [psd, f] = pwelch(params.testSignal, hamming(1024), 512, 1024, params.true.samplingRate);
        plot(f/1e6, 10*log10(psd));
        xlabel('频率 (MHz)');
        ylabel('功率谱密度 (dB/Hz)');
        title('信号功率谱密度');
        grid on;
        hold on;
        plot(carrierFreq/1e6, max(10*log10(psd)), 'ro', 'MarkerSize', 8);
        legend('PSD', '估计载频');
    end
    
    % 6. 自相关函数
    subplot(2,3,6);
    if isfield(params, 'testSignal')
        autocorr = xcorr(real(params.testSignal), 'coeff');
        autocorr = autocorr(length(params.testSignal):end);
        tau = (0:length(autocorr)-1) / params.true.samplingRate * 1e6;
        plot(tau(1:min(1000, end)), autocorr(1:min(1000, end)));
        xlabel('时延 (μs)');
        ylabel('自相关系数');
        title('信号自相关函数');
        grid on;
    end
    
    sgtitle('直扩信号参数估计结果分析', 'FontSize', 14, 'FontWeight', 'bold');
end

function analyzePerformance(carrierFreq, codeRate, spreadingGain, params)
    % 性能分析
    
    fprintf('\n=== 系统性能分析 ===\n');
    
    % 计算处理增益
    processingGain = spreadingGain;
    theoreticalGainDB = 10 * log10(processingGain);
    
    fprintf('处理增益: %.1f dB\n', theoreticalGainDB);
    
    % 抗干扰能力分析
    jammingMargin = theoreticalGainDB - 10;  % 假设需要10dB信噪比
    fprintf('理论抗干扰容限: %.1f dB\n', jammingMargin);
    
    % 多址容量估计
    if spreadingGain > 0
        maxUsers = floor(spreadingGain / 10);  % 简化的容量估计
        fprintf('估计多址容量: %d 用户\n', maxUsers);
    end
    
    % 建议的应用场景
    fprintf('\n--- 建议应用场景 ---\n');
    if codeRate > 50e6
        fprintf('高速数据传输系统\n');
    elseif spreadingGain > 100
        fprintf('高抗干扰军事通信\n');
    else
        fprintf('常规扩频通信\n');
    end
    
    if carrierFreq > 1e9
        fprintf('微波频段通信系统\n');
    else
        fprintf('低频段通信系统\n');
    end
end

参考代码 直扩信号盲识别 www.youwenfan.com/contentcnk/65594.html

系统特点

1. 完整参数估计流程

• 载频估计: 功率谱分析、循环平稳特征、高阶统计量
• 码速率估计: 自相关分析、小波变换、循环自相关
• 扩频增益估计: 周期分析、频谱分析、统计特性

2. 多方法融合

• 多种估计算法并行运行
• 置信度加权融合
• 异常值检测和处理

3. 性能评估

• 精度分析和误差统计
• 置信度评估
• 应用场景建议

4. 可视化分析

• 参数对比显示
• 信号特征可视化
• 误差分布分析