I couldn't resist - here's using the Hilbert Transform to detect the phase change over 100 samples and extrapolating that to 100,000 samples. The input is a file containing 200 samples.
Note that any signal that is a harmonic of the sample rate divided by 100 could be used... so with a sample rate of 100kS/s the signal you are analyzing could be 50kHz, 33.3kHz, 25KHz, 20kHz.. whatever. All that matters is that after 100 samples there is nominally zero phase offset. (oh, and you will need to change the magic fudge factor in the kernel if you are looking for a different frequency).
Assumes that the input data has already been bandpass filtered to just around the frequencies of interest.
/************************************************************
* poc.c : A really bad Proof of concept
*
* Can you extract tiny frequency changes from small datasets
* using the Hilbert Transform? Yes, it seems you can.
*
* Expects a single command arg that is a list of 200 samples
* Processes it assuming the sample rate is 100kS/s, and
* measures the crawl in phase over 100 samples. This is then
* used to extrapolate to the change of phase after 100,000 samples.
*
* This isn't supposed to be the most correct code, just a hack
* to see if it works.
************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
double data_r[200];
double data_i[200];
double kernel[99];
static void build_kernel(void) {
int i;
for(i = -49; i < 50; i++) {
if((i&1)==1) {
// I pulled the scaling factor out of my data
kernel[i+49] = 1.0/i/3345*2047;
}
}
}
static void hilbert(double *data, double *result) {
double r = 0;
int i;
// Apply the kernel
for(i = -49; i < 50; i++) {
r += data[i]*kernel[i+49];
}
*result = r;
}
int main(int argc, char *argv[]) {
int i;
FILE *f;
build_kernel();
if(2 != argc) {
fprintf(stderr,"Only supply data file name\n");
exit(0);
}
f = fopen(argv[1],"rb");
if(NULL == f) {
fprintf(stderr,"Unable to open file\n");
exit(0);
}
for(i = 0; i < 200; i++) {
if(1 != fscanf(f,"%lf",data_r+i)) {
fprintf(stderr,"Error reading data\n");
exit(0);
}
data_i[i] = 0.0;
}
fprintf(stderr, "Data read\n");
// Just calculate the transform at two places
hilbert(data_r+49, data_i+49);
hilbert(data_r+149, data_i+149);
double phase49,phase149,change;
phase49 = atan2(data_r[49],data_i[49]);
phase149 = atan2(data_r[149],data_i[149]);
change = (phase149-phase49)/100;
/* change' needs to be wrapped into +/- PI, but I haven't done it */
printf("Angle at sample 49 is %10.6f\n",phase49);
printf("Angle at sample 149 is %10.6f\n",phase149);
printf("Change in cycles after 100,000 samples %10.6f\n", change*100000/(2*M_PI));
return 0;
}
Here's the output when the data is 0.1 Hz faster:
Data read
Angle at sample 49 is -0.657761
Angle at sample 149 is -0.658373
Change in cycles after 100,000 samples -0.097432
And 0.1 Hz slower
Data read
Angle at sample 49 is -0.657381
Angle at sample 149 is -0.656791
Change in cycles after 100,000 samples 0.093872
Here's the hack I was using to generate test data
#include <stdio.h>
#include <math.h>
#define SAMPLE_RATE (100000.0)
#define FREQUENCY (9990.0)
#define SCALE (2047)
int main(int argc, char *argv[]) {
int i;
for(i = 0; i < 200; i++) {
double d = sin(i/(SAMPLE_RATE/FREQUENCY)*2*M_PI)*SCALE;
printf("%i\n",(int)(d));
}
return 0;
}
Home
Help
Search
Login
Register
