I have been given a digital bandpass filter - what is it?

I have been given a code snippet for a digital bandpass filter. I have no idea where the author got it from or how he derived it.

% Creates y being the BandPass filtering of input signal x.  The BandPass will
% accept components of wavelength *period* plus or minus about *delta*%.   
% Others should be rejected.
period = 30;
delta = 0.20;
beta = cos(2*pi()/20);
gamma = 1/cos(4*pi()*delta/period);
alpha = gamma - sqrt(gamma*gamma - 1);
y(1,1)=0;
y(2,1)=0;
for i = 3:3000 % We only filter first 3000 rows of x
    y(i,1) = 0.5*(1-alpha)*(x(i,1)-x(i-2,1)) + beta*(1+alpha)*y(i-1,1) - alpha*y(i-2,1);
end

If used with Period=30 and Delta=0.20 it should capture about 20% of the cyclic behaviour around a wavelength of 30 data points.

Do any of you know the mathematical derivation or family that he has used? For example is this a Butterworth, Chebyshev Type I Filter etc?

Answer 1

[Wayne King]

At first glance, it doesn't look familiar, but do you have the Signal Processing Toolbox? You can easily design bandpass filters using fdesign.bandpass

Answer 2

[Daniel Shub]

It is not a notation I am familiar with. But it looks like beta, gamma, and alpha are scalars and Q is unused. If that is correct then letting y = BP and x = Price gives:

y = ax-ax[2]+by[1]-cy[2]

That almost looks like a filter, but I don't understand what ax means, are you sure it is not Price[1]-Price[2]? Also I would expect it to look like:

y[n] = 0*x[n]+ax[n-1]-ax[n-2]-by[n-1]-c[n-2]

[Stewart Charles]

Thank you for your answer. Sorry, I have just realised that there was excess/unclear details in this.

Yes, Q is unused sorry.

To clear up the format of the line in question:

y = .5*(1 - alpha)*(x[n] - x[n-2]) + beta*(1 + alpha)*y[n-1] - alpha*y[n-2];

[Daniel Shub]

Are you sure it is not y[n] = ...? Look at the documentation of the filter command. You should be able to easily get this into a format that MATLAB is happy with and then you can plot its magnitude and phase properties.

[Stewart Charles]

Hi Daniel,

Yes, it is y[n]=?.

I'm happy enough getting it into MATLAB format. I was wondering if anyone recognised what type of IIR filter this was or how he came up with the coefficients?

[Daniel Shub]

Then I would suggest you edit the question so that it has valid MATLAB code in standard MATLAB format. For example, cos instead of Cosine and filter coefficients given as b and a. Then clearly ask if anyone knows what type of IIR filter it is and where the coefficients came from. That is a much clear question than "does this look familiar to you".

[Stewart Charles]

Many thanks Daniel, I will do that. Apologies for that.

[Stewart Charles]

Daniel,

I have now updated the question after porting the code into Matlab.

Stewart

[Daniel Shub]

It looks much better and is now a much harder question. I will have a think about it.

[Stewart Charles]

Thank you Wayne and Daniel.

I originally asked the question.

I now understand that this is a Chebyshev II bandpass.

I'm currently using this formula for finding cycles in quite noisy data that is regularly discretely sampled. I use the period to set the wavelength of the cycle I'm looking for, eg a 30 bar cycle (if sampled every minute then this would mean a 30 minute cycle, peak to peak). I usually keep the delta at 0.20. If it is too small, the filter rings. If it is too wide, the output isn't particularly cyclic.

If I may, could I ask a couple more questions.

1) Are there any other IIR bandpass filters I could use for this which I could explore? I've looked high and low on the internet and found much theory but very few (plug in the numbers and use it) formulas.

2) I'm interested in taking the output of this filter, which obviously looks roughly sinusoidal, and converting it to polar coordinates. Eg, an estimated amplitude and phase at any given time. A bit like what engineers do with imaginary numbers. I have tried doing this by taking the slope of the output of the bandpass filter, along with the output itself, then using a big of trig to construct a hypotenuse (the amplitude) and angle (phase).*

So, is there a better way to do this? Which leads me to the next question:

3) Digital filters typically have a lag (of phase). This can be a function of the frequency, obviously). However, would it be possible to construct a digital bandpass filter to complement the one above, which produces an output which is cyclic but lagging by 90 degrees? That way I could take both and use pythagoras to estimate the ampitude and trig to get the phase. Obviously such a process only works when you have a decent cycle in the data to actually estimate.

Any thoughts would be valued. I apologies as my engineering maths is really rusty.

Stewart

• Since the slope of a sine wave is the cosine of the same phase. That means the slope of a cycle should give an approximation of the imaginary component of the cycle.

Answer 3

[Wayne King]

Hi Stewart, this is a standard 2nd order linear constant coefficient difference equation. You can mimic this filter easier in MATLAB with

B = [0.5*(1-alpha) -0.5*(1-alpha)];
A = [1 -beta*(1+alpha) alpha];
% view the filter magnitude response
fvtool(B,A);
% filter data
output = filter(B,A,input);

The filter is a bandpass filter with a narrow passband at pi/10 radians/sample. What period that actually corresponds to depends on the sampling frequency of the data.

To answer your question, the design equations look like a Chebyshev type II approximation to me.

You can design a Chebyshev type II filter using fdesign.bandpass and specifying the design method as 'cheby2'

[Daniel Shub]

From Stewart's comments to me and his revised question, it seems he wants to know if this is a Butterworth filer or a Chebshev or some other common IIR filter.

[Wayne King]

@Daniel Thank you for pointing out that omission, I have added that information.

[Stewart Charles]

Thank you Wayne and Daniel.

I originally asked the question.

I now understand that this is a Chebyshev II bandpass.

I'm currently using this formula for finding cycles in quite noisy data that is regularly discretely sampled. I use the period to set the wavelength of the cycle I'm looking for, eg a 30 bar cycle (if sampled every minute then this would mean a 30 minute cycle, peak to peak). I usually keep the delta at 0.20. If it is too small, the filter rings. If it is too wide, the output isn't particularly cyclic.

If I may, could I ask a couple more questions.

1) Are there any other IIR bandpass filters I could use for this which I could explore? I've looked high and low on the internet and found much theory but very few (plug in the numbers and use it) formulas.

2) I'm interested in taking the output of this filter, which obviously looks roughly sinusoidal, and converting it to polar coordinates. Eg, an estimated amplitude and phase at any given time. A bit like what engineers do with imaginary numbers. I have tried doing this by taking the slope of the output of the bandpass filter, along with the output itself, then using a big of trig to construct a hypotenuse (the amplitude) and angle (phase).*

So, is there a better way to do this? Which leads me to the next question:

3) Digital filters typically have a lag (of phase). This can be a function of the frequency, obviously). However, would it be possible to construct a digital bandpass filter to complement the one above, which produces an output which is cyclic but lagging by 90 degrees? That way I could take both and use pythagoras to estimate the ampitude and trig to get the phase. Obviously such a process only works when you have a decent cycle in the data to actually estimate.

Any thoughts would be valued. I apologies as my engineering maths is really rusty.

Stewart

• Since the slope of a sine wave is the cosine of the same phase. That means the slope of a cycle should give an approximation of the imaginary component of the cycle.