arfreq.c

来自「dsp AD公司ADSP21的代码,里面有FFT FIR IIR EQULIZE」· C语言 代码 · 共 59 行

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include "rtdspc.h"

/* ARFREQ.C - take real data in one record and determine the 1st
order AR frequency estimate versus time.  Uses a Hilbert
transform to convert the real signal to complex representation
*/

main()
{
/* 35 point Hilbert transform FIR filter cutoff at 0.02 and 0.48
   +/- 0.5 dB ripple in passband, zeros at 0 and 0.5 */

    static float  fir_hilbert35[35] = {
    0.038135,    0.000000,    0.024179,    0.000000,    0.032403,
    0.000000,    0.043301,    0.000000,    0.058420,    0.000000,
    0.081119,    0.000000,    0.120167,    0.000000,    0.207859,
    0.000000,    0.635163,    0.000000,   -0.635163,    0.000000,
   -0.207859,    0.000000,   -0.120167,    0.000000,   -0.081119,
    0.000000,   -0.058420,    0.000000,   -0.043301,    0.000000,
   -0.032403,    0.000000,   -0.024179,    0.000000,   -0.038135
                          };

    static float hist[34];
    int i,winlen;
    float sig_real,sig_imag,last_real,last_imag;
    float cpi,xr,xi,freq;

    cpi = 1.0/(2.0*PI);
    winlen = 32;

    last_real = 0.0;
    last_imag = 0.0;
    for(;;) {
/* determine the phase difference between sucessive samples */
      xr = 0.0;
      xi = 0.0;
      for(i = 0 ; i < winlen ; i++) {
        sig_imag = fir_filter(getinput(),fir_hilbert35,35,hist);
        sig_real = hist[16];
        xr += sig_real * last_real;
        xr += sig_imag * last_imag;
        xi += sig_real * last_imag;
        xi -= sig_imag * last_real;
        last_real = sig_real;
        last_imag = sig_imag;
      }
/* make sure the result is valid, give 0 if not */
      if(fabs(xr) > 1e-10)
        freq = cpi*atan2(xi,xr);
      else
        freq = 0.0;
      sendout(freq);
    }
}