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

Recent Posts

Dec 21, 2021 - Power Measurement

Power Measurement

Measuring SNR For Real Valued Signals

Measuring SNR for Complex Valued Signals

Signals With Bandwidth

Measuring SNR