We start by verifying the linearity of the convolution operator. Let us consider the example presented in the course, where the system impulse response is given by h[n]=αnu[n], u[n] being the step function.
Let us take a system impulse response with α=0.999, and as input signal one of the tones used in the Fourier example "three tones". We will use the command
y = np.convolve(x,h,'valid')
in order to perform the convolution between the input x and the impulse response h. Notice that such a command uses only the available data without performing zero padding, in contrast to what is done by
y = np.convolve(x,h).
When using only the available data in performing a convolution between a sequence x of length M and a sequence h of length L, the convolution expression using only the available data can be written as:
if M>L, y[n]=∑Ll=1h[l]x[n−l],n=L+1,…,M
if L>M, y[n]=∑Mm=1x[m]h[n−m],n=M+1,…,L
The resulting signal is of length:
max(M−max(0,L−1),0)Therefore we have
%pylab inline
import numpy as np
import math as m
from scipy import constants as c
alpha = 0.999
h = pow(alpha,np.arange(100)) #construct the impulse response of the system of length 100
y = np.linspace(0,999,1000)
x1 = np.sin((y*2*c.pi*40)/1000)
y1 = np.convolve(x1,h,'valid')
Populating the interactive namespace from numpy and matplotlib
You can visualize or play the input x1 and the output y1 to
see the effect of the system. Notice that the dimension of y1 is
of 901
samples.
We now consider another one of the three tones of the Fourier example "three tones", namely
x2 = np.sin((y*2*c.pi*80)/1000)
The output of the system for the input x2 is given by
y2 = np.convolve(x2,h,'valid')
Once again, plotting or playing the output y2 gives us an idea how the input has been modified.
Given that the system is linear, its output when the input is a×x1+b×x2 must be equal to a×y1+b×y2. Let's check this out by taking a=2 and b=−10
a =2
b = -10
xx = a*x1 + b*x2
yy = np.convolve(xx,h,'valid')
Now we compute the error between the output of the system when the input is the linear combination xx=a×x1+b×x2 and the linear combination of the outputs a×y1+b×y2
error = abs(yy - (a*y1 + b*y2))
maxError = max(error)
If now we printout maxError, by simply typing
maxError
2.1316282072803006e-14
we obtain a number in the order of 10−14, a negligible value which is purely due to numerical computation errors. Therefore (ax1+bx2)×h=ay1+by2
Let us now now check the time invariance of a convolution system. Consider the same system as above with input x1 and output y1. We now introduce a time delay in the input signal, say of 100 samples and compute the corresponding output
x1Delayed = x1[100:1000]
y1InputDelayed = np.convolve(x1Delayed,h,'valid')
Given the time invariance of the system, y1InputDelayed
must be
equal to y1
with a time delay of 100 samples, i.e.,
y1InputDelayed=y1[100:901]
. Indeed
y1Delayed = y1[100:901]
error2 = abs(y1InputDelayed -y1Delayed)
maxError2 = max(error2)
By simply typing
print maxError2
0.0
M = 10
lbd = float(M-1)/float(M)
h = (1-lbd)*pow(lbd,np.arange(100)) # constructs the impulse response of the system of length 100
Generate now a Gaussian white noise (sequence of independent random numbers that are generated according to the same Gaussian distribution)
sigma2 = 0.1 #Power of the noise
noise = sigma2*np.random.randn(2000) # Gaussian noise
that we add to the sequence of two pulses x1 and x2 that we have used in the examples "three tones" of the Fourier module, that is
x1 = np.sin((y*2*c.pi*40)/1000)
x2 = np.sin((y*2*c.pi*80)/1000)
x = np.append(x1,x2)
xNoisy = noise + x # Noisy version of x
You can play or plot x
and xNoisy
to hear or see the effect
of the noise
subplot(2,1,1)
plot(np.arange(x.size),x)
ylabel('x')
xlim([825,1175])
subplot(2,1,2)
plot(np.arange(xNoisy.size),xNoisy)
xlabel('Samples')
ylabel('xNoisy')
xlim([825,1175])
show()
Notice that we have plotted here the values of the sample 825 to the sample 1175, where the transition between the two tones occurs.
We now apply the leaky integrator to the noisy signal
yConv = np.convolve(xNoisy,h,'valid')
Once again you can play or plot the output of the signal y to hear and see the effect of the filter. When plotting it, be careful that, due to the convolution of finite length sequences, the output sequence y is of length 1901 and the transition between the two tones happens at the sample number 901 (and not number 1000 as of the input signal xNoisy
)
plot(np.arange(yConv.size),yConv)
xlabel('Samples')
ylabel('yConv')
xlim([750,1050])
(750, 1050)
As you can see, the output is "smoothed", that is, denoised. You can also remark the different effect on the two tones: The tone x2=sin(2π80n/1000), n=1,…,1000, has been attenuated with respect to the tone x1=sin(2π40n/1000), n=1,…,1000. The different attenuation of the two tones will appear clear as soon we will discuss the convolution theorem.
So let's now apply our denoising method to a human voice signal.
Load the file jingle.mat
into Python using the command
import scipy.io as sio
mat_contents = sio.loadmat('jingle.mat')
jingle = mat_contents['jingle']
print mat_contents
{'jingle': array([[ 0. , 0. , 0. , ..., -0.01513672, -0.01254272, -0.01034546]]), '__version__': '1.0', 'Fs': array([[44100]], dtype=uint16), '__header__': 'MATLAB 5.0 MAT-file, Platform: MACI64, Created on: Tue Feb 5 21:22:28 2013', '__globals__': []}
By typing whos
you get the following line as well as many others:
jingle ndarray 1x255000: 255000 elems, type `float64`, 2040000 bytes (1 Mb)
Test it here:
whos
Variable Type Data/Info ------------------------------------- LAMBDA float 0.94 M int 10 a int 2 alpha float 0.999 b int -10 c module <module 'scipy.constants'<...>\constants\__init__.pyc'> dates ndarray 1306L: 1306 elems, type `float64`, 10448 bytes error ndarray 901L: 901 elems, type `float64`, 7208 bytes error2 ndarray 801L: 801 elems, type `float64`, 6408 bytes h ndarray 100L: 100 elems, type `float64`, 800 bytes jingle ndarray 1Lx255000L: 255000 elems, type `float64`, 2040000 bytes (1 Mb) lbd float 0.9 m module <module 'math' (built-in)> mat_contents dict n=5 maxError float64 2.13162820728e-14 maxError2 float64 0.0 noise ndarray 2000L: 2000 elems, type `float64`, 16000 bytes price ndarray 1306L: 1306 elems, type `float64`, 10448 bytes sigma2 float 0.1 sio module <module 'scipy.io' from '<...>s\scipy\io\__init__.pyc'> x ndarray 2000L: 2000 elems, type `float64`, 16000 bytes x1 ndarray 1000L: 1000 elems, type `float64`, 8000 bytes x1Delayed ndarray 900L: 900 elems, type `float64`, 7200 bytes x2 ndarray 1000L: 1000 elems, type `float64`, 8000 bytes xNoisy ndarray 2000L: 2000 elems, type `float64`, 16000 bytes xx ndarray 1000L: 1000 elems, type `float64`, 8000 bytes y ndarray 1000L: 1000 elems, type `float64`, 8000 bytes y1 ndarray 901L: 901 elems, type `float64`, 7208 bytes y1Delayed ndarray 801L: 801 elems, type `float64`, 6408 bytes y1InputDelayed ndarray 801L: 801 elems, type `float64`, 6408 bytes y2 ndarray 901L: 901 elems, type `float64`, 7208 bytes yConv ndarray 1901L: 1901 elems, type `float64`, 15208 bytes yy ndarray 901L: 901 elems, type `float64`, 7208 bytes
you can see that the vector (signal samples) jingle
is
of size 1x255000. Fs
is the variable containing the
frequency at which the song jingle
has been sampled.
By playing it you will probably recognize part of a famous jingle,
but with a swiss yodeling touch!
We now create a white noise signal
sigma = 0.01
noise = sigma * np.random.randn(255000)
to be added to the signal
jingleNoisy = jingle + noise
jingleNoisy = jingleNoisy.reshape(255000,)
Play it so hear the effect of the noise!
We now apply the Leaky Integrator filter on jingleNoisy
jingleConv = np.convolve(jingleNoisy,h,'valid')
Play the denoised signal y. What do you hear? The noise is partially gone by the voice sounds bass .. what happened?
The effect we have seen on the two tones (higher tone attenuated) and on the jingle signal (bass voice) arise the natural question of what is the effect of the Leaky Integrator filter on the frequency content of the signal. Question that finds a straightforward answer from the convolution theorem, as discussed further in these numerical examples.
Let us also filter a random process with a moving average filter seen in the class. The random process x[n] is obtained by sampling independently at random a normal random variable with zero mean and variance 1.
Random = np.random.normal(0,1,100)
plot(np.arange(Random.size),Random)
xlabel('Samples')
ylabel('Random')
show()
The impulse response of the filter is given by
h[n]={1 if 0≤n≤M−10 otherwise,The following code computes the moving average for various values of M by convolving the input signal with the filter impulse response using the function conv(...).
M = np.linspace(10,50,5) # Parameters
x = np.random.randn(100,1)
from scipy import signal as sc
Random_Conv = [0,1,2,3,4]
for i in range(len(M)):
current_M = M[i]
h = 1.0/current_M*np.ones(current_M).reshape(current_M,1)
Random_Conv[i] = sc.convolve2d(x,h,'valid')
First, note that the output y is of size M+L−1 (Recall that L is the size of h and M is the size of x). Indeed, remember that when convolving a signal with an impulse response, we flip the impulse response in time and shift it progressively on the right until there is no more overlap between the two signals. There is first L−1 points that are formed by a partial overlap of h and x, than M−L+1 points that are formed by a complete overlap of h and x and, again, L−1 points that are formed by a partial overlap of h and x. Let us now display the result for various values of M, the length of the filter.
plot(np.arange(Random_Conv[0].size),Random_Conv[0])
ylabel('M = 10')
xlabel('Samples')
show()
plot(np.arange(Random_Conv[1].size),Random_Conv[1])
ylabel('M = 20')
xlabel('Samples')
show()
plot(np.arange(Random_Conv[2].size),Random_Conv[2])
ylabel('M = 30')
xlabel('Samples')
show()
plot(np.arange(Random_Conv[3].size),Random_Conv[3])
ylabel('M = 40')
xlabel('Samples')
show()
plot(np.arange(Random_Conv[4].size),Random_Conv[4])
ylabel('M = 50')
xlabel('Samples')
show()
We observe that, as M increases, a cycle appears in the process. Whereas the original signal is random and very difficult to predict, the filtered signal appears very predictable as it contains a cycle. Imagine that the original process is the return of a stock, have we find a method to predict it? The answer is (unfortunately) no. The cycle is just a spurious effect known as the Slutzky-Yule effect. Although moving average is a commonly used tool to denoise a signal, we should remember that it creates spurious effects for large M or when applying the filter repeatedly.
The convolution theorem tells us that the Fourier transform of the convolution between two sequences is equal to the product of the Fourier transform of the two sequences.
Working with finite sequences, we need to be careful on the choice of the length of the Fourier transform. We will use here the command for a K point fft:
scipy.fftpack.fft(X,K).
When convoluting a sequence x of length N with a sequence h of length L<N
y = np.convolve(x,h,'valid')
the result is of length N−L+1.
Therefore, we will consider N−L+1 point FFT.
Let's take our noisy two tones xNoisy
(length N=2000), the Leaky
Integrator filter h (length L=100), and plot their spectral
content, that is, the absolute value of their Fourier transform:
# Construct the impulse response of the system of length 100
import numpy as np
import math as m
from scipy import constants as c
M = 10
lbd = float((M-1))/float(M)
h = (1-lbd)*pow(lbd,np.arange(100))
# Generate two tones
y = np.linspace(0,999,1000)
x1 = np.sin((y*2*c.pi*40)/1000)
x2 = np.sin((y*2*c.pi*80)/1000)
x = np.append(x1,x2)
# Generate the noise
sigma3 = 0.1
noise = sigma3 * np.random.randn(2000)
# Add noise to signal
xNoisy = noise + x
N = len(xNoisy)
L = len(h)
# DFT/DFS of noisy signal and impulse response
from scipy import fftpack as f
XNoisy = f.fft(xNoisy,N-L+1)
H = f.fft(h,N-L+1)
normFrequ = np.arange(N-L+1,dtype=float)/(float(N-L+1)) # To plot vs the normalized frequencies
subplot(2,1,1)
plot(normFrequ,abs(XNoisy))
xlabel('Normalized Frequencies')
ylabel('|XNoisy|')
subplot(2,1,2)
plot(normFrequ,abs(H))
xlabel('Normalized Frequencies')
ylabel('|H|')
show()
The absolute value of the product (samplewise) of the two Fourier transforms is computed as:
ABS = H * XNoisy
The plot of the absolute value of the product (samplewise) of the two Fourier transforms
plot(normFrequ,abs(ABS))
xlabel('Normalized Frequencies')
ylabel('|H * XNoisy|')
show()
As expected from the shape of |H|, the spectral content of the second tone, which is at higher frequencies than the first tone, is attenuated. The Leaky Integrator filter has a low pass effect!
Of course, the inverse Fourier transform of ABS gives us the denoised two tones in the time domain
y = scipy.fftpack.ifft(ABS)
Let's now see what happens in the frequency domain with the jingle signal
import scipy.io as sio
mat_contents = sio.loadmat('jingle.mat')
jingle = mat_contents['jingle']
jingle = jingle.reshape(jingle.size,)
# Generate the noise
sigma4 = 0.1
noise = sigma4 * np.random.randn(jingle.size)
# Add noise to signal
jingleNoisy = jingle + noise
# Construct the impulse response of the system of length 100
M = 10
lbd = float(M-1)/float(M)
h = (1-lbd) * pow(lbd,np.arange(100))
L = len(h)
N_J = len(jingle)
# DFT/DFS of noisy signal and impulse response
from scipy import fftpack as f
JingleNoisy = f.fft(jingleNoisy,N_J-L+1)
H = f.fft(h,N_J-L+1)
normFrequJ = np.arange(N_J-L+1,dtype = float)/float(N_J-L+1)
Y_J = H * JingleNoisy
Jingle = f.fft(jingle,N_J-L+1)
# To plot around the origin for visualizationease
normFrequJ_2 = np.arange(50000,dtype=float)/float(N_J-L+1)
The plot of the product (samplewise) of the absolute value of two Fourier transforms |H||JingleNoisy| compared to the absolute value of the Fourier transform of the original signal (around the origin for visualization ease) is depicted below:
Y = H*JingleNoisy
Jingle = f.fft(jingle,N_J-L+1)
normFreq = np.arange(1,50000,dtype=float)/float(N_J-L+1)
subplot(2,1,1)
plot(np.arange(1,Jingle.size+1,dtype=float)/float(Jingle.size),abs(Jingle))
ylabel('|Jingle|')
xlim([0,0.2])
subplot(2,1,2)
plot(np.arange(1,Y.size+1,dtype=float)/float(Y.size),abs(Y))
xlabel('Normalized Frequencies')
ylabel('|Y|')
xlim([0,0.2])
show()
Notice: To run this example, you need some data and code files
that we provide here in compressed format exampleFinanceRisk.zip
.
The goal of this example is to apply the concept of filters seen in the lecture to compute commonly used risk measures in finance. Consider the following problem faced by an investor. He is interested in a certain stock whose evolution of the price p[t] over 5 years sampled at daily frequency is obtained using the following code and depicted in the figure below.
# Parameters
M = 100
LAMBDA = 0.94
import scipy.io as sio
mat_contents = sio.loadmat('data.mat')
# Load data
dates = mat_contents['dates_ts']
dates = dates.reshape(dates.size,)
price = mat_contents['price_ts']
price = price.reshape(price.size,)
#for logarithmic plotting use pylab.semilogy
dates_n = np.arange(dates.size)
semilogy(dates_n,price)
ylim([50,250])
xlim([0,1305])
xlabel('Days (From 11/07/2006)')
ylabel('Log-Price')
show()
More than the price, an investor is interested in the return r[t]of his investment, defined as
r[t]=p[t]−p[t−1]p[t],
and obtained using the following code:
# Compute simple return time-series
return_ts = (price[1:price.size] - price[0:(price.size-1)])/price[1:price.size]
dates_ts = dates[2:dates.size]
length_ts = len(return_ts)
#For small point plotting use pylab.plot with '.' option
plot(np.arange(return_ts.size),return_ts,'.')
xlabel('Days (From 11/07/2006)')
ylabel('Return')
xlim([0,1305])
show()
At the same time, an investor would like to assess the risk of his investment. Since the seminal work of Markowitz, this is done (imperfectly) by measuring the volatility or, in signal processing terms, the standard deviation of returns. The volatility is not directly observable and must be estimated. We now see how the paradigm of filters can be applied to compute several estimates of volatility.
Let us come back to the problem of measuring the volatility. A first way of measuring it is to compute its value over the last M observations of returns. Let us first define the average return ˉr[t] over the last M observations
ˉr[t]=1M∑M−1k=0r(t−k)
In other terms, this is the convolution of the return signal with the moving average filter. The volatility is then given by
σ2[t]=1M∑M−1k=0(r[t−k]−ˉr[t])2
σ[t]=√σ2[t]
Again, this formula can be interpreted in terms of filtering. It is the convolution between the moving average filter with the process that is obtained by squaring the difference between contemporaneous return and the M point moving average return.
# Compute the M-point volatility
import numpy as np
volatility = np.zeros(length_ts)
h = 1.0/float(M)*np.ones(M,dtype=float)
for i in range(M,length_ts):
temp = return_ts[(i-M):i]
moving_average = np.dot(h,temp)
volatility[i] = np.dot(h,pow(temp-moving_average,2))
volatility = np.sqrt(volatility)
plot(np.arange(volatility.size),volatility)
xlabel('Days')
ylabel('M Points Volatility')
xlim([100,1305])
show()
In the previous formula, when computing the average return or volatility, the same weight is given to each observation (1/M). Since more recent observations are more relevant to predict the future, we can consider instead the exponentially weighted volatility, given by the formula
ˉr[t]=λˉr[t−1]+(1−λ)r[t]
σ2[t]=λσ2[t−1]+(1−λ)(r[t]−ˉr[t])2
σ[t]=√σ2[t]
We recognize above the recursive formula of the leaky integrator. Thus, this is again a simple filtering by the impulse response of the leaky integrator. The coefficients of the filter decreases exponentially, hence the name exponentially weighted volatility.
# Compute the exponentially weighted volatility
exp_average = 0
exp_volatility = np.zeros(length_ts + 1)
for i in range(length_ts):
exp_average = LAMBDA * exp_average + (1 - LAMBDA) * return_ts[i]
exp_volatility[i + 1] = LAMBDA * exp_volatility[i] + (1 - LAMBDA)*pow(return_ts[i] - exp_average,2)
exp_volatility = np.sqrt(exp_volatility[1:exp_volatility.size])
plot(np.arange(exp_volatility.size),exp_volatility)
xlabel('Days')
ylabel('Exponential Weighted Volatility')
xlim([0,1305])
show()
Notice: To run this example, you need some data and code files that
we provide here in compressed format exampleDTMF.zip
.
Let us come back to the example of a two-tones signal. Such signal seems very simplistic, but it is the signal underlying a former signaling system known as DTMF. The role of a signaling system in a telephone network is, among others, to setup a call, in particular by determining the number dialed by a user. Before the emergence of touch screen mobile phone, push-button fixed phone were ubiquitous, as some of you may remember. When the user presses a button, the phone generates a two-tones signal. The frequencies composing the signal correspond to the frequency of the column and the frequency of the row where the button lies in, as illustrated by the following table.
1209 Hz | 1336 Hz | 1477 Hz | |
---|---|---|---|
697 Hz | 1 | 2 | 3 |
770 Hz | 4 | 5 | 6 |
852 Hz | 7 | 8 | 9 |
941 Hz | * | 0 | # |
For example, upon pressing a 1, the generated signal contains frequencies 1209 Hz and 697 Hz. We obtain the following continuous time signal:
x(t)=sin(2058πt)+sin(1394πt).
Furthermore, the various frequencies were chosen such that none can be obtained as a result of simple operations of two others. This is meant to avoid interference.
The generated signal is transmitted over an analog transmission line to a switching center where the button that was pressed should be determined. This is done by analyzing the frequencies present in the signal. We have seen that this could be done by using Fourier analysis. This is however not the procedure that is used in DTMF, as it is computationally expensive. Instead, the following digital processing system is used, whose block diagram is represented below.
The signal x(t) is first sampled at a frequency fs so as to obtain the digital signal x[n]. The later is filtered by a series of 7 passband filters placed in parallel. Each of them attenuates all frequencies except the one corresponding to one of the frequencies of DTMF f0.
On the one hand, the frequency f0 is a real-world, continuous frequency, expressed in Hz. On the other hand, the filters are digital and their frequency is a digital, unitless frequency ω0. The two are related by the relation
ω0=2πf0fs.
This link between continuous frequency f0 will be made clear in the next lectures on sampling. In the mean time, let us accept it as it. Finally, a block called decoder outputs the digit depending on the energy of the signals placed at its input. We focus now on the design of one of the filters of the above figure.
One possibility consists of designing an IIR filter. In order to let only a certain frequency ω0 pass, the filter has a pole at the desired ω0. The filter also has a pole at −ω0, so as to obtain a real impulse response. The pole-zero diagram for f0 is depicted below:
Since the poles are located on the unit circle, the filter is not stable. This can be a great problem in practice, as a small error can get amplified and lead to a decoding error. As a solution we can place instead the poles at λejw0 and λe−jw0 with 0<λ<1. The poles-zeros diagram looks like
The z-transform of the filter is then given by
H(z)=1(1−λejw0z−1)(1−λe−jw0z−1)
H(z)=11−2λcos(w0)z−1+λ2z−2
This filter is stable, as all the poles are located inside the unit circle. It is called a resonator. In the time-domain, given the input x[n] the output of the filter y[n] is then given by the following iterative formula
y[n]=x[n]−2λcos(w0)y[n−1]+λ2y[n−2]
The filter is implemented using the following code.
F0 = 697.
FS = 2000.
LAMBDA = 0.98
# Generate discrete-time signal x
n = np.linspace(0.0005,0.1,200)
x = np.sin(2*np.pi*697* n) + np.sin(2*c.pi*1209*n)
#Plot pole diagram of the filter
omega0 = 2*np.pi*float(F0)/float(FS)
pole_1 = np.exp(1j*omega0)
pole_2 = np.exp(-1j*omega0)
#Filter Signal
omega0 = 2*np.pi*float(F0)/float(FS)
b = np.array([1.])
a = np.array([1., -2*np.cos(omega0)*LAMBDA**2, LAMBDA**2])
def my_filter(x,b,a):
import numpy as np
# Extract Parameters
nx = len(x)
nb = len(b)
na = len(a)
# Zero Padding
y = np.zeros(nx+na-1)
x = np.append(np.zeros(nb -1),x)
# Filter
for N in range(nx):
if na > 1:
y[N+na-1] = -(a[na-1:0:-1]*y[N:N+na-1]).sum() + (b[::-1]*x[N:N+nb]).sum()
else:
y[N+na-1] = (b[::-1]*x[N:N+nb]).sum()
# remove padding
y = y[na:]
return y
y1 = my_filter(x,b,a)
We have used the function my_filter(...)
, which
implements a filter written in direct form
H(z)=b(1)+b(2)z−1+…+b(nb)znb−11+a(2)z−1+…+a(na)zna−1
We provide you with an implementation of this function. For example, if we place as input of this filter the signal obtained by pressing a '1', sampled at fs=2000 Hz
plot(np.arange(x.size),x)
xlabel('Samples')
ylabel('x')
show()
we obtain the following output
plot(np.arange(y1.size), y1)
xlabel('Samples')
ylabel('y')
show()
Let us look at the magnitude of the DFT/DFS of this signal
plot(np.arange(y1.size)/float(y1.size), abs(fft.fft(y1)))
xlabel('normalized frequency')
ylabel('|Y1|')
<matplotlib.text.Text at 0xf335c50>
This is exactly what we wanted, a filter that lets a certain frequency pass while attenuating other frequencies. Finally, let us count the number of operations to compute the output of the filter. Each iteration requires 2 multiplications and 2 additions, thus a total of 2N multiplications and 2N additions.
Instead of designing an IIR filter, we could instead choose to
design an FIR filter. The Parks-McClellan algorithm is a standard
procedure that finds the optimal FIR filter given a set of
specifications. Here optimal means that the resulting filter
minimizes the weighted maximum error. To implement our filter, we can use the function
scipy.signal.remez(...)
which is part of the signal processing toolbox.
As before, we would like to design a bandpass filter, but this time we specify certain tolerance ds and dp in the stop band and pass band, respectively.
# Parameters
F0 = 697.
FS = 2000.
# Generate discrete-time signal x
n = np.linspace(0.0005,0.1,200)
x = np.sin(2*np.pi*697* n) + np.sin(2*np.pi*1209*n)
# Parks-McClellan
numbands = 3
bands = (1./float(FS))*np.array([0, F0-100, F0-25, F0+25, F0+100, FS/2.])
df = bands[2] - bands[1]
des = np.array([0, 1, 0])
dp = 0.1
ds = 0.05
k = ds/dp
weight = np.array([1, k, 1])
griddensity = 16
numtaps = np.ceil(1-(10*np.log10(dp*ds)+13)/(2.324*2*np.pi*df)) + 1
from scipy import signal as sc
h = sc.remez(int(numtaps), bands, des, weight=weight, type='bandpass', grid_density=griddensity, maxiter=40)
# plot the frequency response of the filter
freq, response = sc.freqz(h)
ampl = np.abs(response)
fig = plt.figure()
ax1 = fig.add_subplot(111)
ax1.semilogy(freq/(2*np.pi), ampl, 'b-') # freq in Hz
title('Response of optimal filter')
xlabel('Normalized frequency')
plt.show()
Applying the filter to the same input signal x, we obtain the following output
#zzz = my_filter(x, h[::-1], np.array([1.]))
y2 = my_filter(x, h, np.array([1.]))
plot(y2)
show()
where one tones has been attenuated, as illustrated by the norm of its Fourier transform
plot(np.arange(y2.size)/float(y2.size), abs(fft.fft(y2)))
xlabel('normalized frequency')
ylabel('|Y2|')
show()