[SciPy-User] scipy.signal.remez producing odd filter effects / not producing desired effect

Sun Jul 10 12:04:54 EDT 2016

On Sun, Jul 10, 2016 at 12:00 PM, Matti Viljamaa <mviljamaa at kapsi.fi> wrote:

> Or why doing
>
> fs = 44100
>
> # Design a low-pass filter using remez.
> cutoff = 2000.0
> transition_width = 200
> bands = np.array([0, cutoff - 0.5*transition_width,
>                   cutoff + 0.5*transition_width, fs/2.0]) / fs
> desired = [0, 1] # <------
> lpf = remez(513, bands, desired)
>
> T = 0.5
> nsamples = int(T*fs)
> t = np.linspace(0, T, nsamples, endpoint=False)
> freq = 440
> data = np.sin(2*np.pi*freq*t)
> data += np.sin(2*np.pi*(cutoff + 5*transition_width)*t)
> data /= 1.01*np.abs(data).max()
>
> filtered_data = lfilter(lpf, 1, data)
>
> w3, h3 = freqz(filtered_data)
>
>

Matti,

Did you see the email that I sent on July 6, and the plot that I attached
to it?  Here's what I said:

The arguments to 'freqz' are the filter coefficients, not the signal.  See
the example script in my previous email, where I have

    w, h = freqz(lpf)

Warren

> plt.figure(3)
> plt.plot(fs*w3/(2*np.pi), 20*np.log10(abs(h3)))
> plt.xlim(0, fs/2)
> plt.xlabel('Frequency (Hz)')
> plt.ylabel('Gain (dB)')
> plt.grid(True)
>
> plt. show()
>
> gives
>
>
> which doesn’t look like a highpass?
>
>
> On 07 Jul 2016, at 13:08, Matti Viljamaa <mviljamaa at kapsi.fi> wrote:
>
> Any idea why this plot shows positive gain?
> Like +50dB around 500Hz?
>
> -Matti
>
> On 06 Jul 2016, at 17:57, Warren Weckesser <warren.weckesser at gmail.com>
> wrote:
>
>
>
> On Wed, Jul 6, 2016 at 7:53 AM, Matti Viljamaa <mviljamaa at kapsi.fi> wrote:
>
>> Are these the proper frequency domain plots+
>>
>> w2, h2 = freqz(data)
>>
>>
>
> The arguments to 'freqz' are the filter coefficients, not the signal.  See
> the example script in my previous email, where I have
>
>     w, h = freqz(lpf)
>
> The attached plot is the frequency response that I get when I run that
> script.
>
> Warren
>
>
>
>
>> plt.figure(2)
>> plt.plot(fs*w2/(2*np.pi), 20*np.log10(abs(h2)))
>> plt.xlim(0, fs/2)
>> plt.xlabel('Frequency (Hz)')
>> plt.ylabel('Gain (dB)')
>> plt.grid(True)
>>
>> w3, h3 = freqz(filtered_data)
>>
>> plt.figure(3)
>> plt.plot(fs*w3/(2*np.pi), 20*np.log10(abs(h3)))
>> plt.xlim(0, fs/2)
>> plt.xlabel('Frequency (Hz)')
>> plt.ylabel('Gain (dB)')
>> plt.grid(True)
>>
>> plt. show()
>>
>> <Screen Shot 2016-07-06 at 14.53.08.png>
>>
>> On 04 Jul 2016, at 18:45, Warren Weckesser <warren.weckesser at gmail.com>
>> wrote:
>>
>> Matti,
>>
>> I don't have your input file, so I can't reproduce your result exactly.
>> For a basic demonstration, it is probably easier to create a signal in the
>> code instead of requiring an input wav file.  Here's an example that does
>> that.  The input signal (called 'data') is the sum of two sine waves, one
>> at 440 Hz and one at 3000 Hz.  The cutoff frequency for the low-pass filter
>> is 2000 Hz. (Actually the transition from the pass-band to the stop-band is
>> from 1900 Hz to 2100 Hz.)    You should be able to see in the plot and hear
>> in the wav files that the high frequency component has been filtered out.
>>
>> Warren
>>
>> -----
>>
>> from __future__ import division
>>
>> from scipy.signal import remez, freqz, lfilter
>> from scipy.io import wavfile
>> import numpy as np
>> import matplotlib.pyplot as plt
>>
>>
>> fs = 44100
>>
>> # Design a low-pass filter using remez.
>> cutoff = 2000.0
>> transition_width = 200
>> bands = np.array([0, cutoff - 0.5*transition_width,
>>                   cutoff + 0.5*transition_width, fs/2.0]) / fs
>> desired = [1, 0]
>> lpf = remez(513, bands, desired)
>>
>> # Plot the frequency response of the filter.
>> w, h = freqz(lpf)
>> plt.figure(1)
>> plt.plot(fs*w/(2*np.pi), 20*np.log10(abs(h)))
>> plt.xlim(0, fs/2)
>> plt.xlabel('Frequency (Hz)')
>> plt.ylabel('Gain (dB)')
>> plt.grid(True)
>>
>> # Create a test signal with two frequencies: one inside the pass-band,
>> # are one far outside the pass-band that should be filtered out.
>> T = 0.5
>> nsamples = int(T*fs)
>> t = np.linspace(0, T, nsamples, endpoint=False)
>> freq = 440
>> data = np.sin(2*np.pi*freq*t)
>> data += np.sin(2*np.pi*(cutoff + 5*transition_width)*t)
>> data /= 1.01*np.abs(data).max()
>>
>> # Filter the input using lfilter. (Alternatively, convolution could be
>> used.)
>> filtered_data = lfilter(lpf, 1, data)
>>
>> # Plot the input and output in the same figure.
>> plt.figure(2)
>> plt.plot(t, data)
>> plt.plot(t, filtered_data)
>>
>> plt.show()
>>
>> # Save the test signal and the filtered signal as wav files.
>> wavfile.write("data.wav", fs, data)
>> wavfile.write("filtered_data.wav", fs, filtered_data)
>>
>>
>> -----
>>
>>
>>
>>
>> On Mon, Jul 4, 2016 at 9:52 AM, Matti Viljamaa <mviljamaa at kapsi.fi>
>> wrote:
>>
>>>
>>> On 04 Jul 2016, at 14:00, Matti Viljamaa <mviljamaa at kapsi.fi> wrote:
>>>
>>> I’m trying to use the scipy.signal.remez filter design function. What
>>> I’m doing is the following (with the help of
>>> http://www.ee.iitm.ac.in/~nitin/teaching/ee5480/firdesign.html):
>>>
>>> import os
>>> import scipy.io.wavfile
>>> from scipy import signal
>>> from pylab import *
>>>
>>> os.chdir("/Users/mviljamaa/Music")
>>> sr, data = scipy.io.wavfile.read(open(“sound.wav", 'r'))
>>>
>>> lpf = signal.remez(21, [0, 0.05, 0.1], [1.0, .0])
>>>
>>> # Plot the magnitude response
>>> w, h = signal.freqz(lpf)
>>> plot(w/(2*pi), 20*log10(abs(h)))
>>> show()
>>>
>>> # Filtered data
>>> sout = signal.lfilter(lpf, 1, data)
>>>
>>> scipy.io.wavfile.write("equalized.wav", sr, sout)
>>>
>>> —
>>>
>>> Now judging by the magnitude response of the filter lpf, it should be a
>>> lowpass of some sort (I’m not sure how to interpret the cutoff frequency
>>> yet).
>>>
>>> The output I get is substantially gained (up to the point where it
>>> sounds distorted) and I cannot hear the lowpass filter.
>>>
>>> What am I doing wrong?
>>>
>>>
>>>
>>> Using a slightly modified code (from http://pastebin.com/LPEtXdzx):
>>>
>>> import os
>>> import scipy.io.wavfile
>>> from scipy import signal
>>> from pylab import *
>>> import numpy as np
>>>
>>> os.chdir("/Users/mviljamaa/Music")
>>> sr, data = scipy.io.wavfile.read(open(“sound.wav", 'r'))
>>>
>>> fs = 44100
>>>
>>> bands = array([0,3500,4000,5500,6000,fs/2.0]) / fs
>>> desired = [1, 0, 0]
>>> lpf = signal.remez(513, bands, desired)
>>>
>>> from scipy.signal import freqz
>>> w, h = signal.freqz(lpf)
>>> plot(w/(2*pi), 20*log10(abs(h)))
>>> show()
>>>
>>> sout = signal.lfilter(lpf, 1, data)
>>>
>>> sout /=  1.05 * max(abs(sout))
>>>
>>> scipy.io.wavfile.write("equalized.wav", sr, sout)
>>>
>>>
>>> —
>>>
>>> This one is able to retain the gain, but I still don’t hear any lowpass
>>> filtering.
>>>
>>> _______________________________________________
>>> SciPy-User mailing list
>>> SciPy-User at scipy.org
>>> https://mail.scipy.org/mailman/listinfo/scipy-user
>>>
>>>
>> _______________________________________________
>> SciPy-User mailing list
>> SciPy-User at scipy.org
>> https://mail.scipy.org/mailman/listinfo/scipy-user
>>
>>
>>
>> _______________________________________________
>> SciPy-User mailing list
>> SciPy-User at scipy.org
>> https://mail.scipy.org/mailman/listinfo/scipy-user
>>
>>
> <freqresp.png>_______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> https://mail.scipy.org/mailman/listinfo/scipy-user
>
>
>
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> https://mail.scipy.org/mailman/listinfo/scipy-user
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20160710/e64530be/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Screen Shot 2016-07-10 at 18.58.09.png
Type: image/png
Size: 81049 bytes
Desc: not available
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20160710/e64530be/attachment.png>