0
$\begingroup$

Filtering noise leads to a change in the overall RMS of the noise as demonstrated by the following code-block. which uses scipy's firwin2 to calculate the taps. The first set of taps (taps1) generates a large notch from 1 to 40 kHz whereas the second set of taps (taps2) generates a small notch from 4 to 7 kHz.

import matplotlib.pyplot as plt
import numpy as np
from scipy import signal

def rms_db(x):
    rms = np.sqrt(np.mean(x ** 2))
    return 20 * np.log10(rms)

ntaps = 10001

freq1 = [0,  1e3, 2e3,  40e3, 41e3, 50e3]
gain1 = [1,  1,   1e-2, 1e-2, 1,    1   ]
taps1 = signal.firwin2(ntaps, freq1, gain1, fs=100e3)

freq2 = [0,  4e3, 5e3,  6e3,  7e3,  50e3]
gain2 = [1,  1,   1e-3, 1e-3, 1,    1   ]
taps2 = signal.firwin2(ntaps, freq2, gain2, fs=100e3)

# Generating uniform noise from -sqrt(3) to sqrt(3) results in RMS value of 1
noise = np.random.uniform(-np.sqrt(3), np.sqrt(3), size=100000)

# Trim off the first ntaps to eliminate boundary conditions
noise1 = signal.lfilter(taps1, 1, noise)[ntaps:]
noise2 = signal.lfilter(taps2, 1, noise)[ntaps:]
noise = noise[ntaps:]

print('Unfiltered noise RMS (dB re 1Vrms): {:.2f}'.format(rms_db(noise)))
print('Filtered noise1 RMS (dB re 1Vrms): {:.2f}'.format(rms_db(noise1)))
print('Filtered noise2 RMS (dB re 1Vrms): {:.2f}'.format(rms_db(noise2)))

As expected, filtering the noise leads to a reduction in RMS of approx. -6.65 dB for the large notch filter and approx. -0.19 dB for the small notch filter. Running the code gives us the following values.

Unfiltered noise RMS (dB re 1Vrms): 0.01
Filtered noise1 RMS (dB re 1Vrms): -6.65
Filtered noise2 RMS (dB re 1Vrms): -0.19

Since noise is inherently random, the actual reduction, in dB, of the noise will vary each time I run the code. I'd like to estimate the theoretical reduction in RMS that results from a particular when it's applied to a broadband signal (i.e., a signal containing equal power at all frequencies).

I thought I might be able to use freqz to estimate this:

w, h1 = signal.freqz(taps1, fs=fs)
h1 = np.abs(h1)

w, h2 = signal.freqz(taps2, fs=fs)
h2 = np.abs(h2)

plt.plot(w, 20 * np.log10(h1))
plt.plot(w, 20 * np.log10(h2))
plt.xlabel('Frequency (Hz)')
plt.ylabel('Gain (dB)')

average_db1 = 20 * np.log10(h1.mean())
average_db2 = 20 * np.log10(h2.mean())

print('Expected change in RMS for filter 1 (dB): {:.2f}'.format(average_db1))
print('Expected change in RMS for filter 2 (dB): {:.2f}'.format(average_db2))

However, this gives me an expected change in RMS that does not match the empirically-calculated change in RMS above.

Expected change in RMS for filter 1: -12.85
Expected change in RMS for filter 2: -0.35

Here is a plot of the frequency response of each filter.

Frequency response of the filters

Update

Per @Hilmar, I was doing the averaging incorrectly. I needed to average power, not amplitude. The correct calculation:

average_db1 = 20 * np.log10(np.sqrt(np.square(h1).mean()))
average_db2 = 20 * np.log10(np.sqrt(np.square(h2).mean()))
$\endgroup$
4
  • $\begingroup$ Re your final update, just simply do average_db1 = 10 * np.log10(np.square(h1).mean()) . No need to sqrt and then multiply by 20. (minor but it jumps out at me like saying 4 is the square of the square root of 4 :) ). ( @TimWescott you were searching for examples where you would use 10log) $\endgroup$ Jan 6 at 5:48
  • $\begingroup$ Thanks! A very valid point @DanBoschen. $\endgroup$
    – Brad
    Jan 6 at 18:28
  • $\begingroup$ I hope I didn't come across as being nit-picky! $\endgroup$ Jan 6 at 20:00
  • 1
    $\begingroup$ No. It's a good reminder about the difference between power and amplitude. Thanks! $\endgroup$
    – Brad
    Jan 7 at 14:45

1 Answer 1

0
$\begingroup$

You are off by a factor of two since you are averaging amplitudes. You should average the power instead (i.e. the square of the amplitude).

Off course, this only works if the input signal is white, i.e. equal energy at all freqeuncies.

$\endgroup$
1
  • $\begingroup$ Thanks! I updated the example above to reflect the correct calculation. $\endgroup$
    – Brad
    Jan 4 at 19:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.