mirror
/
RS-tracker
kopia lustrzana https://github.com/rs1729/RS
1
0
Forkuj 0
Kod Zgłoszenia Projekty Wydania Wiki Aktywność

C50: speed-up

pull/53/head
Zilog80 2022-12-03 10:06:04 +01:00
rodzic 015137ccfc
commit 86f0996fab
1 zmienionych plików z 75 dodań i 155 usunięć
Pokaż statystyki Ściągnij plik aktualizacji Ściągnij plik porównania Expand all files Collapse all files

230
c34/c50dft.c
Wyświetl plik

@ -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,83 +528,71 @@ 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;
cb0 = dft_k(k0);
cb1 = dft_k(k1);
for (j = 0; j < N; j++) {
xn[j] = Hann[j]*buffer[(ptr + j + 1)%N];
}
bit = (cb0 > cb1) ? 0 : 1;
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
}
if (bit != bit0) {
if (bit != bit0) {
pos0 = pos;
pos = sample_count; //sample_count-(N-1)/2
pos0 = pos;
pos = sample_count; //sample_count-(N-1)/2
len = (pos-pos0+bitlen/2)/bitlen; //(pos-pos0)/(float)bitlen + 0.5;
for (i = 0; i < len; i++) {
inc_bufpos();
buf[bufpos] = 0x30 + bit0;
len = (pos-pos0+bitlen/2)/bitlen; //(pos-pos0)/bitlen + 0.5;
for (i = 0; i < len; i++) {
inc_bufpos();
buf[bufpos] = 0x30 + bit0;
if (!header_found) {
if (compare() >= HEADLEN) {
header_found = 1;
for (bitpos = 0; bitpos < HEADLEN; bitpos++) bits[bitpos] = header[bitpos] & 0x1;
}
}
else {
bits[bitpos] = bit0;
bitpos++;
if (bitpos >= LEN_BITFRAME) {
bits2bytes8N1(bits, bytes, bitpos/BITS);
if (option_raw) {
printRaw(bitpos/BITS);
}
else {
evalBytes2();
}
bitpos = 0;
header_found = 0;
}
if (!header_found) {
if (compare() >= HEADLEN-1) {
header_found = 1;
for (bitpos = 0; bitpos < HEADLEN; bitpos++) bits[bitpos] = header[bitpos] & 0x1;
}
}
else {
bits[bitpos] = bit0;
bitpos++;
if (bitpos >= LEN_BITFRAME) {
bits2bytes8N1(bits, bytes, bitpos/BITS);
if (option_raw) {
printRaw(bitpos/BITS);
}
else {
evalBytes2();
}
bitpos = 0;
header_found = 0;
}
}
bit0 = bit;
}
bit0 = bit;
}
}
printf("\n");