match: FFT K-window

pull/4/head
Zilog80 2017-12-27 01:48:40 +01:00
rodzic 6e553631fd
commit c7b6b9640d
5 zmienionych plików z 73 dodań i 43 usunięć

Wyświetl plik

@ -21,7 +21,7 @@ typedef unsigned int ui32_t;
typedef short i16_t;
typedef int i32_t;
//#include "demod.h"
#include "demod_dft.h"
static unsigned int sample_in, sample_out, delay;
@ -40,20 +40,23 @@ static float *xs = NULL,
*qs = NULL;
static float dc_ofs = 0.0;
static float dc = 0.0;
/* ------------------------------------------------------------------------------------ */
#include <complex.h>
static int LOG2N, N_DFT;
static double complex *ew;
static float complex *ew;
static double complex *Fm, *X, *Z, *cx;
static float complex *Fm, *X, *Z, *cx;
static float *xn;
static void dft_raw(double complex *Z) {
static void dft_raw(float complex *Z) {
int s, l, l2, i, j, k;
double complex w1, w2, T;
float complex w1, w2, T;
j = 1;
for (i = 1; i < N_DFT; i++) {
@ -73,8 +76,8 @@ static void dft_raw(double complex *Z) {
for (s = 0; s < LOG2N; s++) {
l2 = 1 << s;
l = l2 << 1;
w1 = (double complex)1.0;
w2 = ew[s]; // cexp(-I*M_PI/(double)l2)
w1 = (float complex)1.0;
w2 = ew[s]; // cexp(-I*M_PI/(float)l2)
for (j = 1; j <= l2; j++) {
for (i = j; i <= N_DFT; i += l) {
k = i + l2;
@ -87,30 +90,31 @@ static void dft_raw(double complex *Z) {
}
}
static void dft(float *x, double complex *Z) {
static void dft(float *x, float complex *Z) {
int i;
for (i = 0; i < N_DFT; i++) Z[i] = (double complex)x[i];
for (i = 0; i < N_DFT; i++) Z[i] = (float complex)x[i];
dft_raw(Z);
}
static void Nidft(double complex *Z, double complex *z) {
static void Nidft(float complex *Z, float complex *z) {
int i;
for (i = 0; i < N_DFT; i++) z[i] = conj(Z[i]);
dft_raw(z);
// idft():
// for (i = 0; i < N_DFT; i++) z[i] = conj(z[i])/(double)N_DFT; // hier: z reell
// for (i = 0; i < N_DFT; i++) z[i] = conj(z[i])/(float)N_DFT; // hier: z reell
}
/* ------------------------------------------------------------------------------------ */
int getCorrDFT(int K, unsigned int pos, float *maxv, unsigned int *maxvpos) {
//int K = N/2; // N+N/2=960 < 1024 = N_DFT/2
int getCorrDFT(int abs, int K, unsigned int pos, float *maxv, unsigned int *maxvpos) {
int i;
int mp = -1;
float mx = 0.0;
float xnorm = 1;
unsigned int mpos = 0;
dc = 0.0;
if (N + K > N_DFT/2 - 2) return -1;
if (sample_in < delay+N+K) return -2;
@ -122,15 +126,27 @@ int getCorrDFT(int K, unsigned int pos, float *maxv, unsigned int *maxvpos) {
dft(xn, X);
dc = get_bufmu(pos-sample_out); //oder: dc = creal(X[0])/N_DFT;
for (i = 0; i < N_DFT; i++) Z[i] = X[i]*Fm[i];
Nidft(Z, cx);
for (i = N; i < N+K; i++) {
if (creal(cx[i]) > mx) { // imag(cx)=0
mx = creal(cx[i]);
mp = i;
if (abs) {
for (i = N; i < N+K; i++) {
if (fabs(creal(cx[i])) > fabs(mx)) { // imag(cx)=0
mx = creal(cx[i]);
mp = i;
}
}
}
else {
for (i = N; i < N+K; i++) {
if (creal(cx[i]) > mx) { // imag(cx)=0
mx = creal(cx[i]);
mp = i;
}
}
}
if (mp == N || mp == N+K-1) return -4; // Randwert
@ -244,6 +260,11 @@ float get_bufvar(int ofs) {
return var;
}
float get_bufmu(int ofs) {
float mu = xs[(sample_out+M + ofs) % M]/Nvar;
return mu;
}
int f32buf_sample(FILE *fp, int inv, int cm) {
float s = 0.0;
float xneu, xalt;
@ -252,7 +273,7 @@ int f32buf_sample(FILE *fp, int inv, int cm) {
if (f32read_sample(fp, &s) == EOF) return EOF;
if (inv) s = -s;
bufs[sample_in % M] = s;
bufs[sample_in % M] = s - dc_ofs;
xneu = bufs[(sample_in ) % M];
xalt = bufs[(sample_in+M - Nvar) % M];
@ -315,11 +336,14 @@ static int read_bufbit(int symlen, char *bits, unsigned int mvp, int reset) {
return 0;
}
int headcmp(int symlen, char *hdr, int len, unsigned int mvp) {
int headcmp(int symlen, char *hdr, int len, unsigned int mvp, int inv, int option_dc) {
int errs = 0;
int pos;
int step = 1;
char sign = 0;
if (symlen != 1) step = 2;
if (inv) sign=1;
for (pos = 0; pos < len; pos += step) {
read_bufbit(symlen, rawbits+pos, mvp+1-(int)(len*samples_per_bit), pos==0);
@ -327,10 +351,14 @@ int headcmp(int symlen, char *hdr, int len, unsigned int mvp) {
rawbits[pos] = '\0';
while (len > 0) {
if (rawbits[len-1] != hdr[len-1]) errs += 1;
if ((rawbits[len-1]^sign) != hdr[len-1]) errs += 1;
len--;
}
if (option_dc && errs < 3) {
dc_ofs += dc;
}
return errs;
}
@ -407,8 +435,9 @@ int init_buffers(char hdr[], int hLen, int shape) {
float *m = NULL;
N = hLen * samples_per_bit;
M = 3*N; // >= N
N = hLen * samples_per_bit + 0.5;
M = 3*N;
if (samples_per_bit < 6) M = 6*N;
Nvar = N; //N/2; // = N/k
bufs = (float *)calloc( M+1, sizeof(float)); if (bufs == NULL) return -100;
@ -459,15 +488,15 @@ int init_buffers(char hdr[], int hLen, int shape) {
}
delay = N/8;
delay = N/16;
sample_in = 0;
K = N/2 - delay;
K = M-N - delay; //N/2 - delay; // N+K < M
LOG2N = 2 + (int)(log(N)/log(2));
LOG2N = 2 + (int)(log(N+K)/log(2));
N_DFT = 1 << LOG2N;
if (N + K > N_DFT/2 - 2) {
while (N + K > N_DFT/2 - 2) {
LOG2N += 1;
N_DFT <<= 1;
}
@ -475,15 +504,15 @@ int init_buffers(char hdr[], int hLen, int shape) {
xn = calloc(N_DFT+1, sizeof(float)); if (xn == NULL) return -1;
ew = calloc(LOG2N+1, sizeof(complex double)); if (ew == NULL) return -1;
Fm = calloc(N_DFT+1, sizeof(complex double)); if (Fm == NULL) return -1;
X = calloc(N_DFT+1, sizeof(complex double)); if (X == NULL) return -1;
Z = calloc(N_DFT+1, sizeof(complex double)); if (Z == NULL) return -1;
cx = calloc(N_DFT+1, sizeof(complex double)); if (cx == NULL) return -1;
ew = calloc(LOG2N+1, sizeof(complex float)); if (ew == NULL) return -1;
Fm = calloc(N_DFT+1, sizeof(complex float)); if (Fm == NULL) return -1;
X = calloc(N_DFT+1, sizeof(complex float)); if (X == NULL) return -1;
Z = calloc(N_DFT+1, sizeof(complex float)); if (Z == NULL) return -1;
cx = calloc(N_DFT+1, sizeof(complex float)); if (cx == NULL) return -1;
for (n = 0; n < LOG2N; n++) {
k = 1 << n;
ew[n] = cexp(-I*M_PI/(double)k);
ew[n] = cexp(-I*M_PI/(float)k);
}
m = calloc(N_DFT+1, sizeof(float)); if (m == NULL) return -1;

Wyświetl plik

@ -3,9 +3,10 @@ float read_wav_header(FILE*, float);
int f32buf_sample(FILE*, int, int);
int read_sbit(FILE*, int, int*, int, int, int, int);
int getCorrDFT(int, unsigned int, float *, unsigned int *);
int headcmp(int, char*, int, unsigned int);
int getCorrDFT(int, int, unsigned int, float *, unsigned int *);
int headcmp(int, char*, int, unsigned int, int, int);
float get_bufvar(int);
float get_bufmu(int);
int init_buffers(char*, int, int);
int free_buffers(void);

Wyświetl plik

@ -691,7 +691,7 @@ int main(int argc, char **argv) {
k += 1;
if (k >= K-4) {
mv0_pos = mv_pos;
mp = getCorrDFT(K, 0, &mv, &mv_pos);
mp = getCorrDFT(0, K, 0, &mv, &mv_pos);
k = 0;
}
else {
@ -703,10 +703,10 @@ int main(int argc, char **argv) {
if (mv_pos > mv0_pos) {
header_found = 0;
herrs = headcmp(symlen, rawheader, headerlen, mv_pos); // symlen=2
herrs = headcmp(symlen, rawheader, headerlen, mv_pos, mv<0, 0); // symlen=2
herr1 = 0;
if (herrs <= 3 && herrs > 0) {
herr1 = headcmp(symlen, rawheader, headerlen, mv_pos+1);
herr1 = headcmp(symlen, rawheader, headerlen, mv_pos+1, mv<0, 0);
if (herr1 < herrs) {
herrs = herr1;
herr1 = 1;

Wyświetl plik

@ -1172,7 +1172,7 @@ int main(int argc, char *argv[]) {
k += 1;
if (k >= K-4) {
mv0_pos = mv_pos;
mp = getCorrDFT(K, 0, &mv, &mv_pos);
mp = getCorrDFT(0, K, 0, &mv, &mv_pos);
k = 0;
}
else {
@ -1184,10 +1184,10 @@ int main(int argc, char *argv[]) {
if (mv_pos > mv0_pos) {
header_found = 0;
herrs = headcmp(symlen, header, headerlen, mv_pos); // symlen=1
herrs = headcmp(symlen, header, headerlen, mv_pos, mv<0, 0); // symlen=1
herr1 = 0;
if (herrs <= 3 && herrs > 0) {
herr1 = headcmp(symlen, header, headerlen, mv_pos+1);
herr1 = headcmp(symlen, header, headerlen, mv_pos+1, mv<0, 0);
if (herr1 < herrs) {
herrs = herr1;
herr1 = 1;

Wyświetl plik

@ -1404,7 +1404,7 @@ int main(int argc, char *argv[]) {
k += 1;
if (k >= K-4) {
mv0_pos = mv_pos;
mp = getCorrDFT(K, 0, &mv, &mv_pos);
mp = getCorrDFT(0, K, 0, &mv, &mv_pos);
k = 0;
}
else {
@ -1416,10 +1416,10 @@ int main(int argc, char *argv[]) {
if (mv_pos > mv0_pos) {
header_found = 0;
herrs = headcmp(symlen, rawheader, headerlen, mv_pos); // symlen=2
herrs = headcmp(symlen, rawheader, headerlen, mv_pos, mv<0, 0); // symlen=2
herr1 = 0;
if (herrs <= 3 && herrs > 0) {
herr1 = headcmp(symlen, rawheader, headerlen, mv_pos+1);
herr1 = headcmp(symlen, rawheader, headerlen, mv_pos+1, mv<0, 0);
if (herr1 < herrs) {
herrs = herr1;
herr1 = 1;