1
$\begingroup$

I want to plot a specific frequency (e.g., 1 kHz or 100 Hz) of an audio file (sampling rate 44.1 kHz). I have done fft and spectrogram analysis on the file, but I am now interested in plotting only one frequency (as a sine wave) of the signal versus time. Here is also an example signal I created with only two frequencies and tried to pass only one of them, but I failed.

import numpy as np
from scipy.signal import butter, freqz
import matplotlib.pyplot as plt

def bandpass_filter(signal, lowcut, highcut, fs, order=5):
    nyquist = 0.5 * fs
    low = lowcut / nyquist
    high = highcut / nyquist
    b, a = butter(order, [low, high], btype='band')
    w, h = freqz(b, a, fs=fs)
    filtered_signal = np.convolve(signal, b / a, mode='same')
    return filtered_signal, w, h

# Example usage
fs = 44100  # Sample rate of the signal
f0 = 5   # Desired frequency to pass
lowcut = f0 - 1  # Lower cutoff frequency
highcut = f0 + 1 # Upper cutoff frequency

# Generate example signal
t = np.linspace(0, 1, num=fs, endpoint=False)
signal = np.sin(2 * np.pi * 5 * t) +  np.sin(2 * np.pi * 10 * t) # Signal with noise

t = t[:44100]
signal = signal[:44100]

# Apply bandpass filter
filtered_signal, w, h = bandpass_filter(signal, lowcut, highcut, fs)

# Plot original signal, filtered signal, and frequency response of the filter
#plt.figure(figsize=(12, 6))
fig, axes = plt.subplots(figsize=(12, 12), nrows=4, ncols=1, gridspec_kw={'hspace': 0.5, 'wspace': 2})
plt.subplot(4, 1, 1)
plt.plot(t, signal)
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.title('Original Signal')

plt.subplot(4, 1, 2)
plt.plot(t, filtered_signal)
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.title('Filtered Signal')

plt.subplot(4, 1, 3)
plt.plot(w, np.abs(h))
plt.xlabel('Frequency (Hz)')
plt.ylabel('Magnitude')
plt.title('Frequency Response of the Bandpass Filter')

plt.subplot(4, 1, 4)
plt.plot(t, signal-filtered_signal)
plt.xlabel('Time (s)')
plt.ylabel('Magnitude')
plt.title('Specific signal')
plt.show()
$\endgroup$
5
  • 1
    $\begingroup$ Do you know what band-pass filters are? $\endgroup$
    – Jdip
    Apr 11, 2023 at 19:42
  • $\begingroup$ Hi Jdip. Yes, I have created one. Here is a code that I have for the bandpass. The signal I generated only has two frequencies (5 Hz, and 10 Hz), and I set the high and low cuts to 5+1 and 5-1 to only allow the 5Hz pass. But the filtered signal amplitude is in the order of 1e-21. What I expected was to see a sin wave with a frequency of 5 Hz for the passed signal and 10 Hz for the blocked signal. I added the code to the question! $\endgroup$
    – pbena
    Apr 11, 2023 at 19:49
  • $\begingroup$ Chances are high that the filter frequency doesn't match the signal frequency due to some kind of unit confusion, although it doesn't jump out at me $\endgroup$
    – user253751
    Apr 11, 2023 at 21:03
  • $\begingroup$ There's a lot going on here. When I have time I'll try and provide a detailed answer. Start off with something easier. For now, try with sines at frequencies 2000 and 3000 (and cut-off 1600 and 2400). Your filter requirements are too stringent for such low frequencies and fs = 44100 $\endgroup$
    – Jdip
    Apr 11, 2023 at 21:14
  • $\begingroup$ I've written up something. Let me know once you've tried (do option 1 first) if you have more questions or if anything is un-clear. $\endgroup$
    – Jdip
    Apr 12, 2023 at 5:04

1 Answer 1

3
$\begingroup$

Note the very low pass-band amplitude of your filter. You're trying to implement a Butterworth band-pass filter with bandwidth $\texttt{BW} = 2 \texttt{Hz}$ at a sampling frequency $f_s = 44100 \texttt{Hz}$. That's a $\texttt{BW}/f_s$ ratio of $\approx 0.00005$, which is too stringent on your order-5 Butterworth filter.

To be a little more general than your use-case (for which you could simply reduce $f_s$ drastically and use a low-pass filter), I'll assume we have a real-world signal with large bandwidth and we're trying to isolate one component (which we'll keep at $5 \texttt{Hz}$).

There are a few options here, I'll go through 2 of them, I'll leave the coding to you for option 2.

Option 1: High order FIR

Easiest if you're not worried about processing time and want to keep your parameters as is, simply use a high-order FIR filter:

import numpy as np
from scipy.signal import freqz, firwin
import matplotlib.pyplot as plt

def bandpass_filter(signal, lowcut, highcut, fs, order=20000):
    nyquist = 0.5 * fs
    low = lowcut / nyquist
    high = highcut / nyquist
    b = firwin(order, [low, high], pass_zero=False)
    w, h = freqz(b, 1, worN=2*fs, fs=fs)
    filtered_signal = np.convolve(signal, b, mode='same')
    return filtered_signal, w, h

# Example usage
fs = 44100  # Sample rate of the signal
f0 = 5   # Desired frequency to pass
f1 = 10
lowcut = f0 - 1  # Lower cutoff frequency
highcut = f0 + 1 # Upper cutoff frequency

# Generate example signal
t = np.linspace(0, 1, num=fs, endpoint=False)
signal = np.sin(2 * np.pi * f0 * t) +  np.sin(2 * np.pi * f1 * t) 

# Apply bandpass filter
filtered_signal, w, h = bandpass_filter(signal, lowcut, highcut, fs)

# Plot frequency response of the filter
plt.figure(figsize=(12, 6))
plt.plot(w, np.abs(h))
plt.xlabel('Frequency (Hz)')
plt.ylabel('Magnitude')
plt.xlim(0,10)
plt.title('Frequency Response of the Bandpass Filter')
plt.show()

Option 2: Decimate then filter

One way to increase our $\texttt{BW}/f_s$ ratio is to decrease $f_s$ (and/or increase $\texttt{BW}$ of course), which in our case we can do since our signal's frequency component of interest sits at $5 \texttt{Hz}$ (remember the Nyquist rule $f_s \geq 2f_{\texttt{max}}$).
Choose $f_s = 100 \texttt{Hz}$ and $\texttt{BW} = 0.1\texttt{Hz}$ for example, increasing $\texttt{BW}/f_s$ to $0.001$.
We can then easily use either a lower-order FIR (see option 1 link for design) or IIR filter (I suggest you take a look at Robert Bristow-Johnson's Audio EQ Cookbook if you go the IIR route, look for the BPF section).

If you go with this option, do not forget to low-pass your signal before down-sampling, otherwise you'll get aliasing. Follow the link provided above for a decimation function that takes care of that for you.
(for your specific use-case, since your signal's highest component is at $10\texttt{Hz}$, you don't need to low-pass filter, you can simply set $f_s$ somewhere above $20\texttt{Hz}$).

Remarks

  • A second order IIR such as one designed using the cookbook can work with your $f_s$ as is, I’d reduce the passband bandwidth $\texttt{BW}$ though.

  • For IIR filtering, I would use filtfilt

$\endgroup$
6
  • $\begingroup$ Many thanks for your time. I tried the first option, and it works fine on low frequencies. But When I go to audio frequencies (44100 Hz), everything masses up, even if I set taps to 1001. $\endgroup$
    – pbena
    Apr 12, 2023 at 22:47
  • $\begingroup$ What do you mean if you go to “audio frequencies”? What frequency are you trying to isolate? Please edit your question by adding the problem you’re facing now. $\endgroup$
    – Jdip
    Apr 12, 2023 at 22:50
  • 2
    $\begingroup$ I vote for option 2! Nice work $\endgroup$ Apr 13, 2023 at 0:04
  • 1
    $\begingroup$ @DanBoschen That's where I would go to also, but for someone starting off in DSP, I thought I'd mention an easier path as well... $\endgroup$
    – Jdip
    Apr 13, 2023 at 0:17
  • $\begingroup$ Thanks for your good answer @Jdip. I have audio data, and I am trying to isolate only one frequency (e.g., 1 kHz). I tried your option 1, and it works beautifully for more straightforward single-frequency audio signals (100 Hz + 250 Hz, blocking 100 Hz). But when I try the real audio file and I change the ban-pass frequency, the passed signal shape does not change. $\endgroup$
    – pbena
    Apr 13, 2023 at 14:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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