add 60Hz/50Hz notch filter

Dependencies:   mbed

Files at this revision

API Documentation at this revision

Comitter:
nagasm
Date:
Tue Dec 16 08:28:20 2014 +0000
Child:
1:8fa2f521009a
Commit message:
add 60Hz/50Hz notch filter

Changed in this revision

FIR_LPF.hpp Show annotated file Show diff for this revision Revisions of this file
main.cpp Show annotated file Show diff for this revision Revisions of this file
mbed.bld Show annotated file Show diff for this revision Revisions of this file
notch.hpp Show annotated file Show diff for this revision Revisions of this file
sub.hpp Show annotated file Show diff for this revision Revisions of this file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/FIR_LPF.hpp	Tue Dec 16 08:28:20 2014 +0000
@@ -0,0 +1,72 @@
+// from CQpub0 Mikami / FIR_LPF_Direct
+//  http://developer.mbed.org/users/CQpub0Mikami/code/FIR_LPF_Direct
+//
+//                                   Band1       Band2
+// Lower band edge frequency (kHz) 0.000000    0.600000
+// Upper band edge frequency (kHz) 0.500000    5.000000
+// Gain                            1.000000    0.000000
+// Weight                          1.000000    1.000000
+// Deviation                       0.009734    0.009734
+// Deviation [dB]                  0.084139  -40.234188
+
+const float fc[order+1] = {   
+    -3.566292E-03f,  2.335185E-03f,  1.917338E-03f,  1.681921E-03f,  1.522689E-03f,
+     1.364066E-03f,  1.157961E-03f,  8.803014E-04f,  5.296940E-04f,  1.236180E-04f,
+    -3.047488E-04f, -7.123405E-04f, -1.051992E-03f, -1.279654E-03f, -1.360180E-03f,
+    -1.273725E-03f, -1.018960E-03f, -6.155554E-04f, -1.029986E-04f,  4.621914E-04f,
+     1.013316E-03f,  1.480161E-03f,  1.798055E-03f,  1.916527E-03f,  1.805744E-03f,
+     1.463100E-03f,  9.143371E-04f,  2.130115E-04f, -5.640853E-04f, -1.326518E-03f,
+    -1.978121E-03f, -2.430574E-03f, -2.613837E-03f, -2.486610E-03f, -2.043447E-03f,
+    -1.317959E-03f, -3.802110E-04f,  6.680631E-04f,  1.704594E-03f,  2.600376E-03f,
+     3.235972E-03f,  3.515095E-03f,  3.379365E-03f,  2.817894E-03f,  1.872141E-03f,
+     6.322899E-04f, -7.684087E-04f, -2.167773E-03f, -3.393601E-03f, -4.284042E-03f,
+    -4.705797E-03f, -4.577073E-03f, -3.875274E-03f, -2.648668E-03f, -1.012065E-03f,
+     8.616354E-04f,  2.758435E-03f,  4.447251E-03f,  5.706155E-03f,  6.350451E-03f,
+     6.257734E-03f,  5.386359E-03f,  3.786968E-03f,  1.602121E-03f, -9.435526E-04f,
+    -3.565645E-03f, -5.948320E-03f, -7.781761E-03f, -8.797770E-03f, -8.807067E-03f,
+    -7.725336E-03f, -5.594939E-03f, -2.585050E-03f,  1.010693E-03f,  4.809328E-03f,
+     8.361672E-03f,  1.121444E-02f,  1.295107E-02f,  1.324389E-02f,  1.191081E-02f,
+     8.934171E-03f,  4.488148E-03f, -1.062401E-03f, -7.182842E-03f, -1.321242E-02f,
+    -1.841506E-02f, -2.204769E-02f, -2.343209E-02f, -2.202840E-02f, -1.749386E-02f,
+    -9.730157E-03f,  1.092382E-03f,  1.453105E-02f,  2.989524E-02f,  4.629194E-02f,
+     6.269474E-02f,  7.802411E-02f,  9.123583E-02f,  1.014078E-01f,  1.078185E-01f,
+     1.100078E-01f,  1.078185E-01f,  1.014078E-01f,  9.123583E-02f,  7.802411E-02f,
+     6.269474E-02f,  4.629194E-02f,  2.989524E-02f,  1.453105E-02f,  1.092382E-03f,
+    -9.730157E-03f, -1.749386E-02f, -2.202840E-02f, -2.343209E-02f, -2.204769E-02f,
+    -1.841506E-02f, -1.321242E-02f, -7.182842E-03f, -1.062401E-03f,  4.488148E-03f,
+     8.934171E-03f,  1.191081E-02f,  1.324389E-02f,  1.295107E-02f,  1.121444E-02f,
+     8.361672E-03f,  4.809328E-03f,  1.010693E-03f, -2.585050E-03f, -5.594939E-03f,
+    -7.725336E-03f, -8.807067E-03f, -8.797770E-03f, -7.781761E-03f, -5.948320E-03f,
+    -3.565645E-03f, -9.435526E-04f,  1.602121E-03f,  3.786968E-03f,  5.386359E-03f,
+     6.257734E-03f,  6.350451E-03f,  5.706155E-03f,  4.447251E-03f,  2.758435E-03f,
+     8.616354E-04f, -1.012065E-03f, -2.648668E-03f, -3.875274E-03f, -4.577073E-03f,
+    -4.705797E-03f, -4.284042E-03f, -3.393601E-03f, -2.167773E-03f, -7.684087E-04f,
+     6.322899E-04f,  1.872141E-03f,  2.817894E-03f,  3.379365E-03f,  3.515095E-03f,
+     3.235972E-03f,  2.600376E-03f,  1.704594E-03f,  6.680631E-04f, -3.802110E-04f,
+    -1.317959E-03f, -2.043447E-03f, -2.486610E-03f, -2.613837E-03f, -2.430574E-03f,
+    -1.978121E-03f, -1.326518E-03f, -5.640853E-04f,  2.130115E-04f,  9.143371E-04f,
+     1.463100E-03f,  1.805744E-03f,  1.916527E-03f,  1.798055E-03f,  1.480161E-03f,
+     1.013316E-03f,  4.621914E-04f, -1.029986E-04f, -6.155554E-04f, -1.018960E-03f,
+    -1.273725E-03f, -1.360180E-03f, -1.279654E-03f, -1.051992E-03f, -7.123405E-04f,
+    -3.047488E-04f,  1.236180E-04f,  5.296940E-04f,  8.803014E-04f,  1.157961E-03f,
+     1.364066E-03f,  1.522689E-03f,  1.681921E-03f,  1.917338E-03f,  2.335185E-03f,
+    -3.566292E-03f
+};
+
+float FIR_calc1(float data){
+    int i;
+    float acc = 0;
+    xn1[0] = data;
+    for (i=0; i<=order; i++) acc = acc + fc[i]*xn1[i];
+    for (i=order; i>0; i--) xn1[i] = xn1[i-1];
+    return(acc);
+}
+
+float FIR_calc2(float data){
+    int i;
+    float acc = 0;
+    xn2[0] = data;
+    for (i=0; i<=order; i++) acc = acc + fc[i]*xn2[i];
+    for (i=order; i>0; i--) xn2[i] = xn2[i-1];
+    return(acc);
+}
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main.cpp	Tue Dec 16 08:28:20 2014 +0000
@@ -0,0 +1,66 @@
+#include "mbed.h"
+#include "sub.hpp"
+#include "FIR_LPF.hpp"
+#include "notch.hpp"
+
+int main(){
+    common_setup();
+    xbee.baud(38400);
+    xbee.attach(&rx_fifoset, xbee.RxIrq);
+    timer_setup.attach_us(&timer_interrupt, 5); // 5usec
+    coef_set(60);   // Ham notch filter set
+    while(1){
+        tx_fifo_check();
+        if(timer_value[1] > 19){ // 0.1msec sampling
+            timer_value[1] = 0;
+            float data1 = (float)gain * (analog_value2.read() - 0.5f);  // A/D in (3)
+            if (data1 < 0) data1 = -data1; // always detection ON
+            if (fir_lpf != 0) data1 = FIR_calc1(data1);      // FIR calc (1) call
+            data1 = notch_filter1(data1);
+            if(timer_value[2] > 2999){  // 15msec
+                timer_value[2] = 0;
+                if(max_count != 0) data1 = move_mean_calc1(data1);
+                tx_message((uint16_t)((data1 + 1.0f) * 2047)<<4);
+            }
+        }
+        if(timer_value[3] > 19){ // 0.1msec sampling
+            timer_value[3] = 0;
+            float data2 = (float)gain * (analog_value3.read() - 0.5f);  // A/D in (4)
+            if (data2 < 0) data2 = -data2; // always detection ON
+            if (fir_lpf != 0) data2 = FIR_calc2(data2);      // FIR calc (1) call
+            data2 = notch_filter2(data2);
+            if(timer_value[4] > 2999){  // 15msec
+                timer_value[4] = 0;
+                if(max_count != 0) data2 = move_mean_calc2(data2);
+                tx_message(0x400000 + ((uint16_t)((data2 + 1.0f) * 2047)<<4));
+            }
+        }
+        if(timer_value[0] > 199999){ // 1000msec
+            timer_value[0] = 0;
+            myled = !myled;
+        }
+        if(rx_fifo_check() == 1){
+            int sum = 0;
+            for (int i=0; i<6; i++) sum += conv_hex(raw_data[i])<<(4*(5-i));
+            tx_message(sum); // Echo Back 
+            if(sum>>16 == 0x80){
+                switch((sum & 0xff00)>>8){
+                    case 0x00:
+                        fir_lpf = sum & 0x01;
+                        break;
+                    case 0x01:
+                        gain = sum & 0x0f;
+                        break;
+                    case 0x02:
+                        max_count = sum & 0x7f;
+                        if (max_count>100) max_count = 100;
+                        sum_clear();
+                        break;
+                    case 0x03:
+                        coef_set(sum & 0x3f);
+                        break;
+                }
+            }
+        }        
+    }
+}
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/mbed.bld	Tue Dec 16 08:28:20 2014 +0000
@@ -0,0 +1,1 @@
+http://mbed.org/users/mbed_official/code/mbed/builds/4fc01daae5a5
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/notch.hpp	Tue Dec 16 08:28:20 2014 +0000
@@ -0,0 +1,68 @@
+/*
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+int main(){
+    int fs = 10000;
+    double  B0, a1, a2, b1, c0, F0, T;
+    F0 = 60;    // or 50
+    B0 = 100.;
+    T  = 1./fs;
+    a1 = 2. * exp(-M_PI * B0 * T) * cos(2. * M_PI * F0 * T);
+    a2 = -exp(-2. * M_PI * B0 * T);
+    b1 = -2. * cos(2. * M_PI * F0 * T);
+    c0 = (1-a1-a2)/(2+b1);
+    printf("a1 = %f\n", a1);
+    printf("a2 = %f\n", a2);
+    printf("b1 = %f\n", b1);
+    printf("c0 = %f\n", c0);
+    return 0;
+}
+60Hz
+    a1 = 1.936768
+    a2 = -0.939101
+    b1 = -1.998579
+    c0 = 1.642174
+50Hz
+    a1 = 1.937188
+    a2 = -0.939101
+    b1 = -1.999013
+    c0 = 1.938304
+*/
+
+const float _60Hz_a1 = 1.936768, _60Hz_a2 = -0.939101, _60Hz_b1 = -1.998579, _60Hz_c0 = 1.642174;
+const float _50Hz_a1 = 1.937188, _50Hz_a2 = -0.939101, _50Hz_b1 = -1.999013, _50Hz_c0 = 1.938304;
+float a1, a2, b1, c0;
+float y1[3], y2[3];
+
+void coef_set(int herz){
+    if(herz < 55){
+        a1 = _50Hz_a1;
+        a2 = _50Hz_a2;
+        b1 = _50Hz_b1;
+        c0 = _50Hz_c0;
+    }
+    else{
+        a1 = _60Hz_a1;
+        a2 = _60Hz_a2;
+        b1 = _60Hz_b1;
+        c0 = _60Hz_c0;
+    }
+    for (int i=0; i<3; i++) y1[i] = y2[i] = 0;
+}
+
+float notch_filter1(float data){
+    y1[0] = data + a1*y1[1] + a2*y1[2];
+    float reault = y1[0] + b1*y1[1] + y1[2];
+    y1[2] = y1[1];
+    y1[1] = y1[0];
+    return(reault);
+}
+
+float notch_filter2(float data){
+    y2[0] = data + a1*y2[1] + a2*y2[2];
+    float reault = y2[0] + b1*y2[1] + y2[2];
+    y2[2] = y2[1];
+    y2[1] = y2[0];
+    return(reault);
+}
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/sub.hpp	Tue Dec 16 08:28:20 2014 +0000
@@ -0,0 +1,107 @@
+unsigned char rxFIFO[256], txFIFO[256], raw_data[6];
+unsigned char rx_top, rx_end, tx_top, tx_end, phase;
+float mean_sum1, mean_sum2, ad_data1[101], ad_data2[101];
+int timer_value[6], max_count, fir_lpf, gain, ad_pointer1, ad_pointer2;
+const int order = 200;
+float xn1[order+1], xn2[order+1];
+
+RawSerial xbee(PA_2, PA_3);
+Ticker timer_setup;
+AnalogIn analog_value0(A0);
+AnalogIn analog_value1(A1);
+AnalogIn analog_value2(A2);
+AnalogIn analog_value3(A3);
+DigitalOut myled(LED1);
+
+void sum_clear(){
+    int i;
+    for (i=0; i<101; i++) ad_data1[i] = ad_data2[i] = 0;
+    ad_pointer1 = ad_pointer2 = mean_sum1 = mean_sum2 = 0;
+}
+
+void common_setup(){
+    int i;
+    rx_top = rx_end = tx_top = tx_end = phase = 0; 
+    max_count = 100;
+    fir_lpf = 0;
+    gain = 5;
+    for (i=0; i<3; i++) timer_value[i] = 0;
+    timer_value[3] = 10;
+    timer_value[4] = 1500;
+    for (i=0; i<=order; i++) xn1[i] = xn2[i] = 0.0;
+    sum_clear();
+}
+
+void timer_interrupt(){
+    int i;
+    for (i=0; i<6; i++) timer_value[i]++;
+}
+
+void tx_fifo_check(){
+    if(xbee.writeable() == 1){
+        if(tx_top != tx_end){
+            xbee.putc(txFIFO[tx_end]);
+            ++tx_end &= 255;
+        }
+    }
+}
+
+int rx_fifo_check(){
+    unsigned char data;
+    if(rx_top != rx_end){
+        data = rxFIFO[rx_end];
+        ++rx_end &= 255;
+        if (data < 33){
+            phase = 0;
+            return(1);
+        }
+        raw_data[phase] = data;
+        if(++phase > 5) phase = 0;
+        return(0);
+    }
+    return(0);
+}
+
+void rx_fifoset(void){
+    rxFIFO[rx_top] = xbee.getc();
+    ++rx_top &= 255;
+}
+
+void tx_fifoset(unsigned char data){
+    txFIFO[tx_top] = data;
+    ++tx_top &= 255;
+}
+
+unsigned char hex_conv(unsigned char data){
+    data &= 15;
+    if(data < 10) return(data+48);
+    else return(data+55);
+}
+
+unsigned char conv_hex(unsigned char data){
+    if((data > 47) && (data < 58)) return(data-48);
+    else if((data > 64) && (data < 71)) return(data-55);
+    return(0);
+}
+
+void tx_message(int data){
+    int i;
+    for (i=0; i<6; i++) tx_fifoset(hex_conv((data>>(4*(5-i))) & 15));
+    tx_fifoset(13);
+}
+
+float move_mean_calc1(float data){
+    mean_sum1 = mean_sum1 - ad_data1[ad_pointer1] + data;
+    ad_data1[ad_pointer1] = data;
+    ad_pointer1++;
+    if(ad_pointer1 == max_count) ad_pointer1 = 0;
+    return(mean_sum1 / (float)max_count);
+}
+
+float move_mean_calc2(float data){
+    mean_sum2 = mean_sum2 - ad_data2[ad_pointer2] + data;
+    ad_data2[ad_pointer2] = data;
+    ad_pointer2++;
+    if(ad_pointer2 == max_count) ad_pointer2 = 0;
+    return(mean_sum2 / (float)max_count);
+}