Engee documentation
Notebook

Bandwidth analysis of 99% of the signal power: comparison of two methods

In signal processing and telecommunications, it is often necessary to estimate the effective bandwidth of a signal. One of the popular approaches is to calculate a band containing 99% of the signal power. Special attention is paid to the possibility of limiting the analysis to a certain frequency range, which is especially important when working with narrowband signals or when it is necessary to exclude noise and interference outside the band of interest. This metric is especially useful for:

  1. Designing telecommunication systems — to determine the required bandwidth of communication channels

  2. RF/microwave system engineers — in the design of filters and amplifiers

  3. Signal processing specialists — to analyze the spectral efficiency of modulation methods

  4. Researchers in the field of acoustics and vibrations — to characterize the frequency composition of signals

In this article, we compare two approaches to calculating 99% of the power band.:

  1. Direct FFT method using basic digital signal processing (DSP) functions

  2. Method using a spectrum analyzer from the EngeeDSP library

In [ ]: 
using EngeeDSP
using LinearAlgebra
using Statistics

Next, let's move on to the implementation itself, we have developed 3 main functions:

  1. calculate_99_percent_bandwidth - Calculates the bandwidth of 99% of the signal power. Returns three parameters: total width, symmetrical width, and center frequency. Added a parameter band_range to limit the range of analysis.

  2. bandwidth_99_direct - Direct calculation of the FFT using the Welch method to estimate the power spectrum. Uses the parameter band_range to filter the frequency range before calculating.

  3. bandwidth_99_spectrumanalyzer - Uses the built-in EngeeDSP spectrum analyzer. Added parameters band_range (filtering the results) and freq_span (setting the analyzer range).

Function calculate_99_percent_bandwidth

The function calculates the bandwidth of the frequency containing 99% of the total signal power. Accepts an array of frequencies frequencies, the corresponding power array power and an optional parameter band_range to limit the analyzed range.

The algorithm of operation:

  1. Filtering by range: If specified band_range, the function leaves only the frequencies and power values within the specified limits. If there is no data in the range, an error is returned.

  2. Power normalization: The total power is calculated, the values are normalized so that their sum is 1. This allows you to work with percentages.

  3. Calculation of the basic bandwidth:

    • Power is sorted in descending order
    • There is a minimum set of components, the cumulative sum of which reaches 99%
  • The bandwidth is calculated as the difference between the maximum and minimum frequency in this set

  1. Calculation of the symmetrical bandwidth:

    • Frequencies are sorted in ascending order
    • The frequencies corresponding to 0.5% and 99.5% of the cumulative power are found
    • The difference between these frequencies gives a symmetrical bandwidth

  2. Calculation of the centroid: The weighted average frequency is determined, where the relative power values act as weights.

Return values: total bandwidth, symmetrical bandwidth, frequency centroid.

In [ ]: 
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("There is no data in the specified frequency range [$min_freq, $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
Out[0]:
calculate_99_percent_bandwidth (generic function with 1 method)

Function bandwidth_99_direct

The function calculates 99% of the signal power band using the Welch method (periodogram averaging). Accepts a signal, a sampling frequency, and an optional parameter band_range to limit the frequency range.

The algorithm of operation:

  1. Processing parameters:

    • FFT size: 1024 points or the length of the signal (if smaller)
  • The signal is divided into non-overlapping segments for averaging

  1. Segment processing:

    • Segments shorter than 64 counts are ignored
    • The Hanna window is used to reduce spectral leakage
    • The window segment FFT is being performed
    • A periodogram is calculated with normalization to the sampling frequency and the energy of the window
  • A one-sided (positive frequencies) spectrum is extracted from the two-sided spectrum with doubling of amplitudes to preserve power

  1. Averaging: The periodograms of all segments are averaged

  2. Frequency axis creation: Frequencies from 0 to Nyquist frequency are generated

  3. Filtering by range: If specified band_range, only the frequencies and PSD values in the specified range are left

  4. Parameter calculation: The function is called calculate_99_percent_bandwidth to calculate the band width, symmetrical width, and centroid

Return values: total bandwidth, symmetrical bandwidth, centroid, frequency array, PSD array.

In [ ]: 
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
Out[0]:
bandwidth_99_direct (generic function with 1 method)

Function bandwidth_99_spectrumanalyzer

The function calculates 99% of the signal power band using the built-in EngeeDSP spectrum analyzer. Accepts a signal, a sampling frequency, and two optional parameters: band_range to filter the results and freq_span to adjust the range of the analyzer.

The algorithm of operation:

  1. Configuring the analyzer:

    • If set freq_span, the analyzer is configured to operate in the specified range
  • Otherwise, the full range (0 - fs/2)
    is used - The following methods are used: Welch method, Hanna window, automatic RBW setting

  1. Signal analysis: The signal is transmitted to the analyzer for processing

  2. Data acquisition: Arrays of frequencies and spectrum values are extracted

  3. Filtering by range: If specified band_range, only the frequencies and spectrum values within the specified boundaries are left

  4. Parameter calculation: The function is called calculate_99_percent_bandwidth to calculate the band width, symmetric width, and centroid

  5. Resource release: The analyzer is released after use

Return values: total bandwidth, symmetrical bandwidth, centroid, frequency array, spectrum array.

In [ ]: 
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("Couldn't get spectrum data")
    end
end
Out[0]:
bandwidth_99_spectrumanalyzer (generic function with 1 method)

Testing of the developed algorithms

This test demonstrates a comparison of two approaches to calculating the 99% bandwidth of a signal using the example of one test signal consisting of three sinusoidal components (1000 Hz, 2500 Hz, 4000 Hz) with the addition of white noise.

The test sequentially performs four analysis options:

  1. The basic version is a direct calculation of the FFT without frequency restrictions
  2. Direct band-limited method - using the parameter band_range for filtering the frequency range
  3. Limited spectrum analyzer method - using the built-in results filtering tool
  4. Spectrum analyzer method with range adjustment - using the parameter freq_span for the configuration of the analyzer itself

Each option returns three key parameters (bandwidth, symmetrical width, center frequency), which allows you to visually compare the impact of different methods and parameters on the final results.

In [ ]: 
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. Full frequency range:")
bw1, bw1_sym, centroid1, freqs1, power1 = bandwidth_99_direct(signal, fs)
println("   Bandwidth: $(round(bw1, digits=2)) Hz")
println("   Symmetrical width: $(round(bw1_sym, digits=2)) Hz")
println("   Center frequency: $(round(centroid1, digits=2)) Hz")
println("2. Limited range [500, 3000] Hz:")
bw2, bw2_sym, centroid2, freqs2, power2 = bandwidth_99_direct(signal, fs; band_range=(500.0, 3000.0))
println("   Bandwidth: $(round(bw2, digits=2)) Hz")
println("   Symmetrical width: $(round(bw2_sym, digits=2)) Hz")
println("   Center frequency: $(round(centroid2, digits=2)) Hz")
println("3. SpectrumAnalyzer with range [1000, 3500] Hz:")
bw3, bw3_sym, centroid3, freqs3, power3 = bandwidth_99_spectrumanalyzer(signal, fs; band_range=(1000.0, 3500.0))
println("   Bandwidth: $(round(bw3, digits=2)) Hz")
println("   Symmetrical width: $(round(bw3_sym, digits=2)) Hz")
println("   Center frequency: $(round(centroid3, digits=2)) Hz")
println("4. SpectrumAnalyzer with freq_span setting [500, 4500] Hz:")
bw4, bw4_sym, centroid4, freqs4, power4 = bandwidth_99_spectrumanalyzer(signal, fs; freq_span=(500.0, 4500.0))
println("   Bandwidth: $(round(bw4, digits=2)) Hz")
println("   Symmetrical width: $(round(bw4_sym, digits=2)) Hz")
println("   Center frequency: $(round(centroid4, digits=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 Гц
In [ ]: 
p1 = plot(freqs1, 10*log10.(power1 .+ 1e-10), 
         label="Full range", 
         linewidth=2,
         xlabel="Frequency (Hz)", 
         ylabel="Power (dB)",
         title="Direct FFT - Full Range",
         grid=true)
display(p1)
p2 = plot(freqs2, 10*log10.(power2 .+ 1e-10), 
         label="[500, 3000] Hz", 
         linewidth=2,
         xlabel="Frequency (Hz)", 
         ylabel="Power (dB)",
         title="Direct FFT - Limited Range",
         grid=true)
display(p2)
p3 = plot(freqs3, 10*log10.(power3 .+ 1e-10), 
         label="SpectrumAnalyzer [1000, 3500] Hz", 
         linewidth=2,
         xlabel="Frequency (Hz)", 
         ylabel="Power (dB)",
         title="SpectrumAnalyzer with band_range filtering",
         grid=true)
display(p3)
p4 = plot(freqs4, 10*log10.(power4 .+ 1e-10), 
         label="SpectrumAnalyzer with freq_span [500, 4500] Hz", 
         linewidth=2,
         xlabel="Frequency (Hz)", 
         ylabel="Power (dB)",
         title="SpectrumAnalyzer with the freq_span setting",
         grid=true)
display(p4)
In [ ]: 
p5 = bar(["Full", "[500,3000]", "SA [1000,3500]", "SA span [500,4500]"], 
         [bw1_sym, bw2_sym, bw3_sym, bw4_sym],
         xlabel="Method",
         ylabel="Bandwidth (Hz)",
         title="Comparison of symmetrical bandwidth 99%",
         color=[:blue :red :green :purple],
         legend=false,
         xtickfontsize=8)
Out[0]:

A comparison of the four analysis methods shows a significant effect of frequency range limitation on the measured bandwidth parameters. The full frequency range (test 1) demonstrates a maximum bandwidth of 3020 Hz, covering all three sinusoidal components of the test signal.

Limiting the range of analysis using the parameter band_range (tests 2 and 3) leads to a significant decrease in the measured bandwidth (up to 1520 Hz and 1604 Hz, respectively), since frequency components outside the specified boundaries are excluded. In this case, the center frequency is shifted depending on the selected range.

Using a spectrum analyzer (tests 3 and 4) gave more accurate and detailed results compared to the direct FFT method. Test 4 with the parameter is especially significant. freq_span, which provided the most complete characteristic of the signal (3778 Hz bandwidth), which is closer to the theoretically expected value for a three-component signal. The spectrum analyzer also demonstrates the difference between total and symmetrical bandwidth, which indicates a more subtle consideration of power distribution in the frequency domain.

Conclusion

The presented methods make it possible to effectively estimate the bandwidth of 99% of the signal power for both broadband and narrowband signals. Parameter band_range It provides flexibility in analysis, allowing you to focus only on the part of the spectrum of interest, which is especially important in the presence of noise and interference.