[SciPy-User] scipy.signal.remez producing odd filter effects / not producing desired effect
Matti Viljamaa
mviljamaa at kapsi.fi
Sun Jul 10 12:08:49 EDT 2016
> On 10 Jul 2016, at 19:04, Warren Weckesser <warren.weckesser at gmail.com> wrote:
>
>
>
> On Sun, Jul 10, 2016 at 12:00 PM, Matti Viljamaa <mviljamaa at kapsi.fi <mailto: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
But what if I want to plot the effect of the lpf on the signal, rather than the filter magnitude response?
-Matti
>
> 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
>
> <Screen Shot 2016-07-10 at 18.58.09.png>
>
> which doesn’t look like a highpass?
>
>
>> On 07 Jul 2016, at 13:08, Matti Viljamaa <mviljamaa at kapsi.fi <mailto: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 <mailto:warren.weckesser at gmail.com>> wrote:
>>>
>>>
>>>
>>> On Wed, Jul 6, 2016 at 7:53 AM, Matti Viljamaa <mviljamaa at kapsi.fi <mailto: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 <mailto: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 <http://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 <mailto:mviljamaa at kapsi.fi>> wrote:
>>>>
>>>>> On 04 Jul 2016, at 14:00, Matti Viljamaa <mviljamaa at kapsi.fi <mailto: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 <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 <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 <mailto:SciPy-User at scipy.org>
>>>> https://mail.scipy.org/mailman/listinfo/scipy-user <https://mail.scipy.org/mailman/listinfo/scipy-user>
>>>>
>>>>
>>>> _______________________________________________
>>>> SciPy-User mailing list
>>>> SciPy-User at scipy.org <mailto:SciPy-User at scipy.org>
>>>> https://mail.scipy.org/mailman/listinfo/scipy-user <https://mail.scipy.org/mailman/listinfo/scipy-user>
>>>
>>>
>>> _______________________________________________
>>> SciPy-User mailing list
>>> SciPy-User at scipy.org <mailto:SciPy-User at scipy.org>
>>> https://mail.scipy.org/mailman/listinfo/scipy-user <https://mail.scipy.org/mailman/listinfo/scipy-user>
>>>
>>>
>>> <freqresp.png>_______________________________________________
>>> SciPy-User mailing list
>>> SciPy-User at scipy.org <mailto:SciPy-User at scipy.org>
>>> https://mail.scipy.org/mailman/listinfo/scipy-user <https://mail.scipy.org/mailman/listinfo/scipy-user>
>>
>
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org <mailto:SciPy-User at scipy.org>
> https://mail.scipy.org/mailman/listinfo/scipy-user <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/167e28c1/attachment.html>
More information about the SciPy-User
mailing list