Quantcast

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

classic Classic list List threaded Threaded
20 messages Options
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

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

Matti Viljamaa
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?

_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: scipy.signal.remez producing odd filter effects / not producing desired effect

Matti Viljamaa

On 04 Jul 2016, at 14:00, Matti Viljamaa <[hidden email]> 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
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: scipy.signal.remez producing odd filter effects / not producing desired effect

Warren Weckesser-2
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 <[hidden email]> wrote:

On 04 Jul 2016, at 14:00, Matti Viljamaa <[hidden email]> 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
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user



_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: scipy.signal.remez producing odd filter effects / not producing desired effect

Matti Viljamaa
Are these the proper frequency domain plots+

w2, h2 = freqz(data)

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()


On 04 Jul 2016, at 18:45, Warren Weckesser <[hidden email]> 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 <[hidden email]> wrote:

On 04 Jul 2016, at 14:00, Matti Viljamaa <[hidden email]> 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
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user


_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user


_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: scipy.signal.remez producing odd filter effects / not producing desired effect

Warren Weckesser-2


On Wed, Jul 6, 2016 at 7:53 AM, Matti Viljamaa <[hidden email]> 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()


On 04 Jul 2016, at 18:45, Warren Weckesser <[hidden email]> 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 <[hidden email]> wrote:

On 04 Jul 2016, at 14:00, Matti Viljamaa <[hidden email]> 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
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user


_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user


_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user



_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user

freqresp.png (64K) Download Attachment
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: scipy.signal.remez producing odd filter effects / not producing desired effect

Matti Viljamaa
Any idea why this plot shows positive gain?
Like +50dB around 500Hz?

-Matti

On 06 Jul 2016, at 17:57, Warren Weckesser <[hidden email]> wrote:



On Wed, Jul 6, 2016 at 7:53 AM, Matti Viljamaa <[hidden email]> 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 <[hidden email]> 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 <[hidden email]> wrote:

On 04 Jul 2016, at 14:00, Matti Viljamaa <[hidden email]> 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
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user


_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user


_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user


<freqresp.png>_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user


_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: scipy.signal.remez producing odd filter effects / not producing desired effect

Matti Viljamaa
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)

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 <[hidden email]> wrote:

Any idea why this plot shows positive gain?
Like +50dB around 500Hz?

-Matti

On 06 Jul 2016, at 17:57, Warren Weckesser <[hidden email]> wrote:



On Wed, Jul 6, 2016 at 7:53 AM, Matti Viljamaa <[hidden email]> 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 <[hidden email]> 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 <[hidden email]> wrote:

On 04 Jul 2016, at 14:00, Matti Viljamaa <[hidden email]> 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
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user


_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user


_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user


<freqresp.png>_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user



_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: scipy.signal.remez producing odd filter effects / not producing desired effect

Warren Weckesser-2


On Sun, Jul 10, 2016 at 12:00 PM, Matti Viljamaa <[hidden email]> 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 <[hidden email]> wrote:

Any idea why this plot shows positive gain?
Like +50dB around 500Hz?

-Matti

On 06 Jul 2016, at 17:57, Warren Weckesser <[hidden email]> wrote:



On Wed, Jul 6, 2016 at 7:53 AM, Matti Viljamaa <[hidden email]> 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 <[hidden email]> 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 <[hidden email]> wrote:

On 04 Jul 2016, at 14:00, Matti Viljamaa <[hidden email]> 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
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user


_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user


_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user


<freqresp.png>_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user



_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user



_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: scipy.signal.remez producing odd filter effects / not producing desired effect

Matti Viljamaa

On 10 Jul 2016, at 19:04, Warren Weckesser <[hidden email]> wrote:



On Sun, Jul 10, 2016 at 12:00 PM, Matti Viljamaa <[hidden email]> 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 <[hidden email]> wrote:

Any idea why this plot shows positive gain?
Like +50dB around 500Hz?

-Matti

On 06 Jul 2016, at 17:57, Warren Weckesser <[hidden email]> wrote:



On Wed, Jul 6, 2016 at 7:53 AM, Matti Viljamaa <[hidden email]> 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 <[hidden email]> 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 <[hidden email]> wrote:

On 04 Jul 2016, at 14:00, Matti Viljamaa <[hidden email]> 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
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user


_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user


_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user


<freqresp.png>_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user



_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user


_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user


_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: scipy.signal.remez producing odd filter effects / not producing desired effect

Warren Weckesser-2


On Sun, Jul 10, 2016 at 12:08 PM, Matti Viljamaa <[hidden email]> wrote:

On 10 Jul 2016, at 19:04, Warren Weckesser <[hidden email]> wrote:



On Sun, Jul 10, 2016 at 12:00 PM, Matti Viljamaa <[hidden email]> 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?



For that, you want to plot the spectral content of the original and filtered data.  You can use a Fourier transform and plot the magnitude of the Fourier coefficients against the frequencies, or you can use a function such as `scipy.signal.periodogram` or `scipy.signal.welch` which take care of the FFT details for you.

Here's a new version of my script.  It now creates a third figure containing plots of the periodogram of the original data and the filtered data.  I used `scipy.signal.periodogram`, but I recommend experimenting with `welch` also.

Warren

----------

from __future__ import division

from scipy.signal import remez, freqz, lfilter, periodogram
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 three frequencies: two 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) + 0.5*np.cos(2*np.pi*(2*freq)*t)
data += 0.75*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, 'b', label='original')
plt.plot(t, filtered_data, 'g', label='filtered')
plt.xlabel('Time (sec)')
plt.ylim(-1.1, 1.1)
plt.legend()

plt.figure(3)
freqs, data_spec = periodogram(data, fs=fs)
freqs, filtered_data_spec = periodogram(filtered_data, fs=fs)
plt.subplot(2, 1, 1)
plt.plot(freqs, data_spec, 'b', label='original')
plt.xlim(0, 3*cutoff)
plt.axvline(cutoff, color='k', alpha=0.15)
plt.legend()
plt.subplot(2, 1, 2)
plt.plot(freqs, filtered_data_spec, 'g', label='filtered')
plt.xlim(0, 3*cutoff)
plt.axvline(cutoff, color='k', alpha=0.15)
plt.legend()
plt.xlabel('Frequency (Hz)')

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)



----------

 
-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 <[hidden email]> wrote:

Any idea why this plot shows positive gain?
Like +50dB around 500Hz?

-Matti

On 06 Jul 2016, at 17:57, Warren Weckesser <[hidden email]> wrote:



On Wed, Jul 6, 2016 at 7:53 AM, Matti Viljamaa <[hidden email]> 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 <[hidden email]> 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 <[hidden email]> wrote:

On 04 Jul 2016, at 14:00, Matti Viljamaa <[hidden email]> 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
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user


_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user


_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user


<freqresp.png>_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user



_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user


_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user


_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user



_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user

figure_3.png (38K) Download Attachment
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: scipy.signal.remez producing odd filter effects / not producing desired effect

Warren Weckesser-2


On Sun, Jul 10, 2016 at 12:57 PM, Warren Weckesser <[hidden email]> wrote:


On Sun, Jul 10, 2016 at 12:08 PM, Matti Viljamaa <[hidden email]> wrote:

On 10 Jul 2016, at 19:04, Warren Weckesser <[hidden email]> wrote:



On Sun, Jul 10, 2016 at 12:00 PM, Matti Viljamaa <[hidden email]> 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?



For that, you want to plot the spectral content of the original and filtered data.  You can use a Fourier transform and plot the magnitude of the Fourier coefficients against the frequencies, or you can use a function such as `scipy.signal.periodogram` or `scipy.signal.welch` which take care of the FFT details for you.



Also, matplotlib has its own tools for this.  See, for example, the `psd` function demonstrated in this example: http://matplotlib.org/examples/pylab_examples/psd_demo.html

Warren


 
Here's a new version of my script.  It now creates a third figure containing plots of the periodogram of the original data and the filtered data.  I used `scipy.signal.periodogram`, but I recommend experimenting with `welch` also.

Warren

----------

from __future__ import division

from scipy.signal import remez, freqz, lfilter, periodogram
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 three frequencies: two 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) + 0.5*np.cos(2*np.pi*(2*freq)*t)
data += 0.75*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, 'b', label='original')
plt.plot(t, filtered_data, 'g', label='filtered')
plt.xlabel('Time (sec)')
plt.ylim(-1.1, 1.1)
plt.legend()

plt.figure(3)
freqs, data_spec = periodogram(data, fs=fs)
freqs, filtered_data_spec = periodogram(filtered_data, fs=fs)
plt.subplot(2, 1, 1)
plt.plot(freqs, data_spec, 'b', label='original')
plt.xlim(0, 3*cutoff)
plt.axvline(cutoff, color='k', alpha=0.15)
plt.legend()
plt.subplot(2, 1, 2)
plt.plot(freqs, filtered_data_spec, 'g', label='filtered')
plt.xlim(0, 3*cutoff)
plt.axvline(cutoff, color='k', alpha=0.15)
plt.legend()
plt.xlabel('Frequency (Hz)')

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)



----------

 
-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 <[hidden email]> wrote:

Any idea why this plot shows positive gain?
Like +50dB around 500Hz?

-Matti

On 06 Jul 2016, at 17:57, Warren Weckesser <[hidden email]> wrote:



On Wed, Jul 6, 2016 at 7:53 AM, Matti Viljamaa <[hidden email]> 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 <[hidden email]> 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 <[hidden email]> wrote:

On 04 Jul 2016, at 14:00, Matti Viljamaa <[hidden email]> 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
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user


_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user


_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user


<freqresp.png>_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user



_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user


_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user


_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user




_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: scipy.signal.remez producing odd filter effects / not producing desired effect

Matti Viljamaa
In reply to this post by Warren Weckesser-2

On 10 Jul 2016, at 19:57, Warren Weckesser <[hidden email]> wrote:



On Sun, Jul 10, 2016 at 12:08 PM, Matti Viljamaa <[hidden email]> wrote:

But what if I want to plot the effect of the lpf on the signal, rather than the filter magnitude response?



For that, you want to plot the spectral content of the original and filtered data.  You can use a Fourier transform and plot the magnitude of the Fourier coefficients against the frequencies, or you can use a function such as `scipy.signal.periodogram` or `scipy.signal.welch` which take care of the FFT details for you.

Here's a new version of my script.  It now creates a third figure containing plots of the periodogram of the original data and the filtered data.  I used `scipy.signal.periodogram`, but I recommend experimenting with `welch` also.

Warren

----------

from __future__ import division

from scipy.signal import remez, freqz, lfilter, periodogram
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 three frequencies: two 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) + 0.5*np.cos(2*np.pi*(2*freq)*t)
data += 0.75*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, 'b', label='original')
plt.plot(t, filtered_data, 'g', label='filtered')
plt.xlabel('Time (sec)')
plt.ylim(-1.1, 1.1)
plt.legend()

plt.figure(3)
freqs, data_spec = periodogram(data, fs=fs)
freqs, filtered_data_spec = periodogram(filtered_data, fs=fs)
plt.subplot(2, 1, 1)
plt.plot(freqs, data_spec, 'b', label='original')
plt.xlim(0, 3*cutoff)
plt.axvline(cutoff, color='k', alpha=0.15)
plt.legend()
plt.subplot(2, 1, 2)
plt.plot(freqs, filtered_data_spec, 'g', label='filtered')
plt.xlim(0, 3*cutoff)
plt.axvline(cutoff, color='k', alpha=0.15)
plt.legend()
plt.xlabel('Frequency (Hz)')

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)



—————

Thanks for this. However, I think there’s some other problem.

Modifying the above remez to the following:

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, cutoff + 0.5*transition_width + 1200.0, cutoff + 0.5*transition_width + 2400.0,fs/2.0]) / fs
desired = [0, 1, 0]
lpf = remez(513, bands, desired)

gives the following plots:


Why do these parameters lead to that +150dB peak, the peak in the “filtered" frequency plot and the strange full bandwidth burst in the last picture?

-Matti

_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: scipy.signal.remez producing odd filter effects / not producing desired effect

Warren Weckesser-2


On Mon, Jul 11, 2016 at 7:03 AM, Matti Viljamaa <[hidden email]> wrote:

On 10 Jul 2016, at 19:57, Warren Weckesser <[hidden email]> wrote:



On Sun, Jul 10, 2016 at 12:08 PM, Matti Viljamaa <[hidden email]> wrote:

But what if I want to plot the effect of the lpf on the signal, rather than the filter magnitude response?



For that, you want to plot the spectral content of the original and filtered data.  You can use a Fourier transform and plot the magnitude of the Fourier coefficients against the frequencies, or you can use a function such as `scipy.signal.periodogram` or `scipy.signal.welch` which take care of the FFT details for you.

Here's a new version of my script.  It now creates a third figure containing plots of the periodogram of the original data and the filtered data.  I used `scipy.signal.periodogram`, but I recommend experimenting with `welch` also.

Warren

----------

from __future__ import division

from scipy.signal import remez, freqz, lfilter, periodogram
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 three frequencies: two 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) + 0.5*np.cos(2*np.pi*(2*freq)*t)
data += 0.75*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, 'b', label='original')
plt.plot(t, filtered_data, 'g', label='filtered')
plt.xlabel('Time (sec)')
plt.ylim(-1.1, 1.1)
plt.legend()

plt.figure(3)
freqs, data_spec = periodogram(data, fs=fs)
freqs, filtered_data_spec = periodogram(filtered_data, fs=fs)
plt.subplot(2, 1, 1)
plt.plot(freqs, data_spec, 'b', label='original')
plt.xlim(0, 3*cutoff)
plt.axvline(cutoff, color='k', alpha=0.15)
plt.legend()
plt.subplot(2, 1, 2)
plt.plot(freqs, filtered_data_spec, 'g', label='filtered')
plt.xlim(0, 3*cutoff)
plt.axvline(cutoff, color='k', alpha=0.15)
plt.legend()
plt.xlabel('Frequency (Hz)')

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)



—————

Thanks for this. However, I think there’s some other problem.

Modifying the above remez to the following:

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, cutoff + 0.5*transition_width + 1200.0, cutoff + 0.5*transition_width + 2400.0,fs/2.0]) / fs
desired = [0, 1, 0]


You specified three intervals in `bands`: [0, 1900], [2100, 3300] and [4500, 22050] with gains 0, 1, and 0, respectively (i.e. you are designing a bandpass filter).  The remez algorithm designs a filter that is equiripple in the specified bands.  Outside those bands, however, the filter gain is unspecified; typically these intervals are described as the "don't care" intervals.  The remez algorithm makes no promises about the gain in these intervals, and it can--as you have found--result in unacceptable behavior.  The big peak in the filter gain occurs in the "don't care" interval [3300, 4500].

To fix it, you can try reducing the length of that "don't care" interval, or increasing the length of the filter (if that is an option).  It might require several iterations of band adjustments to achieve a desirable result.  You can also try the `weight` argument--perhaps give the stop bands a lower weight than the pass band.  Some experimentation and iteration will be involved.

Warren


lpf = remez(513, bands, desired)

gives the following plots:


Why do these parameters lead to that +150dB peak, the peak in the “filtered" frequency plot and the strange full bandwidth burst in the last picture?

-Matti

_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user



_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: scipy.signal.remez producing odd filter effects / not producing desired effect

Matti Viljamaa
On 11 Jul 2016, at 16:30, Warren Weckesser <[hidden email]> wrote:



On Mon, Jul 11, 2016 at 7:03 AM, Matti Viljamaa <[hidden email]> wrote:

On 10 Jul 2016, at 19:57, Warren Weckesser <[hidden email]> wrote:



On Sun, Jul 10, 2016 at 12:08 PM, Matti Viljamaa <[hidden email]> wrote:

But what if I want to plot the effect of the lpf on the signal, rather than the filter magnitude response?



For that, you want to plot the spectral content of the original and filtered data.  You can use a Fourier transform and plot the magnitude of the Fourier coefficients against the frequencies, or you can use a function such as `scipy.signal.periodogram` or `scipy.signal.welch` which take care of the FFT details for you.

Here's a new version of my script.  It now creates a third figure containing plots of the periodogram of the original data and the filtered data.  I used `scipy.signal.periodogram`, but I recommend experimenting with `welch` also.

Warren

----------

from __future__ import division

from scipy.signal import remez, freqz, lfilter, periodogram
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 three frequencies: two 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) + 0.5*np.cos(2*np.pi*(2*freq)*t)
data += 0.75*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, 'b', label='original')
plt.plot(t, filtered_data, 'g', label='filtered')
plt.xlabel('Time (sec)')
plt.ylim(-1.1, 1.1)
plt.legend()

plt.figure(3)
freqs, data_spec = periodogram(data, fs=fs)
freqs, filtered_data_spec = periodogram(filtered_data, fs=fs)
plt.subplot(2, 1, 1)
plt.plot(freqs, data_spec, 'b', label='original')
plt.xlim(0, 3*cutoff)
plt.axvline(cutoff, color='k', alpha=0.15)
plt.legend()
plt.subplot(2, 1, 2)
plt.plot(freqs, filtered_data_spec, 'g', label='filtered')
plt.xlim(0, 3*cutoff)
plt.axvline(cutoff, color='k', alpha=0.15)
plt.legend()
plt.xlabel('Frequency (Hz)')

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)



—————

Thanks for this. However, I think there’s some other problem.

Modifying the above remez to the following:

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, cutoff + 0.5*transition_width + 1200.0, cutoff + 0.5*transition_width + 2400.0,fs/2.0]) / fs
desired = [0, 1, 0]


You specified three intervals in `bands`: [0, 1900], [2100, 3300] and [4500, 22050] with gains 0, 1, and 0, respectively (i.e. you are designing a bandpass filter).  The remez algorithm designs a filter that is equiripple in the specified bands.  Outside those bands, however, the filter gain is unspecified; typically these intervals are described as the "don't care" intervals.  The remez algorithm makes no promises about the gain in these intervals, and it can--as you have found--result in unacceptable behavior.  The big peak in the filter gain occurs in the "don't care" interval [3300, 4500].

To fix it, you can try reducing the length of that "don't care" interval, or increasing the length of the filter (if that is an option).  It might require several iterations of band adjustments to achieve a desirable result.  You can also try the `weight` argument--perhaps give the stop bands a lower weight than the pass band.  Some experimentation and iteration will be involved.

Warren

So doing e.g.

bands = np.array([0, cutoff - 0.5*transition_width,
                  cutoff + 0.5*transition_width, cutoff + 0.5*transition_width + 2400.0, cutoff + 0.5*transition_width + 2500.0,fs/2.0]) / fs

i.e.

[     0.,   1900.]   [2100.,   4500.]   [4600.,  22050.]

the filter is shaped properly.

However,

do you know why

[     0.,   1900.]   [2100.,   4500.]   [4500.,  22050.]

fails? I.e. when the transition band is 0Hz?

-Matti


lpf = remez(513, bands, desired)

gives the following plots:

<Screen Shot 2016-07-11 at 14.01.40.png>

Why do these parameters lead to that +150dB peak, the peak in the “filtered" frequency plot and the strange full bandwidth burst in the last picture?

-Matti

_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user


_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user


_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: scipy.signal.remez producing odd filter effects / not producing desired effect

Matti Viljamaa
What about the following:

is this (the almost -15dB cut in the transition from <2000Hz to >2000Hz) a “feature" of the filter?


plotted with:

bands = np.array([0, cutoff - 0.5*transition_width,
                  cutoff + 0.5*transition_width, cutoff + 0.5*transition_width + 2400.0, cutoff + 0.5*transition_width + 2500.0,fs/2.0]) / fs
desired = [0.5, 1, 0.3]

-Matti

_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: scipy.signal.remez producing odd filter effects / not producing desired effect

Warren Weckesser-2
In reply to this post by Matti Viljamaa


On Mon, Jul 11, 2016 at 10:13 AM, Matti Viljamaa <[hidden email]> wrote:
On 11 Jul 2016, at 16:30, Warren Weckesser <[hidden email]> wrote:



On Mon, Jul 11, 2016 at 7:03 AM, Matti Viljamaa <[hidden email]> wrote:

On 10 Jul 2016, at 19:57, Warren Weckesser <[hidden email]> wrote:



On Sun, Jul 10, 2016 at 12:08 PM, Matti Viljamaa <[hidden email]> wrote:

But what if I want to plot the effect of the lpf on the signal, rather than the filter magnitude response?



For that, you want to plot the spectral content of the original and filtered data.  You can use a Fourier transform and plot the magnitude of the Fourier coefficients against the frequencies, or you can use a function such as `scipy.signal.periodogram` or `scipy.signal.welch` which take care of the FFT details for you.

Here's a new version of my script.  It now creates a third figure containing plots of the periodogram of the original data and the filtered data.  I used `scipy.signal.periodogram`, but I recommend experimenting with `welch` also.

Warren

----------

from __future__ import division

from scipy.signal import remez, freqz, lfilter, periodogram
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 three frequencies: two 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) + 0.5*np.cos(2*np.pi*(2*freq)*t)
data += 0.75*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, 'b', label='original')
plt.plot(t, filtered_data, 'g', label='filtered')
plt.xlabel('Time (sec)')
plt.ylim(-1.1, 1.1)
plt.legend()

plt.figure(3)
freqs, data_spec = periodogram(data, fs=fs)
freqs, filtered_data_spec = periodogram(filtered_data, fs=fs)
plt.subplot(2, 1, 1)
plt.plot(freqs, data_spec, 'b', label='original')
plt.xlim(0, 3*cutoff)
plt.axvline(cutoff, color='k', alpha=0.15)
plt.legend()
plt.subplot(2, 1, 2)
plt.plot(freqs, filtered_data_spec, 'g', label='filtered')
plt.xlim(0, 3*cutoff)
plt.axvline(cutoff, color='k', alpha=0.15)
plt.legend()
plt.xlabel('Frequency (Hz)')

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)



—————

Thanks for this. However, I think there’s some other problem.

Modifying the above remez to the following:

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, cutoff + 0.5*transition_width + 1200.0, cutoff + 0.5*transition_width + 2400.0,fs/2.0]) / fs
desired = [0, 1, 0]


You specified three intervals in `bands`: [0, 1900], [2100, 3300] and [4500, 22050] with gains 0, 1, and 0, respectively (i.e. you are designing a bandpass filter).  The remez algorithm designs a filter that is equiripple in the specified bands.  Outside those bands, however, the filter gain is unspecified; typically these intervals are described as the "don't care" intervals.  The remez algorithm makes no promises about the gain in these intervals, and it can--as you have found--result in unacceptable behavior.  The big peak in the filter gain occurs in the "don't care" interval [3300, 4500].

To fix it, you can try reducing the length of that "don't care" interval, or increasing the length of the filter (if that is an option).  It might require several iterations of band adjustments to achieve a desirable result.  You can also try the `weight` argument--perhaps give the stop bands a lower weight than the pass band.  Some experimentation and iteration will be involved.

Warren

So doing e.g.

bands = np.array([0, cutoff - 0.5*transition_width,
                  cutoff + 0.5*transition_width, cutoff + 0.5*transition_width + 2400.0, cutoff + 0.5*transition_width + 2500.0,fs/2.0]) / fs

i.e.

[     0.,   1900.]   [2100.,   4500.]   [4600.,  22050.]

the filter is shaped properly.

However,

do you know why

[     0.,   1900.]   [2100.,   4500.]   [4500.,  22050.]

fails? I.e. when the transition band is 0Hz?




You are asking for an ideal transition at 4500 Hz from perfect pass band to perfect stop band.  The remez algorithm doesn't converge.  You'd have to dig into the details of the algorithm for a good answer to "why", so I'll just say that the problem is too hard--possibly even impossible--for the remez algorithm.

Warren


 
-Matti


lpf = remez(513, bands, desired)

gives the following plots:

<Screen Shot 2016-07-11 at 14.01.40.png>

Why do these parameters lead to that +150dB peak, the peak in the “filtered" frequency plot and the strange full bandwidth burst in the last picture?

-Matti

_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user


_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user


_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user



_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: scipy.signal.remez producing odd filter effects / not producing desired effect

Warren Weckesser-2
In reply to this post by Matti Viljamaa


On Mon, Jul 11, 2016 at 10:23 AM, Matti Viljamaa <[hidden email]> wrote:
What about the following:

is this (the almost -15dB cut in the transition from <2000Hz to >2000Hz) a “feature" of the filter?


plotted with:

bands = np.array([0, cutoff - 0.5*transition_width,
                  cutoff + 0.5*transition_width, cutoff + 0.5*transition_width + 2400.0, cutoff + 0.5*transition_width + 2500.0,fs/2.0]) / fs
desired = [0.5, 1, 0.3]



That dip in the gain looks like it is in the "don't care" interval [cutoff - 0.5*transition_width, cutoff + 0.5*transition_width] (i.e. [1900, 2100]).   The remez algorithm doesn't specify the behavior in the "don't care" intervals.  If it turns out that you get something desirable, then sure, it's a feature, but if you tweak the parameters (i.e. the band edges or filter length), the behavior in the "don't care" intervals could be completely different.

Warren


-Matti

_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user



_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: scipy.signal.remez producing odd filter effects / not producing desired effect

Clancy Rowley
In reply to this post by Matti Viljamaa
Just a suggestion: for low-pass and band-pass filters, you might have better luck with Butterworth filters:

from scipy.signal import butter, freqz
fs = 44100
nyquist = 0.5/fs
passband = np.array([2000.,4500.]) * nyquist
order = 4
num, den = butter(order, passband, btype='bandpass')
w, h = freqz(num, den)

-clancy

> On Jul 11, 2016, at 10:23 AM, Matti Viljamaa <[hidden email]> wrote:
>
> What about the following:
>
> is this (the almost -15dB cut in the transition from <2000Hz to >2000Hz) a “feature" of the filter?
>
> <Screen Shot 2016-07-11 at 17.22.02.png>
>
> plotted with:
>
> bands = np.array([0, cutoff - 0.5*transition_width,
>                   cutoff + 0.5*transition_width, cutoff + 0.5*transition_width + 2400.0, cutoff + 0.5*transition_width + 2500.0,fs/2.0]) / fs
> desired = [0.5, 1, 0.3]
>
> -Matti
> _______________________________________________
> SciPy-User mailing list
> [hidden email]
> https://mail.scipy.org/mailman/listinfo/scipy-user

_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: scipy.signal.remez producing odd filter effects / not producing desired effect

Matti Viljamaa
In reply to this post by Warren Weckesser-2
Are there general ways for ensuring that the remez produces a filter that’s “correct”, i.e. no huge undesired boosts or other such “glitches”?

-Matti
_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user
Reply | Threaded
Open this post in threaded view
|  
Report Content as Inappropriate

Re: scipy.signal.remez producing odd filter effects / not producing desired effect

Warren Weckesser-2


On Tue, Jul 12, 2016 at 10:42 AM, Matti Viljamaa <[hidden email]> wrote:
Are there general ways for ensuring that the remez produces a filter that’s “correct”, i.e. no huge undesired boosts or other such “glitches”?



Not that I know of.  I always check the filter response with freqz.

Warren


 
-Matti
_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user


_______________________________________________
SciPy-User mailing list
[hidden email]
https://mail.scipy.org/mailman/listinfo/scipy-user
Loading...