99%信号功率的带宽分析：两种方法的比较

在信号处理和电信中，通常需要估计信号的有效带宽。一种流行的方法是计算包含99％信号功率的频段。特别注意将分析限制在一定频率范围内的可能性，这在处理窄带信号时或在需要排除感兴趣频带外的噪声和干扰时尤其重要。此度量对于以下情况特别有用:

设计电信系统-确定所需的通信信道带宽
射频/微波系统工程师-滤波器和放大器的设计
信号处理专家-分析调制方法的频谱效率
声学和振动领域的研究人员-表征信号的频率组成

在本文中，我们比较了两种计算99％功率带的方法。:

直接FFT方法使用基本数字信号处理(DSP)功能
使用频谱分析仪的方法来自EngeeDSP库

using EngeeDSP
using LinearAlgebra
using Statistics

接下来，让我们继续实现本身，我们开发了3个主要功能:

calculate_99_percent_bandwidth -计算信号功率的99％的带宽。返回三个参数：总宽度、对称宽度和中心频率。添加了一个参数 band_range 以限制分析的范围。
bandwidth_99_direct -使用Welch方法直接计算FFT以估计功率谱。使用参数 band_range 以在计算之前对频率范围进行滤波。
bandwidth_99_spectrumanalyzer -使用内置的EngeeDSP频谱分析仪。添加参数 band_range （过滤结果）和 freq_span （设置分析器范围）。

功能 `calculate_99_percent_bandwidth`

函数计算包含总信号功率的99%的频率的带宽。接受频率数组 frequencies，对应的功率阵列 power 和一个可选参数 band_range 以限制所分析的范围。

运算的算法:

按范围过滤：如果指定 band_range，该函数仅将频率和功率值保留在指定的限值内。如果范围内没有数据，则返回错误。
功率归一化：计算总功率，对值进行归一化，使其总和为1。这允许您使用百分比。
基本带宽的计算:

-功率按降序排序
-有一组最小的组件，其累计总和达到99%
-带宽计算为该集合中最大和最小频率之间的差值
对称带宽的计算:

-频率按升序排序
-对应于累积功率的0.5％和99.5％的频率被发现
-这些频率之间的差异给出了对称的带宽
质心的计算：确定加权平均频率，其中相对功率值充当权重。

返回值：总带宽，对称带宽，频率质心。

function calculate_99_percent_bandwidth(frequencies, power; band_range=nothing)
    if band_range !== nothing
        min_freq, max_freq = band_range
        idx = findall(x -> min_freq <= x <= max_freq, frequencies)
        if isempty(idx)
            error("指定频率范围内没有数据[$min_freq,max max_freq]Hz")
        end
        frequencies = frequencies[idx]
        power = power[idx]
    end
    total_power = sum(power)
    if total_power == 0
        return 0.0, 0.0, 0.0
    end
    normalized_power = power / total_power
    sorted_indices = sortperm(normalized_power, rev=true)
    cumulative = cumsum(normalized_power[sorted_indices])
    idx_99 = findfirst(x -> x >= 0.99, cumulative)
    if idx_99 === nothing
        idx_99 = length(cumulative)
    end
    freq_indices_99 = sorted_indices[1:idx_99]
    freq_values = frequencies[freq_indices_99]
    bandwidth = maximum(freq_values) - minimum(freq_values)
    freq_order = sortperm(frequencies)
    cumulative_by_freq = cumsum(normalized_power[freq_order])
    idx_low = findfirst(x -> x >= 0.005, cumulative_by_freq)
    idx_high = findfirst(x -> x >= 0.995, cumulative_by_freq)
    if idx_low !== nothing && idx_high !== nothing
        freq_sorted = frequencies[freq_order]
        bandwidth_symmetric = freq_sorted[idx_high] - freq_sorted[idx_low]
    else
        bandwidth_symmetric = bandwidth
    end
    centroid = sum(frequencies .* normalized_power) / sum(normalized_power)
    return bandwidth, bandwidth_symmetric, centroid
end

calculate_99_percent_bandwidth (generic function with 1 method)

功能 `bandwidth_99_direct`

该函数使用Welch方法(periodogram averaging)计算99%的信号功率带。接受信号、采样频率和可选参数 band_range 来限制频率范围。

运算的算法:

加工参数:

-FFT大小：1024点或信号的长度（如果较小）
-将信号分成非重叠段进行平均
段处理:

-忽略小于64个计数的段
-汉纳窗口用于减少光谱泄漏
-正在执行窗口段FFT
-计算周期图，并将采样频率和窗口能量归一化
-从双侧频谱中提取单侧（正频）频谱，振幅加倍，以保持功率
平均：对所有段的周期图进行平均
频率轴创建：生成从0到奈奎斯特频率的频率
按范围过滤：如果指定 band_range，只剩下指定范围内的频率和PSD值
参数计算：函数被调用 calculate_99_percent_bandwidth 来计算频带宽度、对称宽度和质心

返回值：总带宽，对称带宽，质心，频率阵列，PSD阵列。

function bandwidth_99_direct(signal, sample_rate; band_range=nothing)
    nfft = min(1024, length(signal))
    segments = max(1, floor(Int, length(signal) / nfft))
    psd_total = zeros(nfft ÷ 2 + 1)
    for i in 1:segments
        start_idx = (i-1)*nfft + 1
        end_idx = min(i*nfft, length(signal))
        segment = signal[start_idx:end_idx]
        if length(segment) < 64
            continue
        end
        window = 0.5 .- 0.5*cos.(2π*(0:length(segment)-1)/(length(segment)-1))
        windowed = segment .* window
        fft_result = EngeeDSP.fft(windowed)
        psd = abs2.(fft_result) / (sample_rate * sum(window.^2))
        if length(psd) % 2 == 0
            idx = 1:length(psd)÷2 + 1
        else
            idx = 1:(length(psd)+1)÷2
        end
        psd_segment = psd[idx]
        if length(psd_segment) > 2
            psd_segment[2:end-1] .*= 2
        end
        if length(psd_segment) <= length(psd_total)
            psd_total[1:length(psd_segment)] .+= psd_segment
        end
    end
    psd_avg = psd_total ./ segments
    frequencies = range(0, sample_rate/2, length=length(psd_avg))
    if band_range !== nothing
        min_freq, max_freq = band_range
        idx = findall(x -> min_freq <= x <= max_freq, frequencies)
        if !isempty(idx)
            frequencies = frequencies[idx]
            psd_avg = psd_avg[idx]
        end
    end
    bandwidth, bandwidth_sym, centroid = calculate_99_percent_bandwidth(frequencies, psd_avg)
    return bandwidth, bandwidth_sym, centroid, frequencies, psd_avg
end

bandwidth_99_direct (generic function with 1 method)

功能 `bandwidth_99_spectrumanalyzer`

该函数使用内置的EngeeDSP频谱分析仪计算99%的信号功率带. 接受信号、采样率和两个可选参数: band_range 过滤结果及 freq_span 来调整分析仪的范围。

运算的算法:

配置分析器:

-如果设置 freq_span，所述分析器被配置为在指定范围内操作
-否则，全范围（0-fs/2）
被使用-使用以下方法：Welch方法，Hanna窗口，自动RBW设置
信号分析：将信号传输到分析仪进行处理
数据采集：提取频率和频谱值数组
按范围过滤：如果指定 band_range，只剩下指定边界内的频率和频谱值
参数计算：函数被调用 calculate_99_percent_bandwidth 来计算频带宽度、对称宽度和质心
资源释放：分析器使用后释放

返回值：总带宽，对称带宽，质心，频率阵列，频谱阵列。

function bandwidth_99_spectrumanalyzer(signal, sample_rate; band_range=nothing, freq_span=nothing)
    if freq_span !== nothing
        span_start, span_end = freq_span
        scope = EngeeDSP.spectrumAnalyzer(
            SampleRate = sample_rate,
            SpectrumType = "power",
            SpectrumUnits = "Watts",
            Method = "welch",
            Window = "Hann",
            RBWSource = "auto",
            FrequencySpan = "span-and-center-frequency",
            Span = span_end - span_start,
            CenterFrequency = (span_start + span_end) / 2,
            PlotAsTwoSidedSpectrum = false,
            AveragingMethod = "exponential",
            ForgettingFactor = 1.0,
            VBWSource = "auto",
            OverlapPercent = 0
        )
    else
        scope = EngeeDSP.spectrumAnalyzer(
            SampleRate = sample_rate,
            SpectrumType = "power",
            SpectrumUnits = "Watts",
            Method = "welch",
            Window = "Hann",
            RBWSource = "auto",
            FrequencySpan = "full",
            PlotAsTwoSidedSpectrum = false,
            AveragingMethod = "exponential",
            ForgettingFactor = 1.0,
            VBWSource = "auto",
            OverlapPercent = 0
        )
    end
    scope(signal)
    spectrum_data = getSpectrumData(scope)
    if haskey(spectrum_data, "spectrum") && haskey(spectrum_data, "frequencies")
        spectrum = spectrum_data["spectrum"]
        frequencies = spectrum_data["frequencies"]
        if ndims(spectrum) > 1
            spectrum = vec(mean(spectrum, dims=2))
        end
        if ndims(frequencies) > 1
            frequencies = vec(frequencies)
        end
        if band_range !== nothing
            min_freq, max_freq = band_range
            idx = findall(x -> min_freq <= x <= max_freq, frequencies)
            if !isempty(idx)
                frequencies = frequencies[idx]
                spectrum = spectrum[idx]
            end
        end
        bandwidth, bandwidth_sym, centroid = calculate_99_percent_bandwidth(frequencies, spectrum)
        release!(scope)
        return bandwidth, bandwidth_sym, centroid, frequencies, spectrum
    else
        release!(scope)
        error("无法获取频谱数据")
    end
end

bandwidth_99_spectrumanalyzer (generic function with 1 method)

开发算法的测试

本测试比较了两种计算信号99%带宽的方法，例如一个测试信号由三个正弦分量（1000hz、2500hz、4000hz）组成，加上白噪声。

测试顺序执行四个分析选项:

基本版本是FFT的直接计算，没有频率限制
直接带限方法-使用参数 band_range 用于过滤频率范围
有限频谱分析仪方法-使用内置的结果过滤工具
带范围调整的频谱分析仪方法-使用参数 freq_span 对于分析器本身的配置

每个选项返回三个关键参数（带宽，对称宽度，中心频率），它允许您直观地比较不同方法和参数对最终结果的影响。

function generate_test_signal(fs, duration)
    t = range(0, duration, step=1/fs)
    f1 = 1000
    f2 = 2500
    f3 = 4000
    signal = 1.0 * EngeeDSP.sin.(2π * f1 * t) +
             0.5 * EngeeDSP.sin.(2π * f2 * t) +
             0.2 * EngeeDSP.sin.(2π * f3 * t) +
             0.03 * randn(length(t))
    return signal
end

fs = 10000
duration = 0.1
signal = generate_test_signal(fs, duration)
println("1. 全频率范围:")
bw1, bw1_sym, centroid1, freqs1, power1 = bandwidth_99_direct(signal, fs)
println("   带宽:∞(round(bw1,digits=2))Hz")
println("   对称宽度：∞(round(bw1_sym,digits=2))Hz")
println("   中心频率：∞（圆形（质心1，数字=2））Hz")
println("2. 有限范围[500,3000]Hz:")
bw2, bw2_sym, centroid2, freqs2, power2 = bandwidth_99_direct(signal, fs; band_range=(500.0, 3000.0))
println("   带宽:∞(round(bw2,digits=2))Hz")
println("   对称宽度：∞(round(bw2_sym,digits=2))Hz")
println("   中心频率：∞（圆形（质心2，数字=2））Hz")
println("3. 光谱分析仪，量程[1000,3500]Hz:")
bw3, bw3_sym, centroid3, freqs3, power3 = bandwidth_99_spectrumanalyzer(signal, fs; band_range=(1000.0, 3500.0))
println("   带宽:∞(round(bw3,digits=2))Hz")
println("   对称宽度:∞(round(bw3_sym,digits=2))Hz")
println("   中心频率：∞（圆形（质心3，数字=2））Hz")
println("4. 具有freq_span设置[500,4500]Hz的SpectrumAnalyzer:")
bw4, bw4_sym, centroid4, freqs4, power4 = bandwidth_99_spectrumanalyzer(signal, fs; freq_span=(500.0, 4500.0))
println("   带宽:∞(round(bw4,digits=2))Hz")
println("   对称宽度:∞(round(bw4_sym,digits=2))Hz")
println("   中心频率：∞（圆形（质心4，数字=2））Hz")

1. Полный диапазон частот:
   Ширина полосы: 3020.0 Гц
   Симметричная ширина: 3020.0 Гц
   Центральная частота: 1385.74 Гц

2. Ограниченный диапазон [500, 3000] Гц:
   Ширина полосы: 1520.0 Гц
   Симметричная ширина: 1520.0 Гц
   Центральная частота: 1299.55 Гц

3. SpectrumAnalyzer с диапазоном [1000, 3500] Гц:
   Ширина полосы: 1604.82 Гц
   Симметричная ширина: 1582.03 Гц
   Центральная частота: 1542.26 Гц

4. SpectrumAnalyzer с настройкой freq_span [500, 4500] Гц:
   Ширина полосы: 3778.65 Гц
   Симметричная ширина: 3322.92 Гц
   Центральная частота: 1380.77 Гц

p1 = plot(freqs1, 10*log10.(power1 .+ 1e-10), 
         label="全系列", 
         linewidth=2,
         xlabel="频率(Hz)", 
         ylabel="功率(dB)",
         title="直接FFT-全范围",
         grid=true)
display(p1)
p2 = plot(freqs2, 10*log10.(power2 .+ 1e-10), 
         label="[500,3000]赫兹", 
         linewidth=2,
         xlabel="频率(Hz)", 
         ylabel="功率(dB)",
         title="直接FFT-有限范围",
         grid=true)
display(p2)
p3 = plot(freqs3, 10*log10.(power3 .+ 1e-10), 
         label="光谱分析仪[1000,3500]赫兹", 
         linewidth=2,
         xlabel="频率(Hz)", 
         ylabel="功率(dB)",
         title="带波段滤波的光谱分析仪",
         grid=true)
display(p3)
p4 = plot(freqs4, 10*log10.(power4 .+ 1e-10), 
         label="带freq_span[500,4500]Hz的光谱分析仪", 
         linewidth=2,
         xlabel="频率(Hz)", 
         ylabel="功率(dB)",
         title="带有freq_span设置的SpectrumAnalyzer",
         grid=true)
display(p4)

p5 = bar(["全", "[500,3000]", "SA [1000,3500]", "SA span [500,4500]"], 
         [bw1_sym, bw2_sym, bw3_sym, bw4_sym],
         xlabel="方法",
         ylabel="带宽(Hz)",
         title="对称带宽的比较99%",
         color=[:blue :red :green :purple],
         legend=false,
         xtickfontsize=8)

四种分析方法的比较显示频率范围限制对测量带宽参数的显着影响。全频率范围(测试1)显示最大带宽为3020Hz，复盖测试信号的所有三个正弦分量。

使用参数限制分析范围 band_range （测试2和3）导致测量带宽显着下降（分别高达1520Hz和1604Hz），因为排除了指定边界之外的频率分量。在这种情况下，中心频率取决于所选范围被移位。

与直接FFT方法相比，使用频谱分析仪（测试3和4）给出了更准确和详细的结果。参数的试验4特别显着。 freq_span，提供了信号最完整的特性（3778赫兹带宽），更接近三分量信号的理论期望值。频谱分析仪还演示了总带宽和对称带宽之间的差异，这表明对频域功率分布的更微妙考虑。

结论

所提出的方法使得能够有效地估计宽带和窄带信号的99%信号功率的带宽。参数 band_range 它提供了分析的灵活性，允许您只关注感兴趣的频谱部分，这在存在噪声和干扰的情况下尤为重要。

99%信号功率的带宽分析：两种方法的比较

功能 calculate_99_percent_bandwidth

功能 bandwidth_99_direct

功能 bandwidth_99_spectrumanalyzer

开发算法的测试

结论

功能 `calculate_99_percent_bandwidth`

功能 `bandwidth_99_direct`

功能 `bandwidth_99_spectrumanalyzer`