kopia lustrzana https://github.com/rs1729/RS
C50: speed-up
rodzic
f20f033472
commit
37ac80ce11
172
c34/c50dft.c
172
c34/c50dft.c
|
@ -1,7 +1,7 @@
|
|||
|
||||
/*
|
||||
C50
|
||||
(empfohlen: sample rate 48kHz)
|
||||
(recommended: sample rate 48kHz)
|
||||
gcc c50dft.c -lm -o c50dft
|
||||
*/
|
||||
|
||||
|
@ -32,7 +32,6 @@ typedef unsigned char ui8_t;
|
|||
static int option_verbose = 0,
|
||||
option_raw = 0,
|
||||
option_ptu = 0,
|
||||
option_dft = 0,
|
||||
option_json = 0,
|
||||
wavloaded = 0;
|
||||
|
||||
|
@ -73,11 +72,13 @@ static int read_wav_header(FILE *fp) {
|
|||
int byte, p=0;
|
||||
|
||||
if (fread(txt, 1, 4, fp) < 4) return -1;
|
||||
if (strncmp(txt, "RIFF", 4)) return -1;
|
||||
if (strncmp(txt, "RIFF", 4) && strncmp(txt, "RF64", 4)) return -1;
|
||||
|
||||
if (fread(txt, 1, 4, fp) < 4) return -1;
|
||||
// pos_WAVE = 8L
|
||||
if (fread(txt, 1, 4, fp) < 4) return -1;
|
||||
if (strncmp(txt, "WAVE", 4)) return -1;
|
||||
|
||||
// pos_fmt = 12L
|
||||
for ( ; ; ) {
|
||||
if ( (byte=fgetc(fp)) == EOF ) return -1;
|
||||
|
@ -115,7 +116,9 @@ static int read_wav_header(FILE *fp) {
|
|||
fprintf(stderr, "bits : %d\n", bits_sample);
|
||||
fprintf(stderr, "channels : %d\n", channels);
|
||||
|
||||
if ((bits_sample != 8) && (bits_sample != 16)) return -1;
|
||||
if ((bits_sample != 8) && (bits_sample != 16) && (bits_sample != 32)) return -1;
|
||||
|
||||
if (sample_rate == 900001) sample_rate -= 1;
|
||||
|
||||
//samples_per_bit = sample_rate/(float)BAUD_RATE;
|
||||
//fprintf(stderr, "samples/bit: %.2f\n", samples_per_bit);
|
||||
|
@ -123,30 +126,33 @@ static int read_wav_header(FILE *fp) {
|
|||
return 0;
|
||||
}
|
||||
|
||||
#define EOF_INT 0x1000000
|
||||
static unsigned int sample_count;
|
||||
|
||||
static int read_signed_sample(FILE *fp) { // int = i32_t
|
||||
int byte, i, ret; // EOF -> 0x1000000
|
||||
static int f32read_sample(FILE *fp, float *s) {
|
||||
int i;
|
||||
unsigned int word = 0;
|
||||
short *b = (short*)&word;
|
||||
float *f = (float*)&word;
|
||||
|
||||
for (i = 0; i < channels; i++) {
|
||||
// i = 0: links bzw. mono
|
||||
byte = fgetc(fp);
|
||||
if (byte == EOF) return EOF_INT;
|
||||
if (i == 0) ret = byte;
|
||||
|
||||
if (bits_sample == 16) {
|
||||
byte = fgetc(fp);
|
||||
if (byte == EOF) return EOF_INT;
|
||||
if (i == 0) ret += byte << 8;
|
||||
if (fread( &word, bits_sample/8, 1, fp) != 1) return EOF;
|
||||
|
||||
if (i == 0) { // i = 0: links bzw. mono
|
||||
//if (bits_sample == 8) sint = b-128; // 8bit: 00..FF, centerpoint 0x80=128
|
||||
//if (bits_sample == 16) sint = (short)b;
|
||||
|
||||
if (bits_sample == 32) {
|
||||
*s = *f;
|
||||
}
|
||||
else {
|
||||
if (bits_sample == 8) { *b -= 128; }
|
||||
*s = *b/(float)(1<<bits_sample);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
if (bits_sample == 8) return ret-128;
|
||||
if (bits_sample == 16) return (short)ret;
|
||||
|
||||
return ret;
|
||||
return 0;
|
||||
}
|
||||
|
||||
/* ------------------------------------------------------------------------------------ */
|
||||
|
@ -169,11 +175,10 @@ static ui8_t bits[LEN_BITFRAME+20] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,
|
|||
0, 1, 1, 1, 1, 1, 1, 1, 1, 1};
|
||||
static ui8_t bytes[LEN_BITFRAME/BITS+2];
|
||||
|
||||
static float x[N];
|
||||
static float complex Z[N], w[N], expw[N][N], ew[N*N];
|
||||
static float complex w[N], ew[N*N];
|
||||
|
||||
static int ptr;
|
||||
static float Hann[N], buffer[N+1], xn[N];
|
||||
static float Hann[N], buffer[N+1];
|
||||
|
||||
|
||||
static void init_dft() {
|
||||
|
@ -186,8 +191,7 @@ static void init_dft() {
|
|||
for (k = 0; k < N; k++) {
|
||||
w[k] = -I*2*M_PI * k / (float)N;
|
||||
for (n = 0; n < N; n++) {
|
||||
expw[k][n] = cexp( w[k] * n );
|
||||
ew[k*n] = expw[k][n];
|
||||
ew[k*n] = cexp( w[k] * n );
|
||||
}
|
||||
}
|
||||
}
|
||||
|
@ -198,78 +202,13 @@ static float dft_k(int k) {
|
|||
float complex Zk;
|
||||
|
||||
Zk = 0;
|
||||
for (n = 0; n < N; n++) {
|
||||
Zk += xn[n] * ew[k*n];
|
||||
for (n = 0; n < WLEN; n++) { //Hann[WLEN:N]=0
|
||||
// Hann[n]*buffer[(ptr + n + 1)%N] * ew[k*n];
|
||||
Zk += Hann[n]*buffer[(ptr + n + 1)&(N-1)] * ew[k*n]; // N=2^6=64
|
||||
}
|
||||
return cabs(Zk);
|
||||
}
|
||||
|
||||
static void dft() {
|
||||
int k, n;
|
||||
|
||||
for (k = 0; k < N/2; k++) { // xn reell, brauche nur N/2 unten
|
||||
Z[k] = 0;
|
||||
for (n = 0; n < N; n++) {
|
||||
Z[k] += xn[n] * ew[k*n];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
static void dft2() {
|
||||
int s, l, l2, i, j, k;
|
||||
float complex w1, w2, T;
|
||||
|
||||
for (i = 0; i < N; i++) {
|
||||
Z[i] = (float complex)xn[i];
|
||||
}
|
||||
|
||||
j = 1;
|
||||
for (i = 1; i < N; i++) {
|
||||
if (i < j) {
|
||||
T = Z[j-1];
|
||||
Z[j-1] = Z[i-1];
|
||||
Z[i-1] = T;
|
||||
}
|
||||
k = N/2;
|
||||
while (k < j) {
|
||||
j = j - k;
|
||||
k = k/2;
|
||||
}
|
||||
j = j + k;
|
||||
}
|
||||
|
||||
for (s = 0; s < LOG2N; s++) {
|
||||
l2 = 1 << s;
|
||||
l = l2 << 1;
|
||||
w1 = (float complex)1.0;
|
||||
w2 = cexp(-I*M_PI/(float)l2);
|
||||
for (j = 1; j <= l2; j++) {
|
||||
for (i = j; i <= N; i += l) {
|
||||
k = i + l2;
|
||||
T = Z[k-1] * w1;
|
||||
Z[k-1] = Z[i-1] - T;
|
||||
Z[i-1] = Z[i-1] + T;
|
||||
}
|
||||
w1 = w1 * w2;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
static int max_bin() {
|
||||
int k, kmax;
|
||||
float max;
|
||||
|
||||
max = 0; kmax = 0;
|
||||
for (k = 0; k < N/2-1; k++) {
|
||||
if (cabs(Z[k]) > max) {
|
||||
max = cabs(Z[k]);
|
||||
kmax = k;
|
||||
}
|
||||
}
|
||||
|
||||
return kmax;
|
||||
}
|
||||
|
||||
static float freq2bin(int f) {
|
||||
return f * N / (float)sample_rate;
|
||||
}
|
||||
|
@ -521,15 +460,15 @@ int main(int argc, char *argv[]) {
|
|||
|
||||
FILE *fp;
|
||||
char *fpname;
|
||||
int sample;
|
||||
int i, j, kmax, k0, k1;
|
||||
int i, k0, k1;
|
||||
int bit = 8, bit0 = 8;
|
||||
int pos = 0, pos0 = 0;
|
||||
int header_found = 0;
|
||||
int bitlen; // sample_rate/BAUD_RATE
|
||||
int len;
|
||||
float k_f0, k_f1, k_df;
|
||||
float k_f0, k_f1;
|
||||
float cb0, cb1;
|
||||
float s = 0.0;
|
||||
int cfreq = -1;
|
||||
|
||||
#ifdef CYGWIN
|
||||
|
@ -557,14 +496,7 @@ int main(int argc, char *argv[]) {
|
|||
else if ( (strcmp(*argv, "--ptu") == 0) ) {
|
||||
option_ptu = 1;
|
||||
}
|
||||
else if ( (strcmp(*argv, "-d1") == 0) || (strcmp(*argv, "--dft1") == 0) ) {
|
||||
option_dft = 1;
|
||||
}
|
||||
else if ( (strcmp(*argv, "-d2") == 0) || (strcmp(*argv, "--dft2") == 0) ) {
|
||||
option_dft = 2;
|
||||
}
|
||||
else if ( (strcmp(*argv, "--json") == 0) ) {
|
||||
option_dft = 0;
|
||||
option_verbose = 1;
|
||||
option_json = 1;
|
||||
}
|
||||
|
@ -578,7 +510,7 @@ int main(int argc, char *argv[]) {
|
|||
else {
|
||||
fp = fopen(*argv, "rb");
|
||||
if (fp == NULL) {
|
||||
fprintf(stderr, "%s konnte nicht geoeffnet werden\n", *argv);
|
||||
fprintf(stderr, "error: open %s\n", *argv);
|
||||
return -1;
|
||||
}
|
||||
wavloaded = 1;
|
||||
|
@ -596,58 +528,46 @@ int main(int argc, char *argv[]) {
|
|||
fclose(fp);
|
||||
return -1;
|
||||
}
|
||||
if ( sample_rate != 48000 ) {
|
||||
fprintf(stderr, "note: sample rate not 48000\n");
|
||||
}
|
||||
|
||||
|
||||
bitlen = sample_rate/BAUD_RATE;
|
||||
k_f0 = freq2bin(4700); // bit0: 4800Hz
|
||||
k_f1 = freq2bin(2900); // bit1: 3000Hz
|
||||
k_df = fabs(k_f0-k_f1)/2.5;
|
||||
k0 = (int)(k_f0+.5);
|
||||
k1 = (int)(k_f1+.5);
|
||||
|
||||
init_dft();
|
||||
|
||||
ptr = -1; sample_count = -1;
|
||||
while ((sample=read_signed_sample(fp)) < EOF_INT) {
|
||||
while (f32read_sample(fp, &s) != EOF) {
|
||||
|
||||
ptr++;
|
||||
sample_count++;
|
||||
if (ptr == N) ptr = 0;
|
||||
buffer[ptr] = sample / (float)(1<<bits_sample);
|
||||
buffer[ptr] = s;
|
||||
|
||||
if (sample_count < N) continue;
|
||||
|
||||
|
||||
for (j = 0; j < N; j++) {
|
||||
xn[j] = Hann[j]*buffer[(ptr + j + 1)%N];
|
||||
}
|
||||
|
||||
if (option_dft) {
|
||||
if (option_dft == 2) dft2();
|
||||
else dft();
|
||||
kmax = max_bin();
|
||||
if (kmax > k_f0-k_df && kmax < k_f0+k_df) bit = 0; // kmax = freq2bin(4800): 4800Hz
|
||||
else if (kmax > k_f1-k_df && kmax < k_f1+k_df) bit = 1; // kmax = freq2bin(3000): 3000Hz
|
||||
}
|
||||
else {
|
||||
cb0 = dft_k(k0);
|
||||
cb1 = dft_k(k1);
|
||||
if ( cb0 > cb1 ) bit = 0; // freq2bin(4800) : 4800Hz
|
||||
else bit = 1; // freq2bin(3000) : 3000Hz
|
||||
}
|
||||
|
||||
bit = (cb0 > cb1) ? 0 : 1;
|
||||
|
||||
if (bit != bit0) {
|
||||
|
||||
pos0 = pos;
|
||||
pos = sample_count; //sample_count-(N-1)/2
|
||||
|
||||
len = (pos-pos0+bitlen/2)/bitlen; //(pos-pos0)/bitlen + 0.5;
|
||||
len = (pos-pos0+bitlen/2)/bitlen; //(pos-pos0)/(float)bitlen + 0.5;
|
||||
for (i = 0; i < len; i++) {
|
||||
inc_bufpos();
|
||||
buf[bufpos] = 0x30 + bit0;
|
||||
|
||||
if (!header_found) {
|
||||
if (compare() >= HEADLEN) {
|
||||
if (compare() >= HEADLEN-1) {
|
||||
header_found = 1;
|
||||
for (bitpos = 0; bitpos < HEADLEN; bitpos++) bits[bitpos] = header[bitpos] & 0x1;
|
||||
}
|
||||
|
|
Ładowanie…
Reference in New Issue