Digital Band pass filter arduino 1100-1300 Hz

Question

For a project I need to design an digital band pass filter. I started of with this in Matlab:

wch=1100*2*pi;
wcl=1300*2*pi;
sys=tf([1/wch 0],[1/(wch*wcl) 1/wch+1/wcl 1])
sysd=c2d(sys,9615,'tustin')
so I got this values:
3.01e-08 z^2 - 3.01e-08/(z^2 + 2 z + 1)

Then I rewrote it to:

y[k] = 3.01*10^-8*u[k]-3.01*10^-8*u[k-2] -y[k-2]-2y[k-1]

and so I put hoped that this arduino script should do the Job:

#include <TimerOne.h>
float out;
float second_filt(float x0) {
 static float y1=0; // y[n-1]
 static float y2=0; // y[n-2]
 static float x1=0; // x[n-1]
 static float x2=0; // x[n-2]
 float y0 = 3.01*10^(-8)*x0-3.01*10^(-8)*x2-2*y1-1*y2 ;
 // Shift values
 y2 = y1; // Shift y[n-1] to y[n-2]
 y1 = y0; // Shift y[n] to y[n-1]
 x2 = x1; // Shift x[n-1] to x[n-2]
 x1 = x0; // Shift x[n] to x[n-1]
 return y0;
}
void setup() {
 Timer1.initialize(32); //Initialize timer with a period of 32 us
 Timer1.pwm(9,0); // Use pin 9 for PWM
}
void loop() {
 if (isTime()) {
 out=second_filt(analogRead(0));
 Timer1.setPwmDuty(9,out); //Sets the duty cycle to the analog value on pin 0
 }
}
bool isTime(){ //Checks if 300 us elapsed since the last call
 static unsigned long t0=0;
 static unsigned long t1=0;
 t1 = micros();
 bool res = (t1-t0)>300;
 if (res) {
 t0 = t1;
 }
 return res;
}

Now I don't know what the error: 'invalid operands of types 'double' and 'double' to binary 'operator^' means and I am not sure if this is actually gonna work. Pleas help me out I have been busy with this for days.

• Chris Stratton: You may have trouble getting this to work on an ATmega-based Arduino, especially if floating point which must be emulated in software. If you encounter difficulty you might consider one of the more recent Teensy series boards which have a more capable processor and a faster ADC as well.

Answer 1

You cannot use ^ as mathematical power symbol in C. Ini C it means the exclusive-or bit operation (hence the error showing it is a bit operation).

For the mathematical power function, you should use the math C library. See library explanation. To include it, use:

#include <math.h>

The power function has the following prototype:

double pow(double x, double y) 
Returns x raised to the power of y.

However if you only need to write a float/double in exponential format, you can use

3.01e-8

Answer 2

You have more than one problem here...

3.01e-08 z^2 - 3.01e-08/(z^2 + 2 z + 1)

I don't quite understand your Mathlab code, but there is something wrong with this transfer function. The denominator can be rewritten as (z + 1)². Thus there is a double pole at z = −1, i.e. at the Nyquist frequency. An infinite gain at Nyquist is not what you expect from a bandpass filter.

y[k] = 3.01*10^-8*u[k] - 3.01*10^-8*u[k-2] - y[k-2] - 2y[k-1]

This is implementing the transfer function 3.01 ×ばつ 10⁻⁸ (z² − 1)/(z² + 2 z + 1), which is not the same as the one written above: mind the parentheses.

float y0 = 3.01*10^(-8)*x0-3.01*10^(-8)*x2-2*y1-1*y2;

As explained in Michel Keijzers’ answer, the constant should be written 3.01e-8. Also, you can save one multiplication by writing it as

float y0 = 3.01e-8*(x0 - x2) - 2*y1 - y2;

if (isTime()) { out=second_filt(analogRead(0)); ... }

You cannot expect to get a uniform sampling time with analogRead(), as your timing will depend on how long your program takes to get trough the loop, which is data dependent and affected by interrupts.

sysd=c2d(sys,9615,'tustin')

I assume 9615 is the expected sampling rate, which translates to 104 μs between samples. This is the time it takes for the ADC to perform one reading. You cannot achieve this rate if you do computations in between the readings.

Here your best option is to configure the ADC in free-running mode. In this mode it will take one sample every 104 μs, and you will achieve the desired sampling rate. You should also do the filtering in parallel with the readings, and have the CPU and the ADC work concurrently instead of sequentially. For implementing this, you will first have to read the datasheet of the ATmega328P powering your Uno.