python specific frequency remove(notch filter)?

Question

#complie by python3 only_test.py

import pyaudio
import numpy as np
import wave
import time
import math
#from pydub import AudioSegment
#from pydub.playback import play
#from scipy.signal import iirfilter
from scipy import signal

RATE = 48000
CHUNK = 4096
WIDTH = 2
volume = 0.0
duration = 1.0
#SHORT_NORMALIZE = (1.0/32768.0)
#INPUT_BLOCK_TIME = 1
#INPUT_BLOCK_PER_BLOCK = int(RATE*INPUT_BLOCK_TIME)

while True:

    #use a blackman window
    window = np.blackman(CHUNK)

    #load audio stream
    p = pyaudio.PyAudio()

    player = p.open(format=pyaudio.paInt16,
                    channels=1,
                    rate=RATE,
                    output=True,
                    frames_per_buffer=CHUNK)

    stream = p.open(format=pyaudio.paInt16,
                    channels=1,
                    rate=RATE,
                    input=True,
                    frames_per_buffer=CHUNK)

    #errorcount = 0

    for i in range(int(20*RATE/CHUNK)):

        sound = stream.read(CHUNK)

        #imp_ff = signal.filtfilt(b,a,sound)

        #playback microphone sound 
        #player.write(np.fromstring(sound,dtype=np.int16),CHUNK)

        #generate samples with return frequency to array
        #samples= (np.sin(2*np.pi*np.arange(RATE*duration)*freq/RATE)).astype(np.int16)

        #inverse frequency samples
        #inverse_samples = -samples

        #return frequency sound stream      
        #player.write(np.fromstring((volume*inverse_samples)\
        ,dtype=np.int16),CHUNK)

        #unpack the data and times by hamming window
        indata = np.array(wave.struct.unpack("%dh"%(len(sound)/WIDTH),\
                                             sound))*window

        #take the fft and square each value
        fftData = abs(np.fft.rfft(indata))*2

        #ifftData = abs(np.fft.irfft(indata))*2

        #find the maxium
        which = fftData[1:].argmax() + 1

        #use quadratic interpolation around the max
        if which != len(fftData)-1:
            y0,y1,y2 = np.log(fftData[which-1:which+2:])
            x1 = (y2-y0)*.5 / (2*y1-y2-y0)

            #find the frequency and output it
            freq = (which+x1)*RATE/CHUNK
            print("the freq is %d hz." % (freq))

        else:
            freq = which*RATE/CHUNK
            print("the freq is %d hz." % (freq))

        #playback the mic sound
        player.write(np.fromstring(sound,dtype=np.int16),CHUNK)

        if freq < 65:
           freq = 0

        #generate samples, note conversion to array
        #samples = 
        (np.sin(2*np.pi*np.arange(RATE*duration)*freq/RATE)).astype(np.int16)

        #invert phase of samples
        #result_samples = samples 

        #playback the invert_mic sound
        #player.write(np.fromstring(result_samples,dtype=np.int16),CHUNK)

    stream.stop_stream()
    stream.close()
    p.terminate()

We are currently processing microphones in real time. It is designed to obtain the frequency through it and to remove the sine wave sound through the notch filter (bandstop filter) for the output frequency. I do not know what code to write to do a notch filter (bandstop filter). Do you have any code or libraries to help?

Answer 1

While ARude's solution is decent, I made this one to be a bit more simple. Also, in the scipy docs the w0 term can be pretty confusing. It says it must be normalized, but it does it for you if you just enter the sampling rate.

I've wrapped such simplifications into the following code, and added a full working example with a fake signal where I remove the 60Hz hum:

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

# Create/view notch filter
samp_freq = 1000  # Sample frequency (Hz)
notch_freq = 60.0  # Frequency to be removed from signal (Hz)
quality_factor = 30.0  # Quality factor
b_notch, a_notch = signal.iirnotch(notch_freq, quality_factor, samp_freq)
freq, h = signal.freqz(b_notch, a_notch, fs = samp_freq)
plt.figure('filter')
plt.plot( freq, 20*np.log10(abs(h)))

# Create/view signal that is a mixture of two frequencies
f1 = 17
f2 = 60
t = np.linspace(0.0, 1, 1_000)
y_pure = np.sin(f1 * 2.0*np.pi*t) + np.sin(f2 * 2.0*np.pi*t) 
plt.figure('result')
plt.subplot(211)
plt.plot(t, y_pure, color = 'r')

# apply notch filter to signal
y_notched = signal.filtfilt(b_notch, a_notch, y_pure)

# plot notch-filtered version of signal
plt.subplot(212)
plt.plot(t, y_notched, color = 'r')

Answer 2

As you are already using scipy.signal, you can use the scipy.signal.iirnotch. Maybe you also want to read some background on IIR filters and for example the Quality Factor.

Use it as follows:

b, a = signal.iirnotch( w0, Q )

w0 is the normalised frequency.

Q is the quality factor that characterizes notch filter -3 dB bandwidth relative to its center frequency.

The function returns the numerator b and denominator a polynomials of the IIR filter.

Example:

fs = 200.0  # Sample frequency (Hz)
f0 = 60.0  # Frequency to be removed from signal (Hz)
Q = 30.0  # Quality factor
w0 = f0 / (fs / 2 )  # Normalized Frequency
b, a = signal.iirnotch( w0, Q )
# Look at frequency response
w, h = signal.freqz( b, a )
freq = w * fs / ( 2 * np.pi )
plt.plot( freq, 20*np.log10( abs( h ) ) )