May 11, 2022 - Scalloping Loss in the FFT

using Gadfly
using ColorSchemes
using Statistics
using Printf
using FFTW
set_default_plot_size(20cm, 15cm)
style(background_color=colorant"white");

An astute reader may have noticed something odd when reading the first post on power measurement. The FFT size used was 10,000. Anyone who has learned about the FFT algorithm knows that algorithm is drastically faster when an FFT size that’s a power of 2 is used (eg 2^10 == 1024).

So why use an FFT size of 10,000? Well, there’s a source of loss that wasn’t accounted for in that post and under certain conditions this loss is 0 dB. This loss is called scalloping loss and scalloping loss is zero when the tone in the FFT (10 kHz in this case) falls on an FFT bin boundary.

When computing the frequency domain representation of a signal, the representation is broken up into bins. The sampling rate in the first post was 100 kHz and there were 10,000 bins meaning that each bin represents 10 Hz of bandwidth.

When the tone frequency falls on that 10 Hz bonudary (eg 40,500 Hz), there’s no scalloping loss. However, when the tone frequency doesn’t fall on a boundary (eg 40,509 Hz) the signal is affected.

Scalloping loss comes from the FFT window that’s being used. In this case no window is explicitly applied, which really means that a rectangualar window is being used (since all the samples before and after the signal are 0).

In this post, we’ll examine scalloping loss and how you can compensate for it. First, we’ll setup a signal that’s identical to the one used in the first post. Then, we’ll take 3 different FFTs of that signal and see the affect of scalloping. Next, we develop a method for compensating for scalloping loss. Finally, we plot scalloping loss directly and then the output of the signal with scalloping loss compensated for.

Step 1 - Signal Setup

# Signal setup
duration = 1                        # Signal durations [s]
fs = 100e3                          # Sampling rate, [Hz]
f = 10e3                            # Tone frequency [Hz]
t = 1/fs:1/fs:duration              # Signal time vector [s]
x = exp.(im * 2π * f * t)           # Signal with no noise
noise = randn(ComplexF64, size(x))  # Noise
xₙ = x .+ noise                     # xₙ, signal with noise 
;

Plot the first 100 samples of the signal

n = 100
p_1 = plot(Guide.title("x (complex valued signal)"), Guide.ylabel("Amplitude", orientation=:vertical), Guide.xlabel("Time [s]"))
p_2 = plot(Guide.title("xₙ (complex valued signal)"), Guide.ylabel("Amplitude",  orientation=:vertical), Guide.xlabel("Time [s]"))

push!(p_1, layer(x=t[1:n], y=real(x[1:n]), Geom.line, color=[colorant"blue"]))
push!(p_1, layer(x=t[1:n], y=imag(x[1:n]), Geom.line, color=[colorant"red"]))
push!(p_2, layer(x=t[1:n], y=real(xₙ[1:n]), Geom.line, color=[colorant"blue"]))
push!(p_2, layer(x=t[1:n], y=imag(xₙ[1:n]), Geom.line, color=[colorant"red"]))

stack = vstack(p_1, p_2)

<?xml version=”1.0” encoding=”UTF-8”?>

Step 2 - Three different FFTs of the same signal

In this section we see three different FFTs of the same signal. The only difference is the FFT Size (10000, 8192 and 4096).

The top plot shows the FFT output over the whole bandwidth and the bottom plot shows a zoom-in of the signal peak.

fft_size = 10000
X = (abs.(fft(xₙ[1:fft_size]))).^2 ./ fs ./ fft_size / (fft_size/length(xₙ))
peak_power_10000 = findmax(10 .* log10.(X))[1]
freqs = FFTW.fftfreq(fft_size, fs)

p_1 = plot()
push!(p_1, Guide.title("Frequency Response w/ FFT Size = $fft_size"), Guide.xlabel("Frequency [kHz]"), Guide.ylabel("Power Density [dB/Hz]"))
push!(p_1, Coord.cartesian(xmin=-fs/(2 * 1e3), xmax=fs/(2*1e3), ymin=-80, ymax=10))
push!(p_1, layer(x=freqs ./ 1e3, y=10 .* log10.(X), Geom.line, color=[colorant"blue"]))

p_2 = plot()
push!(p_2, Guide.title("Zoomed-In Frequency Response = $fft_size"), Guide.xlabel("Frequency [kHz]"), Guide.ylabel("Power Density [dB/Hz]"))
push!(p_2, Coord.cartesian(xmin=8, xmax=12, ymin=-5, ymax=5))
push!(p_2, layer(x=freqs ./ 1e3, y=10 .* log10.(X), Geom.line, color=[colorant"blue"]))
push!(p_2, layer(x=[-50, 50], y=[0, 0], Geom.line, color=[colorant"red"]))

stack = vstack(p_1, p_2)

<?xml version=”1.0” encoding=”UTF-8”?>

fft_size = 8192
X = (abs.(fft(xₙ[1:fft_size]))).^2 ./ fs ./ fft_size / (fft_size/length(xₙ))
peak_power_8192 = findmax(10 .* log10.(X))[1];
freqs = FFTW.fftfreq(fft_size, fs)

p_1 = plot()
push!(p_1, Guide.title("Frequency Response w/ FFT Size = 8192"), Guide.xlabel("Frequency [kHz]"), Guide.ylabel("Power Density [dB/Hz]"))
push!(p_1, Coord.cartesian(xmin=-fs/(2 * 1e3), xmax=fs/(2*1e3), ymin=-80, ymax=10))
push!(p_1, layer(x=freqs ./ 1e3, y=10 .* log10.(X), Geom.line, color=[colorant"blue"]))

p_2 = plot()
push!(p_2, Guide.title("Zoomed-In Frequency Response w/ FFT Size = 8192"), Guide.xlabel("Frequency [kHz]"), Guide.ylabel("Power Density [dB/Hz]"))
push!(p_2, Coord.cartesian(xmin=8, xmax=12, ymin=-10, ymax=10))
push!(p_2, layer(x=freqs ./ 1e3, y=10 .* log10.(X), Geom.line, color=[colorant"blue"]))
push!(p_2, layer(x=[-50, 50], y=[0, 0], Geom.line, color=[colorant"red"]))


stack = vstack(p_1, p_2)

<?xml version=”1.0” encoding=”UTF-8”?>

fft_size = 4096
X = (abs.(fft(xₙ[1:fft_size]))).^2 ./ fs ./ fft_size / (fft_size/length(xₙ))
freqs = FFTW.fftfreq(fft_size, fs)

p_1 = plot()
push!(p_1, Guide.title("Frequency Response w/ FFT Size = $fft_size"), Guide.xlabel("Frequency [kHz]"), Guide.ylabel("Power Density [dB/Hz]"))
push!(p_1, Coord.cartesian(xmin=-fs/(2 * 1e3), xmax=fs/(2*1e3), ymin=-80, ymax=10))
push!(p_1, layer(x=freqs ./ 1e3, y=10 .* log10.(X), Geom.line, color=[colorant"blue"]))

p_2 = plot()
push!(p_2, Guide.title("Zoomed-In Frequency Response w/ FFT Size = $fft_size"), Guide.xlabel("Frequency [kHz]"), Guide.ylabel("Power Density [dB/Hz]"))
push!(p_2, Coord.cartesian(xmin=8, xmax=12, ymin=-10, ymax=10))
push!(p_2, layer(x=freqs ./ 1e3, y=10 .* log10.(X), Geom.line, color=[colorant"blue"]))
push!(p_2, layer(x=[-50, 50], y=[0, 0], Geom.line, color=[colorant"red"]))

stack = vstack(p_1, p_2)

<?xml version=”1.0” encoding=”UTF-8”?>

These results show that as we change the FFT Size the signal power drops from 0 dB, to -0.5 dB to -2.5 dB.

Here we can directly see the affects of scalloping. We’re only changing the FFT Size and the power of the signal is decreasing, that doesn’t make sense!

@printf "Peak Power for FFT Size of 10000 = %2.3f dB\n" peak_power_10000
@printf "Peak Power for FFT Size of 8192 = %2.3f dB\n" peak_power_8192
@printf "Peak Power for FFT Size of 4096 = %2.3f dB\n" peak_power_4096
;

Peak Power for FFT Size of 10000 = 0.041 dB
Peak Power for FFT Size of 8192 = -0.552 dB
Peak Power for FFT Size of 4096 = -2.543 dB

Step 3 - Compensate for Scalloping Loss

Here we write a function to compensate for scalloping loss. The frequency domain representation of the rectangular window in a sinc() function. We apply that sinc function at the correct frequency to get the scalloping loss that we need to compensate for.

Note that scalloping loss that we compensate for matches the loss in step 2 for each FFT Size.

function sinc(x::Float64)
    sin(x)/x;
end

function compute_scalloping_loss(Fs::Float64, fft_size::Int64, f::Float64)
    res_bw = Fs/fft_size
    a = π/(res_bw) * mod((f - (res_bw/2)), res_bw) - π/2
    if a == 0
        # Avoid divide by 0 in sinc(x)
        b = 0  
    else
        b = 10 * log10(sinc(a) .^ 2)
    end
    b
end
;

sloss_10000 = compute_scalloping_loss(fs, 10000, f)
sloss_8192 = compute_scalloping_loss(fs, 8192, f)
sloss_4096 = compute_scalloping_loss(fs, 4096, f)

@printf "Scalloping loss for FFT Size = 10000 is %2.3f [dB]\n"  sloss_10000
@printf "Scalloping loss for FFT Size = 8192 is %2.3f [dB]\n"  sloss_8192
@printf "Scalloping loss for FFT Size = 4096 is %2.3f [dB]\n"  sloss_4096

Scalloping loss for FFT Size = 10000 is 0.000 [dB]
Scalloping loss for FFT Size = 8192 is -0.579 [dB]
Scalloping loss for FFT Size = 4096 is -2.420 [dB]

Step 4 - Plot Scalloping Loss and Compensating For It

In this last step a tone is swept from 9980 Hz to 10020 Hz in 0.08 Hz increments. The power of the tone is recorded for each increment. The power of the tone is then plotted as a function of frequency.

l = 500
X_comps = zeros(l)
scalloping_loss = zeros(l)
N = 10000
freq = zeros(l)

for (ii, kk) in enumerate(range(-20, 20, length=l))
    F = f + kk                              # Tone frequency 
    x = exp.(im * 2π * F .* t)              # Signal, a tone at frequency F
    X = fft(x[1:N])                         # FFT of the signal, FFT Size is 10k
    X_comp = abs.(X).^2 / N / fs / (N/fs)   # Signal, compensated for FFT Size, sampling rate and FFT bin size
    
    s_loss = compute_scalloping_loss(fs, N, F)  # Scalloping loss
    scalloping_loss[ii] = s_loss
    
    X_comps[ii] = 10*log10(findmax(abs.(X_comp))[1]) - s_loss  # X, compensated for all of the above and scalloping loss
    freq[ii] = F
end

Here we plot the scalloping loss as a function of frequency. Notice that when the FFT size divides the tone frequency evenly (eg 9990 Hz, 10000 Hz, etc) there is no scalloping loss.

plot(x=freq, 
     y=scalloping_loss, 
     Geom.line, 
     Guide.xlabel("Frequency [Hz]"), 
     Guide.ylabel("Scalloping Loss [dB]"),
     Guide.title("Scalloping Loss as a function of Frequency for FFT Size = $N"))

This is also a good place to discuss wherer the name “scalloped” comes from. It turns out that curtains and other textiles with this shape are referred to as “scalloped”.

A scalloped curtain:

In the last plot of this post, the power of the tone is shown as the tone is swept from 9980 Hz to 10020 Hz. All parameters are compensated for in this plot:

FFT Size (10k)
Sampling rate (100 kHz)
FFT Bin Size
Scalloping loss

With all of these items factored in, the power of the tone stays at a constant 0 dB/Hz for all frequencies examined.

plot(x=freq, 
     y=X_norms,
     Geom.line,  
     Coord.cartesian(ymin=-1, ymax=1),
     Guide.xlabel("Frequency [Hz]"), 
     Guide.ylabel("Peak Power [dB/Hz]"),
     Guide.title("Peak Power of Swept Tone after Compensating for Scalloping Loss"))

Dec 21, 2021 - Power Measurement

Power Measurement

This first post goes over measuring power of signals and computing signal to noise ratios. Differences between real-valued and complex-valued signals are discussed. Finally, signals with appreciable bandwidth are measured.

using Gadfly
using ColorSchemes
using Statistics
using Printf
using FFTW
set_default_plot_size(20cm, 15cm)
style(background_color=colorant"white");

Two functions to measure power, one for real valued vectors and a second for complex vectors.

function power(x::Vector{Float64})::Float64
    """
    power(x::Vector{Float64})::Float64

    Returns the average power (in dB) of the given sequence

    # Examples
    ```julia-repl
    julia> power([1.0, 2.0, 3.0, 4.0])
    8.75
    ```
    """
    10 .* log10.(mean(abs.(x) .^ 2))
end

function power(x::Vector{ComplexF64})::Float64
    """
    power(x::Vector{Float64})::Float64

    Returns the average power (in dB) of the given sequence

    # Examples
    ```julia-repl
    julia> power([1.0+1im, 2.0+2im, 3.0+3im, 4.0+4im])
    11.76
    ```
    """
    10 .* log10.(mean(abs.(x) .^ 2))
end
;

# Parameters for our signals
duration = 1             # Seconds
fs = 100e3               # Hz
f = 10e3                 # Hz
t = 1/fs:1/fs:duration   # Seconds
;

Measuring SNR For Real Valued Signals

In this example the intended SNR is set to 0 dB, a signal is constructed with unity gain and then a noise signal that will give the intended SNR is created. Then the power() function is used to ensure that the signal has the correct SNR.

To find the signal to noise ratio simply subtract the signal power (in dB) from the noise power (in dB). Remember that subtraction in log-space is division in linear-space. Note that the power of each signal is -3 dB and the SNR is 0 dB.

snr = 0                           # dB
snr_linear = 10 .^ (snr ./ 10)    # linear

x = cos.(2π * f .* t)
x_power = power(x)
x_power_linear = 10 .^ (x_power ./ 10)

noise_power_linear = x_power_linear ./ snr_linear
noise_power = 10*log10(noise_power_linear)

noise = sqrt(noise_power_linear) .* randn(Float64, size(x))

xₙ = x + noise
n_power = power(noise)
xₙ_power = power(xₙ)
snr_measured = x_power .- n_power
@printf "Avg Power of x = %2.2f dB\n" x_power
@printf "Avg Power of noise = %2.2f dB\n" n_power
@printf "Measured SNR = %2.2f dB\n" snr_measured
@printf "Expected SNR = %2.2f dB\n" snr

Avg Power of x = -3.01 dB
Avg Power of noise = -3.02 dB
Measured SNR = 0.01 dB
Expected SNR = 0.00 dB

Plot 100 samples of the tone, before noise is added and after.

n = 100  # number of sample to plot
p_1 = plot(Guide.title("x (real valued signal)"), Guide.ylabel("Amplitude", orientation=:vertical), Guide.xlabel("Time [s]"))
p_2 = plot(Guide.title("xₙ (real valued signal)"), Guide.ylabel("Amplitude", orientation=:vertical), Guide.xlabel("Time [s]"))

push!(p_1, layer(x=t[1:n], y=x[1:n], Geom.line, color=[colorant"blue"]))
push!(p_2, layer(x=t[1:n], y=xₙ[1:n], Geom.line, color=[colorant"blue"]))

stack = vstack(p_1, p_2)

<?xml version=”1.0” encoding=”UTF-8”?>

Next, plot the power spectral density. The PSD shows how the power of the signal is concentrated over the frequency band, 50 kHz in this case. Glancing at the PSD it might be tempting to think that the SNR is closer to 35 dB. However, the power of sine wave is concentrated entirely in one bin, where the noise is spread out over all the FFT bins. Since all the noise is spread out in frequency the PSD clearly shows the presence of the tone even at 0 dB SNR.

fft_size = 10000
xf = (abs.(fft(xₙ[1:fft_size]))).^2 / (fs / 2) / fft_size / (fft_size/length(xₙ))
nf = (abs.(fft(noise[1: fft_size]))).^2 / (fs / 2) / fft_size / (fft_size/length(xₙ))
noise_floor = 10 * log10.(mean(nf))

freqs = FFTW.fftfreq(length(xf), fs)
freq_response = plot()
push!(freq_response, Guide.title("PSD"), 
                     Guide.xlabel("Frequency [kHz]"), 
                     Guide.ylabel("Power Density [dB/Hz]", 
                     orientation=:vertical))
push!(freq_response, Coord.cartesian(xmin=-fs/(2 * 1e3), xmax=fs/(2*1e3), 
                                     ymin=-80, ymax=0))
push!(freq_response, layer(x=[-fs/(2 * 1e3), fs/(2*1e3)], 
                           y=[noise_floor, noise_floor], 
                           Geom.line, 
                           color=[colorant"red"]))
push!(freq_response, layer(x=freqs ./ 1e3, 
                           y=10 .* log10.(xf), 
                           Geom.line, 
                           color=[colorant"blue"]))

Compare the theoretical and measured power of the signal and noise.

Computing the measurements is simple. For the noise, take the mean() of the FFT output and for the tone, take the peak of FFT output.

Calculating the the noise power and signal power is a little more challenging. To calculate the expected noise measurement from the FFT, take the noise power and then calculate how it’s spread out over the usable bandwidth. There are two parts to this calculation. The first is two divide (or subtract in log space) by the usable bandwidth, since this is a real valued signal the usable bandwidth is fs / 2.

The second step is to divide (or subtract in log space) by the portion of the signal that the FFT uses. Longer signals have more power due to their length. Since the power() function uses the entire signal and the FFT has length 10000, that ratio needs to be compensated for.

The approach for calculating the theoretical peak power is similar for noise. The difference is that we must compensate for the bandwidth of the bin that tone occupies as well. The reason for this is that the tone is actually narrower than the bin so the power of the tone gets spread out over the bandwidth of the bin.

peak_power = findmax(10 .* log10.(xf))[1]

@printf "Measured Noise Power %2.1f dB/Hz\n" noise_floor
@printf "Theoretical Noise Power %2.1f dB/Hz\n" noise_power  - 10*log10((fs/2)) - 10*log10(fft_size/length(xₙ))
@printf "Measured Peak power %2.2f dB/Hz\n" peak_power
@printf "Theoretical Peak power %2.2f dB/Hz\n" (x_power + noise_power) - 10*log10((fs/2) / fft_size) - 10*log10(fft_size / length(xₙ)) 

Measured Noise Power -40.0 dB/Hz
Theoretical Noise Power -40.0 dB/Hz
Measured Peak power -3.15 dB/Hz
Theoretical Peak power -3.01 dB/Hz

Measuring SNR for Complex Valued Signals

For a complex signal, calculating the power is the same as a real signal, just count both the real and complex parts. Here, the power of the complex signal is 0 dB where the real signal was -3 dB. This is because the complex signal has real (cos) and imaginary (sin) components. Since the complex signal has twice as many components it has twice the power (ie 3 dB more).

snr = 0  # dB
snr_linear = 10 .^ (snr ./ 10)
x = im.* sin.(2π * f .* t) + cos.(2π * f .* t)
x_power = power(x)
x_power_linear = 10 .^ (x_power ./ 10)

noise_power_linear = x_power_linear ./ snr_linear
noise_power = 10*log10(noise_power_linear)

noise = sqrt(noise_power_linear) .* randn(ComplexF64, size(x))

snr_measured = x_power .- n_power

xₙ = x + noise
n_power = power(noise)
xₙ_power = power(xₙ)
snr_measured = x_power .- n_power
@printf "Avg Power of x = %2.2f dB\n" x_power
@printf "Avg Power of noise = %2.2f dB\n" n_power
@printf "SNR Measured = %2.2f dB\n" snr_measured
@printf "SNR Expected = %2.2f dB\n" snr

Avg Power of x = 0.00 dB
Avg Power of noise = 0.02 dB
SNR Measured = -0.02 dB
SNR Expected = 0.00 dB

Plot the time domain signal. Notice that both the real (inphase) and imaginary (quadrature) parts are present.

p_1 = plot(Guide.title("x (complex valued signal)"), Guide.ylabel("Amplitude", orientation=:vertical), Guide.xlabel("Time [s]"))
p_2 = plot(Guide.title("xₙ (complex valued signal)"), Guide.ylabel("Amplitude",  orientation=:vertical), Guide.xlabel("Time [s]"))

push!(p_1, layer(x=t[1:n], y=real(x[1:n]), Geom.line, color=[colorant"blue"]))
push!(p_1, layer(x=t[1:n], y=imag(x[1:n]), Geom.line, color=[colorant"red"]))
push!(p_2, layer(x=t[1:n], y=real(xₙ[1:n]), Geom.line, color=[colorant"blue"]))
push!(p_2, layer(x=t[1:n], y=imag(xₙ[1:n]), Geom.line, color=[colorant"red"]))

stack = vstack(p_1, p_2)

<?xml version=”1.0” encoding=”UTF-8”?>

Plot the PSD. Now there is only one tone instead of two. This is because a complex signal has twice as many samples as a real signal (I and Q vs just I). Since there are twice as many samples over the same time duration, complex signals are able to meet the Nyquist criterion at fs, not fs/2.

fft_size = 10000
xf = (abs.(fft(xₙ[1:fft_size]))).^2 ./ fs ./ fft_size / (fft_size/length(xₙ))
nf = (abs.(fft(noise[1:fft_size]))).^2 ./ fs ./ fft_size / (fft_size/length(xₙ))
freqs = FFTW.fftfreq(length(xf), fs)
noise_floor = 10 * log10.(mean(nf))

freq_response = plot()
push!(freq_response, Guide.title("Frequency Response"), Guide.xlabel("Frequency [kHz]"), Guide.ylabel("Power Density [dB/Hz]"))
push!(freq_response, Coord.cartesian(xmin=-fs/(2 * 1e3), xmax=fs/(2*1e3), ymin=-80, ymax=10))
push!(freq_response, layer(x=[-fs/(2 * 1e3), fs/(2*1e3)], y=[noise_floor, noise_floor], Geom.line, color=[colorant"red"]))
push!(freq_response, layer(x=freqs ./ 1e3, y=10 .* log10.(xf), Geom.line, color=[colorant"blue"]))

peak_power = findmax(10 .* log10.(xf))[1]

@printf "Measured Noise Power %2.1f dB/Hz\n" noise_floor 
@printf "Theoretical Noise Power %2.1f dB/Hz\n" n_power   - 10*log10(fs) - 10*log10(fft_size/length(xₙ))
@printf "Measured Peak power %2.1f dB/Hz\n" peak_power
@printf "Theoretical Peak power %2.1f dB/Hz" (x_power + noise_power) - 10*log10(fs / fft_size) - 10*log10(fft_size / length(xₙ))

Measured Noise Power -40.0 dB/Hz
Theoretical Noise Power -40.0 dB/Hz
Measured Peak power 0.0 dB/Hz
Theoretical Peak power 0.0 dB/Hz

Note here that the real and complex valued noise signals are off by 3 dB, -43 dB/Hz for the real signal and -40 dB/Hz for the complex. This is because in this simulation the noise power is determined the by the signal power and the SNR (0 dB). The complex signal has 3 dB more power so to maintain the 0 dB SNR the noise is 3 dB higher.

Also note that the complex tone has a peak power of 0 dB/Hz and the power of the real tone is -3 dB/Hz. The complex tone has twice the power from having the real and imaginary components.

Signals With Bandwidth

Measuring SNR

Measuring SNR is again the same for signals with bandwith.

function sinc(x)
   sin.(x) ./ (x); 
end

# Create sinc pulse and normalize to 0 dB
x₁ = sinc.(2π * f .* t .- 1000*π);
x₂ = x₁ ./ sqrt.(10 .^ (power(x₁)./10));
x = x₂
 
snr = 0
snr_linear = 10 .^ (snr ./ 10)
x_pwr = power(x)
x_pwr_linear = 10 .^ (x_pwr ./ 10)
noise = sqrt(x_pwr_linear ./ (snr_linear)) .* randn(Float64, size(x))

xₙ = x + noise
n_pwr = power(noise)
xₙ_pwr = power(xₙ)
snr = x_pwr .- n_pwr;

@printf "Avg Power of x = %2.2f dB\n" x_pwr
@printf "Avg Power of noise = %2.2f dB\n" n_pwr
@printf "SNR = %2.2f dB\n" snr

Avg Power of x = 0.00 dB
Avg Power of noise = -0.02 dB
SNR = 0.02 dB

Below, is the time domain plot of the sinc() function

p_1 = plot(Guide.title("x (complex valued signal)"), Guide.ylabel("Amplitude", orientation=:vertical), Guide.xlabel("Time [s]"))
p_2 = plot(Guide.title("xₙ (complex valued signal)"), Guide.ylabel("Amplitude", orientation=:vertical), Guide.xlabel("Time [s]"))

n = 10000
push!(p_1, layer(x=t[1:n], y=real(x[1:n]), Geom.line, color=[colorant"blue"]))
push!(p_1, layer(x=t[1:n], y=imag(x[1:n]), Geom.line, color=[colorant"red"]))
push!(p_2, layer(x=t[1:n], y=real(xₙ[1:n]), Geom.line, color=[colorant"blue"]))
push!(p_2, layer(x=t[1:n], y=imag(xₙ[1:n]), Geom.line, color=[colorant"red"]))

stack = vstack(p_1, p_2)

<?xml version=”1.0” encoding=”UTF-8”?>

The frequency domain plot is more interesting. Here are plots of the noise (wheat-colored) and signal (blue) separately along with the average noise power (red) and the average signal power (in-band, dark red).

The signal power is is normalized to the noise so the average power across the whole band is the same, making the SNR 0 dB. However, the signal power that “in-band,” -10 kHz to 10 kHz is 7 dB higher, making the SNR 7 dB for the “in-band” portion.

The SNR of our signal is 7 dB, meaning that the signal power is 7 dB greater than then noise power. Since the signal rolloff is so sharp almost all of it’s power is contained from -10 kHz to 10 kHz. Since the sample rate is 50 kHz the results are verified by calculating 10 * log10(10e3 / 50e3) = -6.9 dB, which matches the SNR.

# PSD
fft_size = length(x)  # Use and FFT size of the entire signal, ensuring all the signal power is represented
xf = (abs.(fft(x[1:fft_size])) ./ fft_size).^2 ./ (fs)
nf = (abs.(fft(noise[1:fft_size])) ./ fft_size).^2 ./ (fs)
xₙf = (abs.(fft(xₙ[1:fft_size])) ./ fft_size).^2 ./ (fs)
freqs = FFTW.fftfreq(length(xf), fs)
noise_floor = 10 * log10.(mean(nf))
ind1 = Int64(ceil(fft_size/2 - f/(fs) * fft_size/2))
ind2 = Int64(floor(fft_size/2 + f/(fs) * fft_size/2))

xf_inband = fftshift(fft(x[1:fft_size]))
xf_inband_pwr = 10*log10(mean((abs.(xf_inband[ind1: ind2]) / fft_size).^2)  / fs)
xf_allband_pwr = 10*log10(mean(abs.(xf_inband / fft_size).^2) / fs)

freq_response = plot(Coord.cartesian(xmin=-fs/(2 * 1e3), xmax=fs/(2*1e3), ymin=-120, ymax=-70),
                     Guide.title("Frequency Response"), 
                     Guide.xlabel("Frequency [kHz]"), 
                     Guide.ylabel("Power Density [dB/Hz]", orientation=:vertical))

push!(freq_response, layer(x=freqs ./ 1e3, 
                           y=10 .* log10.(xf), 
                           Geom.line, 
                           color=[colorant"blue"]))
push!(freq_response, layer(x=[-fs/(2*1e3), fs/(2*1e3)], 
                           y=[xf_inband_pwr, xf_inband_pwr], 
                           Geom.line,  color=[colorant"darkred"],
                           style(line_width=.5mm, line_style=[:dash])))
push!(freq_response, layer(x=[-fs/(2*1e3), fs/(2*1e3)], 
                           y=[noise_floor, noise_floor], 
                           Geom.line, 
                           color=[colorant"red"]))
push!(freq_response, layer(x=freqs ./ 1e3, y=10 .* log10.(xₙf), 
                           Geom.line, 
                           color=[colorant"wheat"]))

push!(freq_response, Guide.manual_color_key("Legend", ["Signal", "Signal Power", "Noise Power", "Noise"], 
                                                      ["blue", "darkred", "red", "wheat"]))

x_pwr = power(x)
n_pwr = power(noise)
snr = x_pwr - n_pwr
xₙ_pwr = power(xₙ)

@printf "Time Domain Measurements\n"
@printf "Avg noise power %2.1f dB\n" n_pwr
@printf "Avg signal (all band) power %2.1f dB\n" x_pwr
@printf "SNR (all band) = %2.1f dB\n\n" x_pwr - n_pwr

@printf "Frequency Domain Measurements\n"
@printf "Avg noise power (red line) %2.1f dB\n" noise_floor
@printf "Avg in band power (dark red, dashed line) %2.1f dB\n" xf_inband_pwr
@printf "Avg signal power (blue line) %2.1f dB\n" xf_allband_pwr
@printf "SNR (inband) = %2.1f dB\n" xf_inband_pwr - noise_floor
@printf "SNR (allband) = %2.1f dB\n" xf_allband_pwr - noise_floor

Time Domain Measurements
Avg noise power -0.0 dB
Avg signal (all band) power 0.0 dB
SNR (all band) = 0.0 dB

Frequency Domain Measurements
Avg noise power (red line) -100.0 dB
Avg in band power (dark red, dashed line) -93.0 dB
Avg signal power (blue line) -100.0 dB
SNR (inband) = 7.0 dB
SNR (allband) = 0.0 dB

Apr 30, 2021 - Direction Finding with MUSIC

using DSP
using FFTW
using Gadfly
using LinearAlgebra
using Statistics
using ColorSchemes
using Printf

function plot_signals(X::Matrix{ComplexF64}, n::Int64, fs::Float64)
    t = 1/fs:1/fs:n * fs
    n_elements = size(X)[1]
    p_i = plot(Guide.title("In Phase"), Guide.ylabel("Amplitude"), Guide.xlabel("Time [s]"))
    p_q = plot(Guide.title("Quadrature"), Guide.ylabel("Amplitude"), Guide.xlabel("Time [s]"))
    #colors = [colorant"blue", colorant"red", colorant"green", colorant"cyan", colorant"yellow"]
    colors = ColorSchemes.tab10.colors
    for i in 1:n_elements
        push!(p_i, layer(x=t[1:n], y=real(X[i, 1:n]), Geom.line, color=[colors[i]]))
        push!(p_q, layer(x=t[1:n], y=imag(X[i, 1:n]), Geom.line, color=[colors[i]]))
    end
    stack = vstack(p_i, p_q)
end;

Direction Finding with MUSIC

This post simulates direction finding with a 4 element uniform linear array (ULA) using the MUSIC algorithm.

There are three steps:

Generate a signal for each antenna with the appropriate phase delay for its angle of arrivale (AoA)
Apply the MUSIC algorithm
Evaluate the performance of the MUSIC algorithm

Step 1 - Generate Simulated Signals

For this simulation, there are two parts to generating the signals that are needed.

The first is generating four signals that are identical other than their phase. The phase for each signal gets delayed according to the intended angle of arrival (AoA) and the spacing of the array.

The second step is simply adding additive white Gaussian noise (AWGN).

power(X) = 10 .* log10.(mean(abs.(X) .^ 2, dims=2))
function add_noise!(noise_power::Float64, X::Matrix{ComplexF64})
    noise = sqrt((10 .^ (noise_power./10))) * randn(ComplexF64, size(X))
    X + noise
end;

function add_noise(power::Float64, X::Matrix{ComplexF64})
    noise = sqrt((10 .^ (noise_power./10))) * randn(ComplexF64, size(X))
    X + noise
end;

function gen_vectors(ϕ::Float64, fs::Float64, f::Float64, duration::Float64, n_elements::Int64)::Matrix{ComplexF64}
    t = 1/fs:1/fs:duration
    X = zeros(ComplexF64, n_elements, length(t))
    Φ = ϕ .* collect(0:1:n_elements-1)
    for i in 1:n_elements
        temp = im * sin.(2π * f .* t .+ Φ[i]) +  cos.(2π * f .* t .+ Φ[i]);
        X[i, :] = temp
    end
    X;
end;

# Generate the four signals as if they are coming from an emitter that's at 34° 

θ  = 34.0                  # true angle, degrees
ϕ  = π * cos(deg2rad(θ))   # phase difference, radians

fs = 30e3                  # sampling frquency, Hz
f  = 4e3                   # tone frequency, Hz
duration = 1.0             # signal duration, s

c = 299792458.0            # speed of light, m/s
λ = c / f;                 # wavelength, m
n_elements = 4;            # number of elements

X = gen_vectors(ϕ, fs, f, duration, n_elements);

plot_signals(X, 50, fs)

<?xml version=”1.0” encoding=”UTF-8”?>

noise_power = -3.0                                # dB, W
noise_power_linear = 10 .^ (noise_power ./ 10)    # W
snr = 3                                           # SNR, log-scale, unitless
snr_linear = 10 .^ (snr ./ 10)                    # SNR, linear-scale, unitless
X_pwr = power(X)                                  # Signal power, dB
X_pwr_linear = 10 .^ (X_pwr ./ 10)                # Signal power, linear
noise = sqrt.(X_pwr_linear ./ snr_linear) .* randn(ComplexF64, size(X))
Xn = X + noise
n_pwr = power(noise)

Xn_pwr = power(Xn)
snr = X_pwr .- n_pwr
@printf "Avg power of X = [%2.2f %2.2f %2.2f %2.2f] dB\n" X_pwr[1] X_pwr[2] X_pwr[3] X_pwr[4]
@printf "Avg power of n = [%2.2f %2.2f %2.2f %2.2f] dB\n" n_pwr[1] n_pwr[2] n_pwr[3] n_pwr[4]
@printf "Avg power of Xn = [%2.2f %2.2f %2.2f %2.2f] dB\n" Xn_pwr[1] Xn_pwr[2] Xn_pwr[3] Xn_pwr[4]
@printf "SNR = [%2.2f %2.2f %2.2f %2.2f] dB\n" snr[1] snr[2] snr[3] snr[4]
plot_signals(Xn, 50, fs)

Avg power of X = [0.00 0.00 0.00 0.00] dB
Avg power of n = [-3.00 -3.03 -3.03 -2.99] dB
Avg power of Xn = [1.75 1.77 1.74 1.77] dB
SNR = [3.00 3.03 3.03 2.99] dB

<?xml version=”1.0” encoding=”UTF-8”?>

Step 2 - Apply the MUSIC Algorithm

There are three main components to the MUSIC algorith: Compute the covariance matrix, obtain the eigenvalues and eigenvectors, apply the eigenvector with the largest eigenvalue to the array manifold.

The first component is computing the covariance matrix. In this example there are four signals of some length, N. When the covariance matrix is computed we are left with a 4x4 symmetric matrix that shows the covariance between each the signal that each antenna receives.

Next, obtain the eigenvalues and eigenvector pairs. The key to MUSIC is separating the signal subspace and the noise subspace by discarding the eigenvector with the smallest eigenvalue (as this eigenvector corresponds to the noise subspace).

Finally, take the remaining eigenvectors and create the MUSIC pseudospectrum by applying the array manifold to the eigenvectors. The array manifold is what the array expects to receive for a given angle of arrival under ideal circumstances.

function array_manifold_vector(array::Vector{Float64}, θ::Float64)::Vector{ComplexF64}
    exp.(-2π * im * cos(θ) .* array)
end;

function MUSIC(X::Matrix{ComplexF64}, n_targets::Int64, normalized_spacing::Float64)::Vector{Float64}
    n_elements = size(X)[1]
    element_locations = normalized_spacing .* 0.5 .* (n_elements .- 1 .- 2 .* collect(0:(n_elements-1)))
    pseudo_spectrum_length = 1024;
    vii = zeros(ComplexF64, n_elements, pseudo_spectrum_length);

    thetas = range(0, stop=π, length=pseudo_spectrum_length)
    for i in 1:pseudo_spectrum_length
        vii[:, i] = array_manifold_vector(element_locations, thetas[i])
    end
    
    # Get the covariance matrix, R
    R = Statistics.cor(X, X, dims=2);

    # Get the eigenvalues and vectors of R
    eigen_values = eigvals(R);
    eigen_vectors = eigvecs(R);
    
    # Sort the eigenvalues
    i = sortperm(abs.(eigen_values));
    eigen_vectors = eigen_vectors[:, i]

    # Use the largest n_elements - 1 eigenvectors to make U_n
    U_N = eigen_vectors[:, 1:n_elements-1];
    U_N_sq = U_N * U_N';

    # Generate the pseudospectrum
    pseudo_spectrum = zeros(Float64, 1, pseudo_spectrum_length)
    q = 0.0
    for i in 1:pseudo_spectrum_length
        q = vii'[[i], :] * U_N_sq * vii[:, [i]]
        pseudo_spectrum[i] = abs(1.0/real(q[1]))
    end

    # Normalize the pseudospectrum to 0dB
    pseudo_spectrum = 10 .* log10.(pseudo_spectrum/maximum(pseudo_spectrum));
    
    # Return the pseudospectrum
    vec(pseudo_spectrum)
end;

pseudo_spectrum = MUSIC(X, 1, 0.5)
(m, i) = findmax(pseudo_spectrum);
θ̂ = rad2deg(i/1024.0 * π);
error = abs(θ-θ̂)
@printf "θ = %2.2f°, θ̂ = %2.2f°, error = %2.4f°\n" θ θ̂ error

θ = 34.00°, θ̂ = 34.10°, error = 0.1016°

plot(x=range(0, stop=180, length=length(pseudo_spectrum)), 
     y=pseudo_spectrum, 
     Geom.line, 
     Coord.cartesian(xmin=0, xmax=180),
     Guide.title("MUSIC Pseudospectrum"),
     Guide.xlabel("Angle of Arrival (°)"),
     Guide.ylabel("MUSIC Response (dB)"))

Step 3 - MUSIC Evaluation

The cell below sweeps through every angle in 1 degree increments and records the AoA where the array and MUSIC perform the worst (\Delta max). This will likely be around 1 degree or 179 degrees due to the grating lobes from the ULA.

emitter_angles = 1.0:1.0:179.0
Δ = zeros(size(emitter_angles))
n_elements = 4
for (index, θ) in enumerate(emitter_angles)
    
    ϕ  = π * cos(deg2rad(θ))
    
    X = gen_vectors(ϕ, fs, f, duration, n_elements)
    noise = noise = sqrt((10 .^ (noise_power./10))) * randn(ComplexF64, size(X))

    pseudo_spectrum = MUSIC(X, 1, 0.5)

    (m, i) = findmax(pseudo_spectrum);
    θ̂ = rad2deg(i/1024.0 * π);

    Δ[index] = abs(θ - θ̂)
end

(Δmax, Δmax_ind) = findmax(Δ)
Δmax_angle = emitter_angles[Δmax_ind]

println("Δmax = $Δmax at θ = $Δmax_angle")

Δmax = 0.25 at θ = 11.0

plot(x=emitter_angles, y=Δ,  
     Coord.cartesian(xmin=0, xmax=180),
     Geom.line, 
     Guide.title("MUSIC Error as a function of AoA"),
     Guide.xlabel("Angle of Arrival (°)"),
     Guide.ylabel("Error (°)"))

Page: 1 of 1

Recent Posts

May 11, 2022 - Scalloping Loss in the FFT

Step 1 - Signal Setup

Step 2 - Three different FFTs of the same signal

Step 3 - Compensate for Scalloping Loss

Step 4 - Plot Scalloping Loss and Compensating For It

Dec 21, 2021 - Power Measurement

Power Measurement

Measuring SNR For Real Valued Signals

Measuring SNR for Complex Valued Signals

Signals With Bandwidth

Measuring SNR

Apr 30, 2021 - Direction Finding with MUSIC

Direction Finding with MUSIC

Step 1 - Generate Simulated Signals

Step 2 - Apply the MUSIC Algorithm

Step 3 - MUSIC Evaluation