Merge branch 'rs1729:master' into master

pull/43/head
Pavol Gajdoš 2021-07-02 17:53:31 +02:00 zatwierdzone przez GitHub
commit 53790c6b74
Nie znaleziono w bazie danych klucza dla tego podpisu
ID klucza GPG: 4AEE18F83AFDEB23
20 zmienionych plików z 1727 dodań i 437 usunięć

Wyświetl plik

@ -2,40 +2,59 @@
/*
C34
(empfohlen: sample rate 48kHz)
gcc c34dft.c -lm -o c34dft
*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <complex.h>
#include <math.h>
// optional JSON "version"
// (a) set global
// gcc -DVERSION_JSN [-I<inc_dir>] ...
#ifdef VERSION_JSN
#include "version_jsn.h"
#endif
// or
// (b) set local compiler option, e.g.
// gcc -DVER_JSN_STR=\"0.0.2\" ...
typedef unsigned char ui8_t;
int option_verbose = 0,
option_raw = 0,
option_dft = 0,
wavloaded = 0;
static int option_verbose = 0,
option_raw = 0,
option_ptu = 0,
option_dft = 0,
option_json = 0,
wavloaded = 0;
typedef struct {
int frnr;
//int frnr;
int sn;
int jahr; int monat; int tag;
int std; int min; int sek;
double lat; double lon; double h;
float lat; float lon; float alt;
unsigned chk;
float T; float RH;
int jsn_freq; // freq/kHz (SDR)
} gpx_t;
gpx_t gpx;
static gpx_t gpx;
/* ------------------------------------------------------------------------------------ */
#define BAUD_RATE 2400
int sample_rate = 0, bits_sample = 0, channels = 0;
static unsigned int sample_rate = 0;
static int bits_sample = 0, channels = 0;
//float samples_per_bit = 0;
int findstr(char *buff, char *str, int pos) {
static int findstr(char *buff, char *str, int pos) {
int i;
for (i = 0; i < 4; i++) {
if (buff[(pos+i)%4] != str[i]) break;
@ -43,7 +62,7 @@ int findstr(char *buff, char *str, int pos) {
return i;
}
int read_wav_header(FILE *fp) {
static int read_wav_header(FILE *fp) {
char txt[4+1] = "\0\0\0\0";
unsigned char dat[4];
int byte, p=0;
@ -100,8 +119,9 @@ int read_wav_header(FILE *fp) {
}
#define EOF_INT 0x1000000
static unsigned int sample_count;
int read_signed_sample(FILE *fp) { // int = i32_t
static int read_signed_sample(FILE *fp) { // int = i32_t
int byte, i, ret; // EOF -> 0x1000000
for (i = 0; i < channels; i++) {
@ -135,31 +155,31 @@ int read_signed_sample(FILE *fp) { // int = i32_t
#define LEN_BITFRAME 84 // 7*(4+8)
#define HEADLEN 28
char header[] = "1110000000001110111111111110";
char buf[HEADLEN+1] = "x";
int bufpos = -1;
char headerstr[] = "1110 00000000 1110 11111111";
static char header[] = "1110000000001110111111111110";
static char buf[HEADLEN+1] = "x";
static int bufpos = -1;
static char headerstr[] = "1110 00000000 1110 11111111";
int bitpos;
ui8_t bits[LEN_BITFRAME+1] = { 1, 1, 1, 0};
ui8_t bytes[LEN_BITFRAME/BITS];
static int bitpos;
static ui8_t bits[LEN_BITFRAME+1] = { 1, 1, 1, 0};
static ui8_t bytes[LEN_BITFRAME/BITS];
double x[N];
double complex Z[N], w[N], expw[N][N], ew[N*N];
static float x[N];
static float complex Z[N], w[N], expw[N][N], ew[N*N];
int ptr;
double Hann[N], buffer[N+1], xn[N];
static int ptr;
static float Hann[N], buffer[N+1], xn[N];
void init_dft() {
static void init_dft() {
int i, k, n;
for (i = 0; i < N; i++) Hann[i] = 0;
for (i = 0; i < WLEN; i++) Hann[i] = 0.5 * (1 - cos( 2 * M_PI * i / (double)(WLEN-1) ) );
//Hann[i+(N-1-WLEN)/2] = 0.5 * (1 - cos( 2 * M_PI * i / (double)(WLEN-1) ) );
for (i = 0; i < WLEN; i++) Hann[i] = 0.5 * (1 - cos( 2 * M_PI * i / (float)(WLEN-1) ) );
//Hann[i+(N-1-WLEN)/2] = 0.5 * (1 - cos( 2 * M_PI * i / (float)(WLEN-1) ) );
for (k = 0; k < N; k++) {
w[k] = -I*2*M_PI * k / (double)N;
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];
@ -168,9 +188,9 @@ void init_dft() {
}
double dft_k(int k) {
static float dft_k(int k) {
int n;
double complex Zk;
float complex Zk;
Zk = 0;
for (n = 0; n < N; n++) {
@ -179,7 +199,7 @@ double dft_k(int k) {
return cabs(Zk);
}
void dft() {
static void dft() {
int k, n;
for (k = 0; k < N/2; k++) { // xn reell, brauche nur N/2 unten
@ -190,12 +210,12 @@ void dft() {
}
}
void dft2() {
static void dft2() {
int s, l, l2, i, j, k;
double complex w1, w2, T;
float complex w1, w2, T;
for (i = 0; i < N; i++) {
Z[i] = (double complex)xn[i];
Z[i] = (float complex)xn[i];
}
j = 1;
@ -216,8 +236,8 @@ void dft2() {
for (s = 0; s < LOG2N; s++) {
l2 = 1 << s;
l = l2 << 1;
w1 = (double complex)1.0;
w2 = cexp(-I*M_PI/(double)l2);
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;
@ -230,9 +250,9 @@ void dft2() {
}
}
int max_bin() {
static int max_bin() {
int k, kmax;
double max;
float max;
max = 0; kmax = 0;
for (k = 0; k < N/2-1; k++) {
@ -245,22 +265,22 @@ int max_bin() {
return kmax;
}
double freq2bin(int f) {
return f * N / (double)sample_rate;
static float freq2bin(int f) {
return f * N / (float)sample_rate;
}
int bin2freq(int k) {
static int bin2freq(int k) {
return sample_rate * k / N;
}
/* ------------------------------------------------------------------------------------ */
void inc_bufpos() {
static void inc_bufpos() {
bufpos = (bufpos+1) % HEADLEN;
}
int compare() {
static int compare() {
int i=0, j = bufpos;
while (i < HEADLEN) {
@ -272,7 +292,7 @@ int compare() {
return i;
}
int bits2bytes(ui8_t bits[], ui8_t bytes[]) {
static int bits2bytes(ui8_t bits[], ui8_t bytes[]) {
int i, j, byteval=0, d=1;
for (j = 0; j < 7; j++) {
@ -288,26 +308,60 @@ int bits2bytes(ui8_t bits[], ui8_t bytes[]) {
return 0;
}
void printGPX() {
static void printGPX() {
int i;
printf(" %4d-%02d-%02d", gpx.jahr, gpx.monat, gpx.tag);
printf(" %2d:%02d:%02d", gpx.std, gpx.min, gpx.sek);
if (gpx.sn) printf("( %d ) ", gpx.sn);
printf(" %04d-%02d-%02d", gpx.jahr, gpx.monat, gpx.tag);
printf(" %02d:%02d:%02d", gpx.std, gpx.min, gpx.sek);
printf(" ");
printf(" lat: %.5f°", gpx.lat);
printf(" lon: %.5f°", gpx.lon);
printf(" h: %.1fm", gpx.h);
printf(" lat: %.5f", gpx.lat);
printf(" lon: %.5f", gpx.lon);
printf(" alt: %.1f", gpx.alt);
if (option_ptu && (gpx.T > -273.0 || gpx.RH > -0.5)) {
printf(" ");
if (gpx.T > -273.0) printf(" T=%.1fC", gpx.T);
//if (gpx.RH > -0.5) printf(" RH=%.0f%%", gpx.RH);
}
if (option_verbose) {
printf(" # ");
for (i = 0; i < 5; i++) printf("%d", (gpx.chk>>i)&1);
if (option_ptu) for (i = 6; i < 7; i++) printf("%d", (gpx.chk>>i)&1);
}
printf("\n");
}
static void printJSON() {
// UTC or GPS time ?
char *ver_jsn = NULL;
char json_sonde_id[] = "C34-xxxx\0\0\0\0\0\0\0";
if (gpx.sn) {
sprintf(json_sonde_id, "C34-%u", gpx.sn);
}
printf("{ \"type\": \"%s\"", "C34");
printf(", \"id\": \"%s\", ", json_sonde_id);
printf("\"datetime\": \"%04d-%02d-%02dT%02d:%02d:%02dZ\", \"lat\": %.5f, \"lon\": %.5f, \"alt\": %.1f",
gpx.jahr, gpx.monat, gpx.tag, gpx.std, gpx.min, gpx.sek, gpx.lat, gpx.lon, gpx.alt);
if (option_ptu && (gpx.T > -273.0 || gpx.RH > -0.5)) {
if (gpx.T > -273.0) printf(", \"temp\": %.1f", gpx.T);
//if (gpx.RH > -0.5) printf(", \"humidity\": %.1f", gpx.RH);
}
if (gpx.jsn_freq > 0) {
printf(", \"freq\": %d", gpx.jsn_freq);
}
#ifdef VER_JSN_STR
ver_jsn = VER_JSN_STR;
#endif
if (ver_jsn && *ver_jsn != '\0') printf(", \"version\": \"%s\"", ver_jsn);
printf(" }\n");
//printf("\n");
}
// Chechsum Fletcher16
unsigned check2(ui8_t *bytes, int len) {
static unsigned check2(ui8_t *bytes, int len) {
int sum1, sum2;
int i;
@ -323,7 +377,7 @@ unsigned check2(ui8_t *bytes, int len) {
return sum2 | (sum1<<8);
}
/* // equivalent
unsigned check16(ui8_t *bytes, int len) {
static unsigned check16(ui8_t *bytes, int len) {
unsigned sum1, sum2;
int i;
sum1 = sum2 = 0;
@ -336,16 +390,20 @@ unsigned check16(ui8_t *bytes, int len) {
}
*/
double NMEAll(int ll) { // NMEA GGA,GLL: ll/1e4=(D)DDMM.mmmm
static float NMEAll(int ll) { // NMEA GGA,GLL: ll/1e4=(D)DDMM.mmmm
int deg = ll / 1000000;
double min = (ll - deg*1000000)/1e4;
float min = (ll - deg*1000000)/1e4;
return deg+min/60.0;
}
int evalBytes() {
static int evalBytes() {
int i, val = 0;
ui8_t id = bytes[0];
unsigned check;
static unsigned int cnt_dat = -1, cnt_tim = -1,
cnt_lat = -1, cnt_lon = -1, cnt_alt = -1,
cnt_sn = -1,
cnt_t3 = -1, cnt_rh = -1;
check = ((bytes[5]<<8)|bytes[6]) != check2(bytes, 5);
@ -353,44 +411,79 @@ int evalBytes() {
if (id == 0x14 ) { // date
int tag = val / 10000;
int mon = (val-tag*10000) / 100;
int jrz = val % 100;
int mon = (val-tag*10000) / 100;
int jrz = val % 100;
gpx.tag = tag;
gpx.monat = mon;
gpx.jahr = 2000+jrz;
gpx.chk = (gpx.chk & ~(0x1<<0)) | (check<<0);
if (check==0) cnt_dat = sample_count;
}
else if (id == 0x15 ) { // time
int std = val / 10000;
int min = (val-std*10000) / 100;
int sek = val % 100;
int min = (val-std*10000) / 100;
int sek = val % 100;
gpx.std = std;
gpx.min = min;
gpx.sek = sek;
gpx.chk = (gpx.chk & ~(0x1<<1)) | (check<<1);
if (check==0) cnt_tim = sample_count;
}
else if (id == 0x16 ) { // lat: wie NMEA mit Faktor 1e4
gpx.lat = NMEAll(val);
gpx.chk = (gpx.chk & ~(0x1<<2)) | (check<<2);
if (check==0) cnt_lat = sample_count;
}
else if (id == 0x17 ) { // lon: wie NMEA mit Faktor 1e4
gpx.lon = NMEAll(val);
gpx.chk = (gpx.chk & ~(0x1<<3)) | (check<<3);
if (check==0) cnt_lon = sample_count;
}
else if (id == 0x18 ) { // alt: decimeter
gpx.h = val/10.0;
gpx.alt = val/10.0;
gpx.chk = (gpx.chk & ~(0x1<<4)) | (check<<4);
if (check==0) cnt_alt = sample_count;
}
else if (id == 0x64 ) { // serial number
if (check==0) gpx.sn = val; // 16 bit
//gpx.chk = (gpx.chk & ~(0x1<<15)) | (check<<15);
//if (check==0) cnt_sn = sample_count;
}
if (id == 0x18) {
printGPX();
if (option_json && check==0) {
if ( ((cnt_dat|cnt_tim|cnt_lat|cnt_lon)&0x80000000)==0 &&
cnt_alt - cnt_dat < sample_rate &&
cnt_alt - cnt_tim < sample_rate &&
cnt_alt - cnt_lat < sample_rate &&
cnt_alt - cnt_lon < sample_rate )
{
if (cnt_alt - cnt_t3 > sample_rate) gpx.T = -273.15;
if (cnt_alt - cnt_rh > sample_rate) gpx.RH = -1.0;
printJSON();
}
}
}
if (id == 0x18 && !option_raw) printGPX();
// PTU floats
if (id == 0x03) { // temperature
float t = -273.15;
memcpy(&t, &val, 4);
if (t < -273.0 || t > 100.0) t = -273.15;
gpx.T = t;
gpx.chk = (gpx.chk & ~(0x1<<6)) | (check<<6);
if (check==0) cnt_t3 = sample_count;
}
return check;
}
void printRaw() {
static void printRaw() {
int j;
//if ( ((bytes[5]<<8)|bytes[6]) == check2(bytes, 5))
unsigned chkbyt = (bytes[5]<<8) | bytes[6];
unsigned chksum = check2(bytes, 5);
//if (chksum == chkbyt)
{
printf("%s", headerstr);
for (j = 0; j < LEN_BITFRAME; j++) {
@ -402,9 +495,10 @@ void printRaw() {
printf("%02X%02X ", 0x00, 0xFF);
printf("%02X ", bytes[0]);
printf("%02X%02X%02X%02X ", bytes[1], bytes[2], bytes[3], bytes[4]);
printf("%02X%02X", bytes[5], bytes[6]);
printf("%02X%02X", bytes[5], bytes[6]); // chkbyt
if (option_verbose) {
printf(" # %04X", check2(bytes, 5));
printf(" # %04X", chksum);
if (chksum == chkbyt) printf(" [OK]"); else printf(" [NO]");
}
printf("\n");
}
@ -415,16 +509,16 @@ int main(int argc, char *argv[]) {
FILE *fp;
char *fpname;
int sample;
unsigned int sample_count;
int sample;
int i, j, kmax, k0, k1;
int bit = 8, bit0 = 8;
int pos = 0, pos0 = 0;
int header_found = 0;
int bitlen; // sample_rate/BAUD_RATE
int len;
double k_f0, k_f1, k_df;
double cb0, cb1;
float k_f0, k_f1, k_df;
float cb0, cb1;
int cfreq = -1;
fpname = argv[0];
++argv;
@ -442,12 +536,27 @@ int main(int argc, char *argv[]) {
else if ( (strcmp(*argv, "-r") == 0) || (strcmp(*argv, "--raw") == 0) ) {
option_raw = 1;
}
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;
}
else if ( (strcmp(*argv, "--jsn_cfq") == 0) ) {
int frq = -1; // center frequency / Hz
++argv;
if (*argv) frq = atoi(*argv); else return -1;
if (frq < 300000000) frq = -1;
cfreq = frq;
}
else {
fp = fopen(*argv, "rb");
if (fp == NULL) {
@ -461,6 +570,9 @@ int main(int argc, char *argv[]) {
if (!wavloaded) fp = stdin;
gpx.jsn_freq = 0;
if (cfreq > 0) gpx.jsn_freq = (cfreq+500)/1000;
i = read_wav_header(fp);
if (i) {
fclose(fp);
@ -483,7 +595,7 @@ int main(int argc, char *argv[]) {
ptr++;
sample_count++;
if (ptr == N) ptr = 0;
buffer[ptr] = sample / (double)(1<<bits_sample);
buffer[ptr] = sample / (float)(1<<bits_sample);
if (sample_count < N) continue;

Wyświetl plik

@ -23,16 +23,23 @@ receive IF stream from baseband IQ via TCP, default `PORT=1280 (iq_svcl.h)`<br /
&nbsp;&nbsp;&nbsp;&nbsp; `-0.5 < fq < 0.5`: (relative) frequency, `fq=frequency/sr` <br />
&nbsp;&nbsp;&nbsp;&nbsp; `<if_sr>`: IF sample rate <br />
&nbsp;&nbsp;&nbsp;&nbsp; `<bo>=8,16,32`: output/IF bits per (real) sample (u8, s16 or f32) <br />
down-converts up to `MAX_FQ=(4+1) (iq_base.h)` channels/signals. More signals than number of CPUs/cores is not recommended.<br />
down-converts up to `MAX_FQ=(4+1) (iq_base.h)` channels/signals. (On older CPUs, more signals than number of CPU cores is not recommended.)<br />
(Note: If the baseband sample rate has no appropriate factors (e.g. if prime), the IF sample rate might be high and IF-processing slow.)<br />
One channel can be used for scanning, `--fft <fft.txt>` makes FFT (2 seconds average).
The FFT is saved in `<fft.txt>` as `<fq>;<dB>`, approx. 200 Hz per bin.<br />
One channel can be used for scanning, e.g. `./iq_server --fft_avg <m> <fft_avg.csv>` makes *m* rows of avg-FFT (FFT_AVG=2 seconds average), and
this channel will be reused for client FFT requests. Only one channel/thread can be used for FFT/scanning. A client can request a new FFT,
if the last FFT has finished.<br />
There are two kinds of FFTs, `fft_all` and `fft_avg`. `fft_avg` integrates over FFT_AVG=2 seconds and can be used for signal peak scanning.
For waterfall display, `--fft_all <m> <fft_all.csv>` produces *m*\*FFT_AVG\*FFT_FPS/2 rows of FFT (FFT_AVG seconds, FFT_FPS/2 per sec).
The FFT is saved in `<fft.csv>` as<br />
&nbsp;&nbsp; `sec.ms,freq_min,freq_max,Hz/bin,N_bins, db_1,...,db_N`<br />
approx. 200 Hz per bin.
Choose `filename="-"` for `stdout`.<br />
If no output bps is chosen (`--bo [8,16,32]`), the IF bps is equal to the baseband bps. It is recommended to use
`--bo 32` (i.e. float32) output, then no quantization noise is introduced when converting from internal float32 samples.<br />
- Ex.2<br />
[terminal 1]<br />
`T1$ rtl_sdr -f 403.0M -s 1920000 - | ./iq_server --fft fft_server.txt --bo 32 - 1920000 8`<br />
`T1$ rtl_sdr -f 403.0M -s 1920000 - | ./iq_server --fft_avg 1 fft_server.csv --bo 32 - 1920000 8`<br />
[terminal 2]<br />
`T2$ ./iq_client --freq -0.3125 | ./m10mod -c -vv --IQ 0.0 - 48000 32`<br />
[terminal 3]<br />
@ -41,7 +48,17 @@ receive IF stream from baseband IQ via TCP, default `PORT=1280 (iq_svcl.h)`<br /
`T4$ ./iq_client -1` &nbsp;&nbsp; (*close channel 1*)<br />
`T4$ ./iq_client --stop` &nbsp;&nbsp; (*close all clients and stop server*)<br />
The `iq_server` `--fft` option immediately starts reading the IQ stream (so buffering is reduced).<br />
`./iq_client --fft <fft_cl.txt>` can also request FFT.<br />
The IF sample rate `if_sr` is at least 48000 such that the baseband sample rate `sr` is a multiple of `if_sr`.
- The `iq_server` FFT options immediately start reading the IQ stream (so buffering is reduced).<br />
`./iq_client <fft_opt> <m> <fft_cl.txt>`, where `<fft_opt>=--fft_avg_cl` or `--fft_all_cl`, requests FFT from the server.<br />
(`<fft_opt>=--fft_avg_sv/--fft_all_sv` would save the FFT at the server, but only if `./iq_server --enable_clsv_out`.)<br />
The IF sample rate `if_sr` is at least 48000 and such that the baseband sample rate `sr` is a multiple of `if_sr`.
- Ex.3<br />
[terminal 1]<br />
`T1$ rtl_sdr -f 404550k -s 2048000 - | ./iq_server --fft_avg 1 fft_avg.csv --bo 32 - 2048000 8`<br />
(scan FFT: `./scan_fft fft_avg.csv`)<br />
[terminal 3]<br />
`T3$ ./iq_client --fft_all_cl -1 - | python plot_fft_ani.py 3 -`<br />
![FFT image](fft3-1.png "FFT")

Plik binarny nie jest wyświetlany.

Po

Szerokość:  |  Wysokość:  |  Rozmiar: 212 KiB

Wyświetl plik

@ -34,7 +34,10 @@ typedef struct {
double xlt_fq;
float complex *blk;
int used;
//
int fft;
int stop_fft;
int fft_num;
} thd_t;
@ -97,6 +100,7 @@ typedef struct {
thd_t thd;
int fd;
char *fname;
FILE *fpo;
} thargs_t;

Wyświetl plik

@ -4,6 +4,27 @@
*
* gcc -O2 iq_client.c -o iq_client
*
* usage:
* (request IF IQ samples)
* ./iq_client [--ip <ip_adr>] [--port <pn>] --freq <fq> # -0.5 < fq < 0.5
*
* (request FFT)
* ./iq_client <fft_opt> <m> <filename> # FFT csv output
* <fft_opt>:
* --fft_avg_cl # client out
* --fft_all_cl # client out
* --fft_avg_sv # server out (if iq_server --enable_clsv_out)
* --fft_avg_sv # server out (if iq_server --enable_clsv_out)
* <m>:
* _avg_: m avg-FFTs
* _all_: m*FFT_FPS/2 FFTs
* m = -1 : continuous FFT output
* <filename>:
* csv filename ("-": stdout)
*
* ./iq_client -<n> # close client <n>
* ./iq_client --stop # close all clients, stop server
*
* author: zilog80
*/
@ -49,17 +70,24 @@ int main(int argc, char *argv[]) {
}
else serv_port = port;
}
else if (strcmp(*argv, "--fft0") == 0) {
sprintf(sendln, "%s", "--fft0");
}
else if (strcmp(*argv, "--fft") == 0) {
sprintf(sendln, "%s", "--fft");
++argv;
if (*argv) {
fname_fft = *argv;
else if (strncmp(*argv, "--fft", 5) == 0) {
char *arg_fft = *argv;
int fft_num = 0;
int opt_fft_cl = 0;
if (strncmp(arg_fft+5, "_avg", 4) == 0) opt_fft_cl = OPT_FFT_AVG;
else if (strncmp(arg_fft+5, "_all", 4) != 0) return -1;
if (strncmp(arg_fft+5+4, "_cl", 3) == 0) {
opt_fft_cl |= OPT_FFT_CLNT;
re = 2;
}
else if (strncmp(arg_fft+5+4, "_sv", 3) == 0) opt_fft_cl |= OPT_FFT_SERV;
else return -1;
re = 2;
++argv;
if (*argv) fft_num = atoi(*argv); else return -1;
++argv;
if (*argv) fname_fft = *argv; else return -1;
//sprintf(sendln, "%s_%s_%s", "--fft", (opt_fft_cl & OPT_FFT_AVG) ? "avg" : "all", (opt_fft_cl & OPT_FFT_CLNT) ? "cl" : "sv");
sprintf(sendln, "%s %d %s", arg_fft, fft_num, fname_fft);
}
else if (strcmp(*argv, "--freq") == 0) {
++argv;
@ -155,7 +183,9 @@ int main(int argc, char *argv[]) {
else if ( re == 2 )
{
// fft data
FILE *fpo = fopen(fname_fft, "wb");
FILE *fpo = NULL;
if (fname_fft[0] == '-') fpo = stdout;
else fpo = fopen(fname_fft, "wb");
if (fpo != NULL) {
memset(recvln, 0, LINELEN+1);

Wyświetl plik

@ -180,7 +180,7 @@ static float write_wav_header(pcm_t *pcm) {
}
fwrite("data", 1, 4, fp);
data = 0; // datasize
data = 0xFFFFFFFF; // datasize unknown
fwrite(&data, 1, 4, fp);
return 0;

Wyświetl plik

@ -7,6 +7,24 @@
*
* (gcc -O2 iq_client.c -o iq_client)
*
* usage:
* ./iq_server [--enable_clsv_out] [--port <pn>] [iq_baseband.wav]
* no wav-file: stdin
*
* ./iq_server [--bo <b>] - <sr> <bs> [iq_baseband.raw]
* <b>=8,16,32 bit client/IF output
*
* ./iq_server <fft_opt> <m> <filename> [iq_baseband.wav] # FFT csv output
* <fft_opt>:
* --fft_avg
* --fft_all
* <m>:
* _avg: m avg-FFTs
* _all: m*FFT_FPS/2 FFTs
* m = -1 : continuous FFT output
* <filename>:
* csv filename ("-": stdout)
*
* author: zilog80
*/
@ -27,13 +45,13 @@
#include "iq_svcl.h"
#include "iq_base.h"
#define FFT_AVG 2 // fft_avg: integrate FFT_AVG seconds, 2*FFT_FPS FFTs
#define FFT_FPS 16 // fft_all: output (ca.) FFT_FPS/2 per sec
#define FPOUT stderr
#define OPT_FFT_SERV 1 // server
#define OPT_FFT_CLSV 2 // server (client request)
#define OPT_FFT_CLNT 3 // server -> client
static int option_dbg = 0;
static int option_clsv_out = 0;
static int tcp_eof = 0;
@ -275,21 +293,82 @@ static void *thd_IF(void *targs) { // pcm_t *pcm, double xlt_fq
return NULL;
}
#define FFT_SEC 2
#define FFT_FPS 20
static int fft_txt_prn(FILE *fpo, dft_t *dft, float *db) {
int j;
fprintf(fpo, "# <freq/sr>;<dB> ## sr:%d , N:%d\n", dft->sr, dft->N);
for (j = dft->N/2; j < dft->N/2 + dft->N; j++) {
fprintf(fpo, "%+11.8f;%7.2f\n", bin2fq(dft, j % dft->N), db[j % dft->N]);
}
return 0;
}
static int fft_txt_tcp(int fd, dft_t *dft, float *db) {
char sendln[LINELEN+1];
int sendln_len;
int j, l;
snprintf(sendln, LINELEN, "# <freq/sr>;<dB> ## sr:%d , N:%d\n", dft->sr, dft->N);
sendln_len = strlen(sendln);
l = write(fd, sendln, sendln_len);
for (j = dft->N/2; j < dft->N/2 + dft->N; j++) {
memset(sendln, 0, LINELEN+1);
snprintf(sendln, LINELEN, "%+11.8f;%7.2f\n", bin2fq(dft, j % dft->N), db[j % dft->N]);
sendln_len = strlen(sendln);
l = write(fd, sendln, sendln_len);
}
return 0;
}
static int fft_csv_prn(FILE *fpo, dft_t *dft, float *db, double t_sec) {
int j;
fprintf(fpo, "%7.3f, ", t_sec);
fprintf(fpo, "%d, %d, ", (int)bin2freq(dft, dft->N/2), (int)bin2freq(dft, dft->N/2 - 1));
fprintf(fpo, "%.2f, ", dft->sr/(double)dft->N);
fprintf(fpo, "%d, ", dft->N);
for (j = dft->N/2; j < dft->N/2 + dft->N; j++) {
fprintf(fpo, "%7.2f%c", db[j % dft->N], j < dft->N/2 + dft->N-1 ? ',' : '\n');
}
return 0;
}
static int fft_csv_tcp(int fd, dft_t *dft, float *db, double t_sec) {
char sendln[LINELEN+1];
int sendln_len;
int j, l;
snprintf(sendln, LINELEN, "%7.3f, %d, %d, %.2f, %d, ",
t_sec, (int)bin2freq(dft, dft->N/2), (int)bin2freq(dft, dft->N/2 - 1), dft->sr/(double)dft->N, dft->N);
sendln_len = strlen(sendln);
l = write(fd, sendln, sendln_len);
for (j = dft->N/2; j < dft->N/2 + dft->N; j++) {
memset(sendln, 0, LINELEN+1);
snprintf(sendln, LINELEN, "%7.2f%c", db[j % dft->N], j < dft->N/2 + dft->N-1 ? ',' : '\n');
sendln_len = strlen(sendln);
l = write(fd, sendln, sendln_len);
}
return 0;
}
static void *thd_FFT(void *targs) {
thargs_t *tharg = targs;
pcm_t *pcm = &(tharg->pcm);
FILE *fpo = NULL;
char *fname_fft = "db_fft.txt";
int k;
int bitQ = 0;
float complex *z = NULL;
float *db = NULL;
float *all_rZ = NULL;
float *avg_rZ = NULL;
float *avg_db = NULL;
@ -312,7 +391,6 @@ static void *thd_FFT(void *targs) {
dsp.bps_out = pcm->bps_out;
//(dsp.thd)->fft = 1;
if (option_dbg) {
fprintf(stderr, "init FFT buffers\n");
}
@ -326,6 +404,8 @@ static void *thd_FFT(void *targs) {
z = calloc(dsp.decM+1, sizeof(float complex)); if (z == NULL) goto exit_thread;
db = calloc(dsp.DFT.N+1, sizeof(float)); if (db == NULL) goto exit_thread;
all_rZ = calloc(dsp.DFT.N+1, sizeof(float)); if (all_rZ == NULL) goto exit_thread;
avg_rZ = calloc(dsp.DFT.N+1, sizeof(float)); if (avg_rZ == NULL) goto exit_thread;
avg_db = calloc(dsp.DFT.N+1, sizeof(float)); if (avg_db == NULL) goto exit_thread;
@ -334,12 +414,16 @@ static void *thd_FFT(void *targs) {
int len = dsp.DFT.N / dsp.decM;
int mlen = len*dsp.decM;
int sum_n = 0;
int sec = FFT_SEC;
int fft_step = dsp.sr_base/(dsp.DFT.N*FFT_FPS);
int sum_fft = 0;
int avg_sec = FFT_AVG;
int fft_step = (int)(dsp.sr_base/(double)(dsp.DFT.N*FFT_FPS) + 0.5);
int n_fft = 0;
int th_used = 0;
int readSamples = 1;
int n_out = 0;
bitQ = 0;
while ( bitQ != EOF )
{
@ -367,65 +451,95 @@ static void *thd_FFT(void *targs) {
n++;
if (n == len) { // mlen = len * decM <= DFT.N
n_fft += 1;
if ((dsp.thd)->fft && sum_n*n_fft*mlen < sec*dsp.sr_base && n_fft >= fft_step)
if ( (dsp.thd)->fft )
{
for (j = 0; j < mlen; j++) {
dsp.DFT.Z[j] *= dsp.DFT.win[j];
if ((dsp.thd)->fft_num == 0)
{
if ( tharg->fpo ) { fclose(tharg->fpo); tharg->fpo = NULL; }
else if ( tharg->fd > STDIN_FILENO ) { close(tharg->fd); tharg->fd = -1; }
(dsp.thd)->fft = 0;
}
while (j < dsp.DFT.N) dsp.DFT.Z[j++] = 0.0; // dft(Z[...]) != 0
raw_dft(&(dsp.DFT), dsp.DFT.Z);
n_fft += 1;
for (j = 0; j < dsp.DFT.N; j++) avg_rZ[j] += cabs(dsp.DFT.Z[j]);
if (sum_n*n_fft*mlen < avg_sec*dsp.sr_base && n_fft >= fft_step) {
n_fft = fft_step;
sum_n++;
n_fft = 0;
}
if (sum_n*fft_step*mlen >= sec*dsp.sr_base) {
for (j = 0; j < dsp.DFT.N; j++) avg_rZ[j] /= dsp.DFT.N*(float)sum_n;
for (j = 0; j < dsp.DFT.N; j++) avg_db[j] = 20.0*log10(avg_rZ[j]+1e-20);
pthread_mutex_lock( (dsp.thd)->mutex );
fprintf(FPOUT, "<%d: FFT>\n", (dsp.thd)->tn);
pthread_mutex_unlock( (dsp.thd)->mutex );
if ( (dsp.thd)->fft == OPT_FFT_CLNT ) { // send FFT data to client
char sendln[LINELEN+1];
int sendln_len;
int l;
snprintf(sendln, LINELEN, "# <freq/sr>;<dB> ## sr:%d , N:%d\n", dsp.DFT.sr, dsp.DFT.N);
sendln_len = strlen(sendln);
l = write(tharg->fd, sendln, sendln_len);
for (j = dsp.DFT.N/2; j < dsp.DFT.N/2 + dsp.DFT.N; j++) {
memset(sendln, 0, LINELEN+1);
snprintf(sendln, LINELEN, "%+11.8f;%7.2f\n", bin2fq(&(dsp.DFT), j % dsp.DFT.N), avg_db[j % dsp.DFT.N]);
sendln_len = strlen(sendln);
l = write(tharg->fd, sendln, sendln_len);
if (sum_fft == 0) {
pthread_mutex_lock( (dsp.thd)->mutex );
fprintf(FPOUT, "<%d: FFT_START>\n", (dsp.thd)->tn);
pthread_mutex_unlock( (dsp.thd)->mutex );
}
}
else { // save FFT at server
if ( (dsp.thd)->fft == OPT_FFT_SERV ) fname_fft = tharg->fname;
else /* OPT_FFT_CLSV */ fname_fft = "db_fft_cl.txt";
fpo = fopen(fname_fft, "wb");
if (fpo != NULL) {
fprintf(fpo, "# <freq/sr>;<dB> ## sr:%d , N:%d\n", dsp.DFT.sr, dsp.DFT.N);
for (j = dsp.DFT.N/2; j < dsp.DFT.N/2 + dsp.DFT.N; j++) {
fprintf(fpo, "%+11.8f;%7.2f\n", bin2fq(&(dsp.DFT), j % dsp.DFT.N), avg_db[j % dsp.DFT.N]);
for (j = 0; j < mlen; j++) {
dsp.DFT.Z[j] *= dsp.DFT.win[j];
}
while (j < dsp.DFT.N) dsp.DFT.Z[j++] = 0.0; // dft(Z[...]) != 0
raw_dft(&(dsp.DFT), dsp.DFT.Z);
for (j = 0; j < dsp.DFT.N; j++) {
float rZ = cabs(dsp.DFT.Z[j]);
avg_rZ[j] += rZ;
all_rZ[j] += rZ;
}
if ( (sum_fft&1)==1 && ((dsp.thd)->fft & OPT_FFT_AVG) == 0 ) {
double t_sec = sum_fft*n_fft*mlen / (double)dsp.sr_base;
for (j = 0; j < dsp.DFT.N; j++) { // if sum_fft odd,
db[j] = 20.0*log10(0.5*all_rZ[j]/dsp.DFT.N+1e-20); // 0.5: rZ_0+rZ_1
all_rZ[j] = 0.0f;
}
// sec.ms, freq_min, freq_max, Hz/bin, N_bins, db_1, ..., db_N
if ( tharg->fpo ) { // save FFT at server
fft_csv_prn(tharg->fpo, &(dsp.DFT), db, t_sec);
}
else if ( tharg->fd > STDIN_FILENO ) {
fft_csv_tcp(tharg->fd, &(dsp.DFT), db, t_sec);
}
fclose(fpo);
}
else {
fprintf(stderr, "error: open %s\n", fname_fft);
}
}
if ( (dsp.thd)->fft != OPT_FFT_SERV ) close(tharg->fd);
(dsp.thd)->fft = 0;
sum_n = 0;
sum_n++;
sum_fft++;
n_fft = 0;
}
if (sum_n*fft_step*mlen >= avg_sec*dsp.sr_base) {
float nN = 1.0/(dsp.DFT.N*(float)sum_n);
for (j = 0; j < dsp.DFT.N; j++) {
avg_db[j] = 20.0*log10(nN*avg_rZ[j]+1e-20);
avg_rZ[j] = 0.0f;
}
if ( (dsp.thd)->fft & OPT_FFT_AVG ) {
double t_sec = sum_fft*fft_step*mlen / (double)dsp.sr_base;
if ( tharg->fpo ) { // send FFT data to client
fft_csv_prn(tharg->fpo, &(dsp.DFT), avg_db, t_sec);
}
else if ( tharg->fd > STDIN_FILENO ) {
fft_csv_tcp(tharg->fd, &(dsp.DFT), avg_db, t_sec);
}
}
n_out++;
if ((dsp.thd)->fft_num > 0 && n_out >= (dsp.thd)->fft_num)
{
if ( tharg->fpo ) { fclose(tharg->fpo); tharg->fpo = NULL; }
else if ( tharg->fd > STDIN_FILENO ) { close(tharg->fd); tharg->fd = -1; }
(dsp.thd)->fft = 0;
sum_fft = 0;
n_out = 0;
pthread_mutex_lock( (dsp.thd)->mutex );
fprintf(FPOUT, "<%d: FFT_STOP>\n", (dsp.thd)->tn);
pthread_mutex_unlock( (dsp.thd)->mutex );
}
sum_n = 0;
}
}
#ifdef FFT_READ_SINK_MIN
@ -445,12 +559,39 @@ static void *thd_FFT(void *targs) {
if ( (dsp.thd)->used == 0 )
{
(dsp.thd)->fft = 0;
pthread_mutex_lock( (dsp.thd)->mutex );
fprintf(FPOUT, "<%d: CLOSE>\n", (dsp.thd)->tn);
pthread_mutex_unlock( (dsp.thd)->mutex );
break;
}
if ( (dsp.thd)->stop_fft > 0 )
{
if ( tharg->fpo ) { fclose(tharg->fpo); tharg->fpo = NULL; }
else if ( tharg->fd > STDIN_FILENO ) { close(tharg->fd); tharg->fd = -1; }
(dsp.thd)->fft = 0;
(dsp.thd)->stop_fft = 0;
sum_fft = 0;
sum_n = 0;
n_out = 0;
pthread_mutex_lock( (dsp.thd)->mutex );
fprintf(FPOUT, "<%d: STOP_FFT>\n", (dsp.thd)->tn);
pthread_mutex_unlock( (dsp.thd)->mutex );
}
}
if ((dsp.thd)->fft_num < 0)
{
if ( tharg->fpo ) { fclose(tharg->fpo); tharg->fpo = NULL; }
else if ( tharg->fd > STDIN_FILENO ) { close(tharg->fd); tharg->fd = -1; }
(dsp.thd)->fft = 0;
pthread_mutex_lock( (dsp.thd)->mutex );
fprintf(FPOUT, "<%d: FFT_STOP>\n", (dsp.thd)->tn);
pthread_mutex_unlock( (dsp.thd)->mutex );
}
if (bitQ == EOF) {
@ -465,6 +606,8 @@ static void *thd_FFT(void *targs) {
exit_thread:
if (z) { free(z); z = NULL; }
if (db) { free(db); db = NULL; }
if (all_rZ) { free(all_rZ); all_rZ = NULL; }
if (avg_rZ) { free(avg_rZ); avg_rZ = NULL; }
if (avg_db) { free(avg_db); avg_db = NULL; }
@ -496,6 +639,8 @@ int main(int argc, char **argv) {
char tcp_buf[TCPBUF_LEN];
int th_used = 0;
int tn_fft = -1;
int opt_fft = 0;
int fft_num = 0;
pcm_t pcm = {0};
@ -510,11 +655,15 @@ int main(int argc, char **argv) {
for (k = 0; k < MAX_FQ; k++) base_fqs[k] = 0.0;
// server options
++argv;
while ((*argv) && (!wavloaded)) {
if (strcmp(*argv, "--dbg") == 0) {
option_dbg = 1;
}
else if (strcmp(*argv, "--enable_clsv_out") == 0) {
option_clsv_out = 1;
}
else if (strcmp(*argv, "--port") == 0) {
int port = 0;
++argv;
@ -524,7 +673,11 @@ int main(int argc, char **argv) {
}
else serv_port = port;
}
else if (strcmp(*argv, "--fft") == 0) {
else if (strncmp(*argv, "--fft", 5) == 0) {
char *arg_fft = *argv;
if (strncmp(arg_fft+5, "_avg", 4) == 0) opt_fft = OPT_FFT_SERV | OPT_FFT_AVG;
else if (strncmp(arg_fft+5, "_all", 4) == 0) opt_fft = OPT_FFT_SERV;
else return -1;
if (xlt_cnt < MAX_FQ) {
base_fqs[xlt_cnt] = 0.0;
rstype[xlt_cnt] = thd_FFT;
@ -532,6 +685,8 @@ int main(int argc, char **argv) {
xlt_cnt++;
}
++argv;
if (*argv) fft_num = atoi(*argv); else return -1;
++argv;
if (*argv) fname_fft = *argv; else return -1;
}
else if (strcmp(*argv, "-") == 0) {
@ -597,8 +752,21 @@ int main(int argc, char **argv) {
for (k = 0; k < xlt_cnt; k++) {
if (k == tn_fft) {
tharg[k].thd.fft = OPT_FFT_SERV;
tharg[k].thd.fft = opt_fft;
tharg[k].fname = fname_fft;
tharg[k].thd.fft_num = fft_num;
tharg[k].thd.stop_fft = 0;
if ( (tharg[k].thd.fft & 0xF) == OPT_FFT_SERV ) { // save FFT at server
if (fname_fft) {
if (fname_fft[0] == '-') tharg[k].fpo = stdout;
else tharg[k].fpo = fopen(fname_fft, "wb");
}
else return -1;
if (tharg[k].fpo == NULL) {
fprintf(stderr, "error: open %s\n", fname_fft);
}
}
}
tharg[k].thd.tn = k;
tharg[k].thd.tn_bit = (1<<k);
@ -686,6 +854,7 @@ int main(int argc, char **argv) {
//if (th_used == 0) break;
// client commands
if ( l > 1 ) {
char *freq = tcp_buf;
while (l > 1 && tcp_buf[l-1] < 0x20) l--;
@ -701,23 +870,76 @@ int main(int argc, char **argv) {
break;
}
else if ( strncmp(tcp_buf, "--fft", 5) == 0 ) {
opt_fft = 0;
if (strncmp(tcp_buf+5, "_avg", 4) == 0) opt_fft = OPT_FFT_AVG;
else if (strncmp(tcp_buf+5, "_all", 4) != 0) return -1;
if (strncmp(tcp_buf+5+4, "_cl", 3) == 0) opt_fft |= OPT_FFT_CLNT;
else if (strncmp(tcp_buf+5+4, "_sv", 3) == 0) {
if (option_clsv_out) opt_fft |= OPT_FFT_SERV;
else {
pthread_mutex_lock( &mutex );
fprintf(FPOUT, "<NO CLSV OUT>\n");
pthread_mutex_unlock( &mutex );
close(conn_fd);
continue;
}
}
else return -1;
if (tcp_buf+5+4+4) fft_num = atoi(tcp_buf+5+4+4); else fft_num = 0;
char *fname_fft_cl = "db_fft_cl.txt";
int opt_fft = strcmp(tcp_buf, "--fft0") == 0 ? OPT_FFT_CLSV : OPT_FFT_CLNT;
//close(conn_fd);
char *pbuf = tcp_buf;
for (pbuf = tcp_buf+5+4+4; *pbuf; pbuf++) {
if (*pbuf == ' ') break;
}
if (*pbuf == ' ') {
fname_fft_cl = pbuf+1;
}
if ( !tcp_eof )
{
if (tn_fft >= 0) {
tharg[tn_fft].thd.fft = opt_fft;
tharg[tn_fft].fname = fname_fft_cl;
tharg[tn_fft].fd = conn_fd;
if (tharg[tn_fft].thd.fft == 0) {
tharg[tn_fft].thd.fft_num = fft_num;
tharg[tn_fft].thd.stop_fft = 0;
tharg[tn_fft].fname = fname_fft_cl;
tharg[tn_fft].fd = conn_fd;
if ( (opt_fft & 0xF) == OPT_FFT_SERV ) { // save FFT at server
if (fname_fft_cl) {
if (fname_fft_cl[0] == '-') tharg[tn_fft].fpo = stdout;
else tharg[tn_fft].fpo = fopen(fname_fft_cl, "wb");
}
else return -1;
if (tharg[tn_fft].fpo == NULL) {
fprintf(stderr, "error: open %s\n", fname_fft_cl);
}
close(tharg[tn_fft].fd); tharg[tn_fft].fd = -1; // tharg[tn_fft].fd == conn_fd
}
else { // (opt_fft & 0xF) == OPT_FFT_CLNT : send FFT to client
tharg[tn_fft].fpo = NULL;
}
tharg[tn_fft].thd.fft = opt_fft;
}
else {
pthread_mutex_lock( &mutex );
fprintf(FPOUT, "<%d: FFT running>\n", tn_fft);
pthread_mutex_unlock( &mutex );
close(conn_fd);
}
}
else {
for (k = 0; k < MAX_FQ; k++) {
if (tharg[k].thd.used == 0) break;
if (tharg[k].thd.used == 0) {
if (k != tn_fft) break;
}
}
if (k < MAX_FQ) {
tharg[k].thd.fft = opt_fft;
tharg[k].thd.fft_num = fft_num;
tharg[k].thd.stop_fft = 0;
tharg[k].fname = fname_fft_cl;
tharg[k].fd = conn_fd;
tn_fft = k;
tharg[k].thd.tn = k;
tharg[k].thd.tn_bit = (1<<k);
@ -726,12 +948,26 @@ int main(int argc, char **argv) {
//tharg[k].thd.lock = &lock;
tharg[k].thd.blk = block_decMB;
tharg[k].thd.xlt_fq = 0.0;
tharg[k].pcm = pcm;
tharg[k].fd = conn_fd;
if ( (opt_fft & 0xF) == OPT_FFT_SERV ) { // save FFT at server
if (fname_fft_cl) {
if (fname_fft_cl[0] == '-') tharg[k].fpo = stdout;
else tharg[k].fpo = fopen(fname_fft_cl, "wb");
}
else return -1;
if (tharg[k].fpo == NULL) {
fprintf(stderr, "error: open %s\n", fname_fft_cl);
}
close(tharg[k].fd); tharg[k].fd = -1; // tharg[k].fd == conn_fd
}
else { // (opt_fft & 0xF) == OPT_FFT_CLNT : send FFT to client
tharg[k].fpo = NULL;
}
rbf1 |= tharg[k].thd.tn_bit;
tharg[k].thd.used = 1;
tharg[k].thd.fft = opt_fft;
pthread_create(&tharg[k].thd.tid, NULL, thd_FFT, &tharg[k]);
@ -753,9 +989,14 @@ int main(int argc, char **argv) {
else if (tcp_buf[0] == '-') { // -<n> : close <n>
int num = atoi(tcp_buf+1);
if (num >= 0 && num < MAX_FQ) {
if (num != tn_fft) {
if (num != tn_fft)
{
tharg[num].thd.used = 0;
}
else
{
tharg[num].thd.stop_fft = 1;
}
}
close(conn_fd);
}
@ -789,6 +1030,7 @@ int main(int argc, char **argv) {
rbf1 |= tharg[k].thd.tn_bit;
tharg[k].thd.used = 1;
tharg[k].thd.fft = 0;
tharg[k].thd.stop_fft = 0;
pthread_create(&tharg[k].thd.tid, NULL, thd_IF, &tharg[k]);

Wyświetl plik

@ -6,6 +6,10 @@
#include <arpa/inet.h>
#define OPT_FFT_SERV 1 // server
#define OPT_FFT_CLNT 2 // server -> client
#define OPT_FFT_AVG 0x100
#define TCPBUF_LEN 1024
#define SERV_BACKLOG 6

Wyświetl plik

@ -8,7 +8,7 @@ import matplotlib.ticker as ticker
if len(sys.argv) < 2:
print("usage:")
print("\tpython %s <fft_file>" % sys.argv[0])
print("\tpython %s <fft_file.csv>" % sys.argv[0])
sys.exit()
fft_file = sys.argv[1]
@ -18,10 +18,17 @@ if not os.path.isfile(fft_file):
sys.exit()
data = np.genfromtxt( fft_file, delimiter=';', names=['fq','db'] , skip_header=1 )
raw_data = np.genfromtxt( fft_file, delimiter=',', max_rows=1)
data1 = raw_data[:] # max_rows=2: raw_data[0,:]
print(data1)
fq = data['fq']
db = data['db']
db = data1[5:]
sr = -2.0*data1[1]
freq_min = data1[1] / sr
freq_max = data1[2] / sr
fq = np.arange(freq_min, freq_max, 1.0/(data1[4]+1))
N = len(db)
m = np.mean(db)

Wyświetl plik

@ -0,0 +1,186 @@
import os
import sys
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
import matplotlib.animation as animation
if len(sys.argv) < 3:
print("usage:")
print("\tpython <1/2/3> %s <fft_all.csv>" % sys.argv[0])
sys.exit()
OPT_FFT_L = 1
OPT_FFT_W = 2
OPT_FFT_B = 3
OPT_FFT = OPT_FFT_B
if (sys.argv[1] == '1'):
OPT_FFT = OPT_FFT_L
elif (sys.argv[1] == '2'):
OPT_FFT = OPT_FFT_W
fft_file = sys.argv[2]
if (fft_file == "-"):
f = sys.stdin
else:
try:
f = open(fft_file)
except IOError:
print("error: open %s" % fft_file)
sys.exit()
FFT_FPS = 16/2
WIN_SEC = 10
win = WIN_SEC*FFT_FPS
line = f.readline()
#line = f.readline()
row = np.fromstring(line, dtype=float, sep=',')
data = row[5:]
l = len(data)
m = np.mean(data)
data5 = np.full([win,l], m)
for x in range(2, l-2):
data5[0,x] = (data[x-2]+data[x-1]+data[x]+data[x+1]+data[x+2])/5.0
min_db = np.min(data5)
max_db = np.max(data5)
sr = -2.0*row[1]
freq_min = row[1]
freq_max = row[2]
N = row[4]
limits = [freq_min/1e3, freq_max/1e3, WIN_SEC, 0.0]
fq = np.arange(freq_min/sr, freq_max/sr, 1.0/(row[4]+1))
################################################################################################
if (OPT_FFT == OPT_FFT_L):
fig = plt.figure(figsize=(12, 5))
ax1 = fig.add_subplot(111)
ax1.set_xlim([fq[0], fq[-1]])
ax1.set_ylim([-110.0, -30.0])
lp, = ax1.plot( fq, data5[0,:], color='g', linewidth=0.4)
count = 0
def animate_lp(i):
line = f.readline()
if line:
row = np.fromstring(line, dtype=float, sep=',')
data = row[5:]
global count
count += 1
global data5
#data5 = np.roll(data5, 1, axis=0)
for x in range(2, l-2):
data5[0,x] = (data[x-2]+data[x-1]+data[x]+data[x+1]+data[x+2])/5.0
lp.set_data(fq, data5[0,:])
return [lp]
ani = animation.FuncAnimation(fig, animate_lp, interval=10, blit=True)
################################################################################################
elif (OPT_FFT == OPT_FFT_W):
fig = plt.figure(figsize=(12, 5))
ax2 = fig.add_subplot(111)
ax2.set_xlabel('Frequency (kHz)')
ax2.set_ylabel('Time (s)')
im = ax2.imshow(data5, vmin=-110.0, vmax=-50.0, extent=limits, animated=True)
ax2.set_aspect('auto')
fig.colorbar(im, orientation='vertical')
count = 0
def animate_im(i):
line = f.readline()
if line:
row = np.fromstring(line, dtype=float, sep=',')
data = row[5:]
global count
count += 1
global data5
data5 = np.roll(data5, 1, axis=0)
for x in range(2, l-2):
data5[0,x] = (data[x-2]+data[x-1]+data[x]+data[x+1]+data[x+2])/5.0
im.set_data(data5)
# update vmin/vmax
if (count % win == 0):
min_db = np.min(data5)
max_db = np.max(data5)
im.set_clim(vmin=min_db, vmax=max_db)
return [im]
ani = animation.FuncAnimation(fig, animate_im, interval=10, blit=True)
################################################################################################
else:
fig = plt.figure(figsize=(12, 8))
ax1 = fig.add_subplot(211)
ax1.set_xlim([fq[0], fq[-1]])
ax1.set_ylim([-110.0, -30.0])
ax2 = fig.add_subplot(212)
ax2.set_xlabel('Frequency (kHz)')
ax2.set_ylabel('Time (s)')
lp, = ax1.plot( fq, data5[0,:], color='g', linewidth=0.4)
im = ax2.imshow(data5, vmin=-110.0, vmax=-50.0, extent=limits, animated=True)
ax2.set_aspect('auto')
count = 0
def animate(i):
line = f.readline()
if line:
row = np.fromstring(line, dtype=float, sep=',')
data = row[5:]
global count
count += 1
global data5
data5 = np.roll(data5, 1, axis=0)
for x in range(2, l-2):
data5[0,x] = (data[x-2]+data[x-1]+data[x]+data[x+1]+data[x+2])/5.0
im.set_data(data5)
lp.set_data(fq, data5[0,:])
# update vmin/vmax
if (count % win == 0):
min_db = np.min(data5)
max_db = np.max(data5)
im.set_clim(vmin=min_db, vmax=max_db)
return [lp,im]
ani = animation.FuncAnimation(fig, animate, interval=10, blit=True)
################################################################################################
plt.show()

Wyświetl plik

@ -0,0 +1,52 @@
import os
import sys
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
if len(sys.argv) < 2:
print("usage:")
print("\tpython %s <wfft_file>" % sys.argv[0])
sys.exit()
fft_file = sys.argv[1]
if not os.path.isfile(fft_file):
print("error: %s not found" % fft_file)
sys.exit()
raw_data = np.genfromtxt( fft_file, delimiter=',')
data = raw_data[:,5:]
dim = np.shape(data)
m = np.mean(data)
data5 = np.full(dim, m)
for y in range(dim[0]):
for x in range(2, dim[1]-2):
data5[y,x] = (data[y,x-2]+data[y,x-1]+data[y,x]+data[y,x+1]+data[y,x+2])/5.0
freq_min = raw_data[0,1] / 1e3
freq_max = raw_data[0,2] / 1e3
N = raw_data[0,4]
tmin = raw_data[dim[0]-1,0]
tmax = raw_data[0,0]
limits = [freq_min, freq_max, tmin, tmax]
fig = plt.figure(figsize=(15, 5))
ax = fig.add_subplot(111)
ax.set_title('Waterfall')
ax.set_xlabel('Frequency (kHz)')
ax.set_ylabel('Time (s)')
plt.imshow(data5, extent=limits)
ax.set_aspect('auto')
plt.colorbar(orientation='vertical')
plt.show()

Wyświetl plik

@ -52,7 +52,7 @@ int main(int argc, char **argv) {
if (argv[1] == NULL) {
fprintf(stderr, "usage:\n");
fprintf(stderr, "\t%s <fft_file>\n", argv[0]);
fprintf(stderr, "\t%s <fft_avg.csv>\n", argv[0]);
return 1;
}
fp = fopen(argv[1], "rb");
@ -61,9 +61,34 @@ int main(int argc, char **argv) {
return 1;
}
memset(line, 0, LINELEN+1);
N = 0;
j = 0;
n = 0;
// sec.ms,freq_min,freq_max,Hz/bin,N_bins, ...
while ( (c = fgetc(fp)) != EOF) {
if (c == '\n') N++;
if (c == '\n') break;
if (c == ' ') continue;
if (c == ',') {
if (n == 1) {
int freq_min = atoi(line);
sr = -2*freq_min;
}
if (n == 4) {
N = atoi(line);
break;
}
n++;
memset(line, 0, LINELEN+1);
j = 0;
}
else {
line[j] = c;
j++;
}
}
db = calloc(N+1, sizeof(float)); if (db == NULL) return 2;
@ -73,20 +98,30 @@ int main(int argc, char **argv) {
intdb = calloc(N+1, sizeof(float)); if (intdb == NULL) return 2;
peak = calloc(N+1, sizeof(float)); if (peak == NULL) return 2;
fseek(fp, 0, SEEK_SET);
pbuf = fgets(line, LINELEN, fp);
p1 = strstr(line, "sr:");
if (p1) sr = atoi(p1+3);
for (n = 0; n < N; n++) {
memset(line, 0, LINELEN+1);
pbuf = fgets(line, LINELEN, fp);
p1 = strstr(line, ";"); //p2 = strstr(p1+1, ";");
if (p1) {
fq[n] = atof(line); //freq[n] = atof(p1+1);
db[n] = atof(p1+1); //atof(p2+1);
// ..., db_1,...,db_N:
memset(line, 0, LINELEN+1);
j = 0;
n = 0;
while ( (c = fgetc(fp)) != EOF) {
if (c == '\n') break;
if (c == ' ') continue;
if (c == ',') {
if (n < N) {
db[n] = atof(line);
fq[n] = -0.5 + n/(float)N;
}
n++;
memset(line, 0, LINELEN+1);
j = 0;
}
else {
line[j] = c;
j++;
}
}
f0 = N/2;
globmin = 0.0;
@ -94,14 +129,12 @@ int main(int argc, char **argv) {
float db_spike3 = 10.0;
int spike_wl3 = 3; //freq2bin(&DFT, 200); // 3 // 200 Hz
int spike_wl5 = 5; //freq2bin(&DFT, 200); // 3 // 200 Hz
//float db_spike1 = 15.0;
//int spike_wl1 = 1; //freq2bin(&DFT, 200); // 3 // 200 Hz
dx = 200.0;
if (sr) dx = sr*(fq[f0+1]-fq[f0]); //freq[f0+1]-freq[f0];
dn = 2*(int)(2400.0/dx)+1; // (odd/symmetric) integration width: 4800+dx Hz
if (option_verbose > 1) fprintf(stderr, "dn = %d\n", dn);
//for (j = 0; j < N; j++) db[j] /= (float)n;
// dc-spike (N-1,)N,0,1(,2): subtract mean/avg
// spikes in general:

Wyświetl plik

@ -15,13 +15,24 @@ simultaneous decoding
&nbsp;&nbsp;&nbsp;&nbsp; `-lm -pthread -o rs_multi`
#### Usage/Examples
`./rs_multi --rs41 <fq0> --dfm <fq1> --m10 <fq2> --lms <fq3> <iq_baseband.wav>` <br />
`./rs_multi --rs41 <fq0> --dfm <fq1> --m10 <fq2> --lms <fq3> - <sr> <bs> <iq_baseband.raw>` <br />
`$ ./rs_multi --rs41 <fq0> --dfm <fq1> --m10 <fq2> --lms <fq3> <iq_baseband.wav>` <br />
`$ ./rs_multi --rs41 <fq0> --dfm <fq1> --m10 <fq2> --lms <fq3> - <sr> <bs> <iq_baseband.raw>` <br />
where <br />
&nbsp;&nbsp;&nbsp;&nbsp; `-0.5 < fq < 0.5`: (relative) frequency, `fq=freq/sr` <br />
&nbsp;&nbsp;&nbsp;&nbsp; `-0.5 < fqX < 0.5`: (relative) frequency, `fq=freq/sr` <br />
&nbsp;&nbsp;&nbsp;&nbsp; `<sr>`: sample rate <br />
&nbsp;&nbsp;&nbsp;&nbsp; `<bs>=8,16,32`: bits per (real) sample (u8, s16 or f32) <br />
decodes up to `MAX_FQ=5 (demod_base.h)` signals. Decoding more signals than number of CPUs/cores is not recommended. <br />
c.f.
https://youtu.be/5YXP9LYUgLs
decodes up to `MAX_FQ=5 (demod_base.h)` signals. Decoding more signals than number of CPUs/cores is not recommended.<br />
Note: If the baseband sample rate has no appropriate factors (e.g. if prime), the IF sample rate might be high and IF-processing slow.<br />
Sending add/remove commands via fifo: <br />
[terminal 1]<br />
`$ rtl_sdr -f 403.0M -s 1920000 - | ./rs_multi --fifo rsfifo --rs41 <fq0> --dfm <fq1> - 1920000 8`<br />
where `<fqX>` is the (relative) frequency of signal `X`<br />
[terminal 2]<br />
*add M10 on `<fq2>`*:<br /> `$ echo "m10 <fq2>" > rsfifo`<br />
*remove `<fq1>`*:<br /> `$ echo "-1" > rsfifo`<br />
If there is no decode for `SEC_NO_SIGNAL=10 (demod_base.c)` seconds, the signal is removed,
unless option `-c` is used.<br />
`--json` output is also possible.

Wyświetl plik

@ -404,9 +404,9 @@ static int f32read_csample(dsp_t *dsp, float complex *z) {
}
static volatile int bufeof = 0; // threads exit
volatile int bufeof = 0; // threads exit
volatile int rbf1;
static volatile int rbf;
extern int rbf1;
static int __f32read_cblock_nocond(dsp_t *dsp) { // blk
@ -417,136 +417,165 @@ static int __f32read_cblock_nocond(dsp_t *dsp) { // blk
if (bufeof) return 0;
pthread_mutex_lock( dsp->thd.mutex );
pthread_mutex_lock( dsp->thd->mutex );
if (rbf == 0)
{
if (dsp->bps == 8) { //uint8
ui8_t u[2*BL];
len = fread( u, dsp->bps/8, 2*BL, dsp->fp) / 2;
for (n = 0; n < len; n++) dsp->thd.blk[n] = (u[2*n]-128)/128.0 + I*(u[2*n+1]-128)/128.0;
for (n = 0; n < len; n++) dsp->thd->blk[n] = (u[2*n]-128)/128.0 + I*(u[2*n+1]-128)/128.0;
}
else if (dsp->bps == 16) { //int16
short b[2*BL];
len = fread( b, dsp->bps/8, 2*BL, dsp->fp) / 2;
for (n = 0; n < len; n++) dsp->thd.blk[n] = b[2*n]/32768.0 + I*b[2*n+1]/32768.0;
for (n = 0; n < len; n++) dsp->thd->blk[n] = b[2*n]/32768.0 + I*b[2*n+1]/32768.0;
}
else { // dsp->bps == 32 //float32
float f[2*BL];
len = fread( f, dsp->bps/8, 2*BL, dsp->fp) / 2;
for (n = 0; n < len; n++) dsp->thd.blk[n] = f[2*n] + I*f[2*n+1];
for (n = 0; n < len; n++) dsp->thd->blk[n] = f[2*n] + I*f[2*n+1];
}
if (len < BL) bufeof = 1;
rbf = rbf1; // set all bits
}
pthread_mutex_unlock( dsp->thd.mutex );
pthread_mutex_unlock( dsp->thd->mutex );
while ((rbf & dsp->thd.tn_bit) == 0) ; // only if #{fqs} leq #{cores} ...
while ((rbf & dsp->thd->tn_bit) == 0) ; // only if #{fqs} leq #{cores} ...
for (n = 0; n < dsp->decM; n++) dsp->decMbuf[n] = dsp->thd.blk[dsp->decM*dsp->blk_cnt + n];
for (n = 0; n < dsp->decM; n++) dsp->decMbuf[n] = dsp->thd->blk[dsp->decM*dsp->blk_cnt + n];
dsp->blk_cnt += 1;
if (dsp->blk_cnt == blk_sz) {
pthread_mutex_lock( dsp->thd.mutex );
rbf &= ~(dsp->thd.tn_bit); // clear bit(tn)
pthread_mutex_lock( dsp->thd->mutex );
rbf &= ~(dsp->thd->tn_bit); // clear bit(tn)
dsp->blk_cnt = 0;
pthread_mutex_unlock( dsp->thd.mutex );
pthread_mutex_unlock( dsp->thd->mutex );
}
return len;
}
static int f32read_cblock(dsp_t *dsp) { // blk_cond
static int f32_cblk(dsp_t *dsp) {
int n;
int BL = dsp->decM * blk_sz;
int len = BL;
float x, y;
if (bufeof) return 0;
if (dsp->bps == 8) { //uint8
ui8_t u[2*BL];
len = fread( u, dsp->bps/8, 2*BL, dsp->fp) / 2;
//for (n = 0; n < len; n++) dsp->thd->blk[n] = (u[2*n]-128)/128.0 + I*(u[2*n+1]-128)/128.0;
// u8: 0..255, 128 -> 0V
for (n = 0; n < len; n++) {
x = (u[2*n ]-128)/128.0;
y = (u[2*n+1]-128)/128.0;
dsp->thd->blk[n] = (x-IQdc.avgIQx) + I*(y-IQdc.avgIQy);
IQdc.sumIQx += x;
IQdc.sumIQy += y;
IQdc.cnt += 1;
if (IQdc.cnt == IQdc.maxcnt) {
IQdc.avgIQx = IQdc.sumIQx/(float)IQdc.maxcnt;
IQdc.avgIQy = IQdc.sumIQy/(float)IQdc.maxcnt;
IQdc.sumIQx = 0; IQdc.sumIQy = 0; IQdc.cnt = 0;
if (IQdc.maxcnt < IQdc.maxlim) IQdc.maxcnt *= 2;
}
}
}
else if (dsp->bps == 16) { //int16
short b[2*BL];
len = fread( b, dsp->bps/8, 2*BL, dsp->fp) / 2;
for (n = 0; n < len; n++) {
x = b[2*n ]/32768.0;
y = b[2*n+1]/32768.0;
dsp->thd->blk[n] = (x-IQdc.avgIQx) + I*(y-IQdc.avgIQy);
IQdc.sumIQx += x;
IQdc.sumIQy += y;
IQdc.cnt += 1;
if (IQdc.cnt == IQdc.maxcnt) {
IQdc.avgIQx = IQdc.sumIQx/(float)IQdc.maxcnt;
IQdc.avgIQy = IQdc.sumIQy/(float)IQdc.maxcnt;
IQdc.sumIQx = 0; IQdc.sumIQy = 0; IQdc.cnt = 0;
if (IQdc.maxcnt < IQdc.maxlim) IQdc.maxcnt *= 2;
}
}
}
else { // dsp->bps == 32 //float32
float f[2*BL];
len = fread( f, dsp->bps/8, 2*BL, dsp->fp) / 2;
for (n = 0; n < len; n++) {
x = f[2*n];
y = f[2*n+1];
dsp->thd->blk[n] = (x-IQdc.avgIQx) + I*(y-IQdc.avgIQy);
IQdc.sumIQx += x;
IQdc.sumIQy += y;
IQdc.cnt += 1;
if (IQdc.cnt == IQdc.maxcnt) {
IQdc.avgIQx = IQdc.sumIQx/(float)IQdc.maxcnt;
IQdc.avgIQy = IQdc.sumIQy/(float)IQdc.maxcnt;
IQdc.sumIQx = 0; IQdc.sumIQy = 0; IQdc.cnt = 0;
if (IQdc.maxcnt < IQdc.maxlim) IQdc.maxcnt *= 2;
}
}
}
if (len < BL) bufeof = 1;
pthread_mutex_lock( dsp->thd.mutex );
return len;
}
static int f32read_cblock(dsp_t *dsp) { // blk_cond
int n;
int len = dsp->decM;
if (bufeof) return 0;
//if (dsp->thd->used == 0) { }
pthread_mutex_lock( dsp->thd->mutex );
if (rbf == 0)
{
if (dsp->bps == 8) { //uint8
ui8_t u[2*BL];
len = fread( u, dsp->bps/8, 2*BL, dsp->fp) / 2;
//for (n = 0; n < len; n++) dsp->thd.blk[n] = (u[2*n]-128)/128.0 + I*(u[2*n+1]-128)/128.0;
// u8: 0..255, 128 -> 0V
for (n = 0; n < len; n++) {
x = (u[2*n ]-128)/128.0;
y = (u[2*n+1]-128)/128.0;
dsp->thd.blk[n] = (x-IQdc.avgIQx) + I*(y-IQdc.avgIQy);
IQdc.sumIQx += x;
IQdc.sumIQy += y;
IQdc.cnt += 1;
if (IQdc.cnt == IQdc.maxcnt) {
IQdc.avgIQx = IQdc.sumIQx/(float)IQdc.maxcnt;
IQdc.avgIQy = IQdc.sumIQy/(float)IQdc.maxcnt;
IQdc.sumIQx = 0; IQdc.sumIQy = 0; IQdc.cnt = 0;
if (IQdc.maxcnt < IQdc.maxlim) IQdc.maxcnt *= 2;
}
}
}
else if (dsp->bps == 16) { //int16
short b[2*BL];
len = fread( b, dsp->bps/8, 2*BL, dsp->fp) / 2;
for (n = 0; n < len; n++) {
x = b[2*n ]/32768.0;
y = b[2*n+1]/32768.0;
dsp->thd.blk[n] = (x-IQdc.avgIQx) + I*(y-IQdc.avgIQy);
IQdc.sumIQx += x;
IQdc.sumIQy += y;
IQdc.cnt += 1;
if (IQdc.cnt == IQdc.maxcnt) {
IQdc.avgIQx = IQdc.sumIQx/(float)IQdc.maxcnt;
IQdc.avgIQy = IQdc.sumIQy/(float)IQdc.maxcnt;
IQdc.sumIQx = 0; IQdc.sumIQy = 0; IQdc.cnt = 0;
if (IQdc.maxcnt < IQdc.maxlim) IQdc.maxcnt *= 2;
}
}
}
else { // dsp->bps == 32 //float32
float f[2*BL];
len = fread( f, dsp->bps/8, 2*BL, dsp->fp) / 2;
for (n = 0; n < len; n++) {
x = f[2*n];
y = f[2*n+1];
dsp->thd.blk[n] = (x-IQdc.avgIQx) + I*(y-IQdc.avgIQy);
IQdc.sumIQx += x;
IQdc.sumIQy += y;
IQdc.cnt += 1;
if (IQdc.cnt == IQdc.maxcnt) {
IQdc.avgIQx = IQdc.sumIQx/(float)IQdc.maxcnt;
IQdc.avgIQy = IQdc.sumIQy/(float)IQdc.maxcnt;
IQdc.sumIQx = 0; IQdc.sumIQy = 0; IQdc.cnt = 0;
if (IQdc.maxcnt < IQdc.maxlim) IQdc.maxcnt *= 2;
}
}
}
if (len < BL) bufeof = 1;
len = f32_cblk(dsp);
rbf = rbf1; // set all bits
pthread_cond_broadcast( dsp->thd.cond );
pthread_cond_broadcast( dsp->thd->cond );
}
while ((rbf & dsp->thd.tn_bit) == 0) pthread_cond_wait( dsp->thd.cond, dsp->thd.mutex );
while ((rbf & dsp->thd->tn_bit) == 0) pthread_cond_wait( dsp->thd->cond, dsp->thd->mutex );
for (n = 0; n < dsp->decM; n++) dsp->decMbuf[n] = dsp->thd.blk[dsp->decM*dsp->blk_cnt + n];
for (n = 0; n < dsp->decM; n++) dsp->decMbuf[n] = dsp->thd->blk[dsp->decM*dsp->blk_cnt + n];
dsp->blk_cnt += 1;
if (dsp->blk_cnt == blk_sz) {
rbf &= ~(dsp->thd.tn_bit); // clear bit(tn)
rbf &= ~(dsp->thd->tn_bit); // clear bit(tn)
dsp->blk_cnt = 0;
}
pthread_mutex_unlock( dsp->thd.mutex );
pthread_mutex_unlock( dsp->thd->mutex );
return len;
}
int reset_blockread(dsp_t *dsp) {
int len = 0;
pthread_mutex_lock( dsp->thd->mutex );
rbf1 &= ~(dsp->thd->tn_bit);
if ( (rbf & dsp->thd->tn_bit) == dsp->thd->tn_bit )
{
len = f32_cblk(dsp);
rbf = rbf1; // set all bits
pthread_cond_broadcast( dsp->thd->cond );
}
pthread_mutex_unlock( dsp->thd->mutex );
return len;
}
// decimate lowpass
static float *ws_dec;
@ -654,6 +683,7 @@ int f32buf_sample(dsp_t *dsp, int inv) {
ui32_t s_reset = dsp->dectaps*dsp->lut_len;
int j;
if ( f32read_cblock(dsp) < dsp->decM ) return EOF;
//if ( f32read_cblock(dsp) < dsp->decM * blk_sz) return EOF;
for (j = 0; j < dsp->decM; j++) {
dsp->decXbuffer[dsp->sample_dec % dsp->dectaps] = dsp->decMbuf[j] * dsp->ex[dsp->sample_dec % dsp->lut_len];
dsp->sample_dec += 1;
@ -749,14 +779,14 @@ int f32buf_sample(dsp_t *dsp, int inv) {
if (inv) s = -s;
dsp->bufs[dsp->sample_in % dsp->M] = s;
/*
xneu = dsp->bufs[(dsp->sample_in ) % dsp->M];
xalt = dsp->bufs[(dsp->sample_in+dsp->M - dsp->Nvar) % dsp->M];
dsp->xsum += xneu - xalt; // + xneu - xalt
dsp->qsum += (xneu - xalt)*(xneu + xalt); // + xneu*xneu - xalt*xalt
dsp->xs[dsp->sample_in % dsp->M] = dsp->xsum;
dsp->qs[dsp->sample_in % dsp->M] = dsp->qsum;
*/
dsp->sample_out = dsp->sample_in - dsp->delay;
@ -870,7 +900,7 @@ int read_slbit(dsp_t *dsp, int *bit, int inv, int ofs, int pos, float l, int spi
}
sample -= dc;
if ( l < 0 || (mid-l < dsp->sc && dsp->sc < mid+l)) sum -= sample;
if (l < 0 || (mid-l < dsp->sc && dsp->sc < mid+l)) sum -= sample;
dsp->sc++;
} while (dsp->sc < bg); // n < dsp->sps
@ -888,9 +918,9 @@ int read_slbit(dsp_t *dsp, int *bit, int inv, int ofs, int pos, float l, int spi
+dsp->bufs[(dsp->sample_out-dsp->buffered+1 + ofs + dsp->M) % dsp->M]);
sample = avg + scale*(sample - avg); // spikes
}
sample -= dc;
sample -= dc;
if ( l < 0 || (mid-l < dsp->sc && dsp->sc < mid+l)) sum += sample;
if (l < 0 || (mid-l < dsp->sc && dsp->sc < mid+l)) sum += sample;
dsp->sc++;
} while (dsp->sc < bg); // n < dsp->sps
@ -902,6 +932,81 @@ int read_slbit(dsp_t *dsp, int *bit, int inv, int ofs, int pos, float l, int spi
return 0;
}
int read_softbit(dsp_t *dsp, hsbit_t *shb, int inv, int ofs, int pos, float l, int spike) {
// symlen==2: manchester2 10->0,01->1: 2.bit
float sample;
float avg;
float ths = 0.5, scale = 0.27;
double sum = 0.0;
double mid;
//double l = 1.0;
double bg = pos*dsp->symlen*dsp->sps;
double dc = 0.0;
ui8_t bit = 0;
if (dsp->opt_dc && dsp->opt_iq < 2) dc = dsp->dc;
if (pos == 0) {
bg = 0;
dsp->sc = 0;
}
if (dsp->symlen == 2) {
mid = bg + (dsp->sps-1)/2.0;
bg += dsp->sps;
do {
if (dsp->buffered > 0) dsp->buffered -= 1;
else if (f32buf_sample(dsp, inv) == EOF) return EOF;
sample = dsp->bufs[(dsp->sample_out-dsp->buffered + ofs + dsp->M) % dsp->M];
if (spike && fabs(sample - avg) > ths) {
avg = 0.5*(dsp->bufs[(dsp->sample_out-dsp->buffered-1 + ofs + dsp->M) % dsp->M]
+dsp->bufs[(dsp->sample_out-dsp->buffered+1 + ofs + dsp->M) % dsp->M]);
sample = avg + scale*(sample - avg); // spikes
}
sample -= dc;
if (l < 0 || (mid-l < dsp->sc && dsp->sc < mid+l)) sum -= sample;
dsp->sc++;
} while (dsp->sc < bg); // n < dsp->sps
}
mid = bg + (dsp->sps-1)/2.0;
bg += dsp->sps;
do {
if (dsp->buffered > 0) dsp->buffered -= 1;
else if (f32buf_sample(dsp, inv) == EOF) return EOF;
sample = dsp->bufs[(dsp->sample_out-dsp->buffered + ofs + dsp->M) % dsp->M];
if (spike && fabs(sample - avg) > ths) {
avg = 0.5*(dsp->bufs[(dsp->sample_out-dsp->buffered-1 + ofs + dsp->M) % dsp->M]
+dsp->bufs[(dsp->sample_out-dsp->buffered+1 + ofs + dsp->M) % dsp->M]);
sample = avg + scale*(sample - avg); // spikes
}
sample -= dc;
if (l < 0 || (mid-l < dsp->sc && dsp->sc < mid+l)) sum += sample;
dsp->sc++;
} while (dsp->sc < bg); // n < dsp->sps
if (sum >= 0) bit = 1;
else bit = 0;
shb->hb = bit;
shb->sb = (float)sum;
return 0;
}
/* -------------------------------------------------------------------------- */
#define IF_TRANSITION_BW (4e3) // 4kHz transition width
@ -952,7 +1057,7 @@ int init_buffers(dsp_t *dsp) {
// lookup table, exp-rotation
int W = 2*8; // 16 Hz window
int d = 1; // 1..W , groesster Teiler d <= W von sr_base
int freq = (int)( dsp->thd.xlt_fq * (double)dsp->sr_base + 0.5);
int freq = (int)( dsp->thd->xlt_fq * (double)dsp->sr_base + 0.5);
int freq0 = freq; // init
double f0 = freq0 / (double)dsp->sr_base; // init
@ -1037,6 +1142,7 @@ int init_buffers(dsp_t *dsp) {
dsp->delay = L/16;
dsp->sample_in = 0;
dsp->last_detect = 0;
p2 = 1;
while (p2 < M) p2 <<= 1;
@ -1059,10 +1165,10 @@ int init_buffers(dsp_t *dsp) {
dsp->bufs = (float *)calloc( M+1, sizeof(float)); if (dsp->bufs == NULL) return -100;
dsp->match = (float *)calloc( L+1, sizeof(float)); if (dsp->match == NULL) return -100;
/*
dsp->xs = (float *)calloc( M+1, sizeof(float)); if (dsp->xs == NULL) return -100;
dsp->qs = (float *)calloc( M+1, sizeof(float)); if (dsp->qs == NULL) return -100;
*/
dsp->rawbits = (char *)calloc( 2*dsp->hdrlen+1, sizeof(char)); if (dsp->rawbits == NULL) return -100;
@ -1143,10 +1249,11 @@ int free_buffers(dsp_t *dsp) {
if (dsp->match) { free(dsp->match); dsp->match = NULL; }
if (dsp->bufs) { free(dsp->bufs); dsp->bufs = NULL; }
if (dsp->rawbits) { free(dsp->rawbits); dsp->rawbits = NULL; }
/*
if (dsp->xs) { free(dsp->xs); dsp->xs = NULL; }
if (dsp->qs) { free(dsp->qs); dsp->qs = NULL; }
if (dsp->rawbits) { free(dsp->rawbits); dsp->rawbits = NULL; }
*/
if (dsp->DFT.xn) { free(dsp->DFT.xn); dsp->DFT.xn = NULL; }
if (dsp->DFT.ew) { free(dsp->DFT.ew); dsp->DFT.ew = NULL; }
if (dsp->DFT.Fm) { free(dsp->DFT.Fm); dsp->DFT.Fm = NULL; }
@ -1196,6 +1303,7 @@ ui32_t get_sample(dsp_t *dsp) {
/* ------------------------------------------------------------------------------------ */
#define SEC_NO_SIGNAL 10
int find_header(dsp_t *dsp, float thres, int hdmax, int bitofs, int opt_dc) {
ui32_t k = 0;
@ -1212,6 +1320,16 @@ int find_header(dsp_t *dsp, float thres, int hdmax, int bitofs, int opt_dc) {
mp = getCorrDFT(dsp); // correlation score -> dsp->mv
//if (option_auto == 0 && dsp->mv < 0) mv = 0;
k = 0;
// signal lost
if ( dsp->thd->used == 0 ||
(!dsp->opt_cnt && dsp->mv_pos - dsp->last_detect > SEC_NO_SIGNAL*dsp->sr) )
{
pthread_mutex_lock( dsp->thd->mutex );
fprintf(stdout, "<%d: close>\n", dsp->thd->tn);
pthread_mutex_unlock( dsp->thd->mutex );
return EOF;
}
}
else {
dsp->mv = 0.0;
@ -1249,6 +1367,8 @@ int find_header(dsp_t *dsp, float thres, int hdmax, int bitofs, int opt_dc) {
herrs = headcmp(dsp, opt_dc);
if (herrs <= hdmax) header_found = 1; // max bitfehler in header
dsp->last_detect = dsp->mv_pos;
if (header_found) return 1;
}
}

Wyświetl plik

@ -36,6 +36,7 @@ typedef struct {
int max_fq;
double xlt_fq;
float complex *blk;
int used;
} thd_t;
@ -80,13 +81,16 @@ typedef struct {
float *bufs;
float mv;
ui32_t mv_pos;
ui32_t last_detect;
//
int N_norm;
int Nvar;
/*
float xsum;
float qsum;
float *xs;
float *qs;
*/
// IQ-data
int opt_iq;
@ -149,7 +153,9 @@ typedef struct {
float *lpFM_buf;
float *fm_buffer;
thd_t thd;
int opt_cnt;
thd_t *thd;
} dsp_t;
@ -167,12 +173,18 @@ typedef struct {
} pcm_t;
typedef struct {
ui8_t hb;
float sb;
} hsbit_t;
typedef struct {
pcm_t pcm;
thd_t thd;
int option_jsn;
int option_dc;
int option_cnt;
int jsn_freq;
} thargs_t;
@ -181,6 +193,7 @@ typedef struct {
float read_wav_header(pcm_t *);
int f32buf_sample(dsp_t *, int);
int read_slbit(dsp_t *, int*, int, int, int, float, int);
int read_softbit(dsp_t *, hsbit_t *, int, int, int, float, int );
int init_buffers(dsp_t *);
int free_buffers(dsp_t *);
@ -193,3 +206,6 @@ int decimate_init(float f, int taps);
int decimate_free(void);
int iq_dc_init(pcm_t *);
int reset_blockread(dsp_t *);

Wyświetl plik

@ -772,21 +772,26 @@ static int print_frame(gpx_t *gpx, dsp_t *dsp) {
if (ret1 == 0 || ret1 > 0) {
frid = dat_out(gpx, block_dat1, ret1);
if (frid == 8) {
pthread_mutex_lock( dsp->thd.mutex );
fprintf(stdout, "<%d> ", dsp->thd.tn);
pthread_mutex_lock( dsp->thd->mutex );
fprintf(stdout, "<%d> ", dsp->thd->tn);
ret1 = print_gpx(gpx);
if (ret1==0) fprintf(stdout, "\n");
pthread_mutex_unlock( dsp->thd.mutex );
pthread_mutex_unlock( dsp->thd->mutex );
}
}
if (ret2 == 0 || ret2 > 0) {
frid = dat_out(gpx, block_dat2, ret2);
if (frid == 8) {
pthread_mutex_lock( dsp->thd.mutex );
fprintf(stdout, "<%d> ", dsp->thd.tn);
pthread_mutex_lock( dsp->thd->mutex );
//fprintf(stdout, "<%d> ", dsp->thd->tn);
fprintf(stdout, "<%d: ", dsp->thd->tn);
fprintf(stdout, "s=%+.4f, ", dsp->mv);
fprintf(stdout, "f=%+.4f", -dsp->thd->xlt_fq);
if (dsp->opt_dc) fprintf(stdout, "%+.6f", dsp->Df/(double)dsp->sr);
fprintf(stdout, "> ");
ret2 = print_gpx(gpx);
if (ret2==0) fprintf(stdout, "\n");
pthread_mutex_unlock( dsp->thd.mutex );
pthread_mutex_unlock( dsp->thd->mutex );
}
}
@ -879,7 +884,7 @@ void *thd_dfm09(void *targs) {
dsp.dectaps = pcm->dectaps;
dsp.decM = pcm->decM;
dsp.thd = tharg->thd;
dsp.thd = &(tharg->thd);
dsp.bps = pcm->bps;
dsp.nch = pcm->nch;
@ -897,17 +902,18 @@ void *thd_dfm09(void *targs) {
dsp.opt_lp = 1;
dsp.lpIQ_bw = 12e3; // IF lowpass bandwidth
dsp.lpFM_bw = 4e3; // FM audio lowpass
dsp.opt_dc = tharg->option_dc;
dsp.opt_dc = tharg->option_dc;
dsp.opt_cnt = tharg->option_cnt;
if ( dsp.sps < 8 ) {
fprintf(stderr, "note: sample rate low\n");
//fprintf(stderr, "note: sample rate low\n");
}
k = init_buffers(&dsp);
if ( k < 0 ) {
fprintf(stderr, "error: init buffers\n");
return NULL;
goto exit_thread;
};
@ -975,6 +981,9 @@ void *thd_dfm09(void *targs) {
free_buffers(&dsp);
exit_thread:
reset_blockread(&dsp);
(dsp.thd)->used = 0;
return NULL;
}

Wyświetl plik

@ -107,12 +107,11 @@ typedef struct {
ui8_t bIn;
ui8_t codeIn;
ui8_t prevState; // 0..M=64
int w; // > 255 : if (w>250): w=250 ?
//float sw;
float w;
} states_t;
typedef struct {
char rawbits[RAWBITFRAME_LEN+OVERLAP*BITS*2 +8];
hsbit_t rawbits[RAWBITFRAME_LEN+OVERLAP*BITS*2 +8];
states_t state[RAWBITFRAME_LEN+OVERLAP +8][M];
states_t d[N];
} VIT_t;
@ -129,7 +128,7 @@ typedef struct {
double lat; double lon; double alt;
double vH; double vD; double vV;
double vE; double vN; double vU;
char blk_rawbits[RAWBITBLOCK_LEN+SYNC_LEN*BITS*2 +9];
hsbit_t blk_rawbits[RAWBITBLOCK_LEN+SYNC_LEN*BITS*2 +9];
ui8_t frame[FRM_LEN]; // = { 0x24, 0x54, 0x00, 0x00}; // dataheader
int frm_pos; // ecc_blk <-> frm_blk
int sf6;
@ -205,12 +204,16 @@ static int vit_initCodes(gpx_t *gpx) {
return 0;
}
static int vit_dist(int c, char *rc) {
return (((c>>1)^rc[0])&1) + ((c^rc[1])&1);
static float vit_dist2(int c, hsbit_t *rc) {
int c0 = 2*((c>>1) & 1)-1; // {0,1} -> {-1,+1}
int c1 = 2*(c & 1)-1;
float d2 = (c0-rc[0].sb)*(c0-rc[0].sb) + (c1-rc[1].sb)*(c1-rc[1].sb);
return d2;
}
static int vit_start(VIT_t *vit, char *rc) {
int t, m, j, c, d;
static int vit_start(VIT_t *vit, hsbit_t *rc) {
int t, m, j, c;
float d;
t = L-1;
m = M;
@ -228,7 +231,7 @@ static int vit_start(VIT_t *vit, char *rc) {
c = vit_code[j];
vit->state[t][j].bIn = j % 2;
vit->state[t][j].codeIn = c;
d = vit_dist( c, rc+2*(t-1) );
d = vit_dist2( c, rc+2*(t-1) );
vit->state[t][j].w = vit->state[t-1][vit->state[t][j].prevState].w + d;
}
m *= 2;
@ -237,7 +240,7 @@ static int vit_start(VIT_t *vit, char *rc) {
return t;
}
static int vit_next(VIT_t *vit, int t, char *rc) {
static int vit_next(VIT_t *vit, int t, hsbit_t *rc) {
int b, nstate;
int j, index;
@ -247,7 +250,7 @@ static int vit_next(VIT_t *vit, int t, char *rc) {
vit->d[nstate].bIn = b;
vit->d[nstate].codeIn = vit_code[nstate];
vit->d[nstate].prevState = j;
vit->d[nstate].w = vit->state[t][j].w + vit_dist( vit->d[nstate].codeIn, rc );
vit->d[nstate].w = vit->state[t][j].w + vit_dist2( vit->d[nstate].codeIn, rc );
}
}
@ -264,11 +267,11 @@ static int vit_next(VIT_t *vit, int t, char *rc) {
static int vit_path(VIT_t *vit, int j, int t) {
int c;
vit->rawbits[2*t] = '\0';
vit->rawbits[2*t].hb = '\0';
while (t > 0) {
c = vit->state[t][j].codeIn;
vit->rawbits[2*t -2] = 0x30 + ((c>>1) & 1);
vit->rawbits[2*t -1] = 0x30 + (c & 1);
vit->rawbits[2*t -2].hb = 0x30 + ((c>>1) & 1);
vit->rawbits[2*t -1].hb = 0x30 + (c & 1);
j = vit->state[t][j].prevState;
t--;
}
@ -276,13 +279,20 @@ static int vit_path(VIT_t *vit, int j, int t) {
return 0;
}
static int viterbi(VIT_t *vit, char *rc) {
static int hbstr_len(hsbit_t *hsbit) {
int len = 0;
while (hsbit[len].hb) len++;
return len;
}
static int viterbi(VIT_t *vit, hsbit_t *rc) {
int t, tmax;
int j, j_min, w_min;
int j, j_min;
float w_min;
vit_start(vit, rc);
tmax = strlen(rc)/2;
tmax = hbstr_len(rc)/2;
for (t = L-1; t < tmax; t++)
{
@ -307,15 +317,15 @@ static int viterbi(VIT_t *vit, char *rc) {
// ------------------------------------------------------------------------
static int deconv(char* rawbits, char *bits) {
static int deconv(hsbit_t *rawbits, char *bits) {
int j, n, bitA, bitB;
char *p;
hsbit_t *p;
int len;
int errors = 0;
int m = L-1;
len = strlen(rawbits);
len = hbstr_len(rawbits);
for (j = 0; j < m; j++) bits[j] = '0';
n = 0;
while ( 2*(m+n) < len ) {
@ -325,10 +335,10 @@ static int deconv(char* rawbits, char *bits) {
bitA ^= (bits[n+j]&1) & (polyA[j]&1);
bitB ^= (bits[n+j]&1) & (polyB[j]&1);
}
if ( (bitA^(p[0]&1))==(polyA[m]&1) && (bitB^(p[1]&1))==(polyB[m]&1) ) bits[n+m] = '1';
else if ( (bitA^(p[0]&1))==0 && (bitB^(p[1]&1))==0 ) bits[n+m] = '0';
if ( (bitA^(p[0].hb&1))==(polyA[m]&1) && (bitB^(p[1].hb&1))==(polyB[m]&1) ) bits[n+m] = '1';
else if ( (bitA^(p[0].hb&1))==0 && (bitB^(p[1].hb&1))==0 ) bits[n+m] = '0';
else {
if ( (bitA^(p[0]&1))!=(polyA[m]&1) && (bitB^(p[1]&1))==(polyB[m]&1) ) bits[n+m] = 0x39;
if ( (bitA^(p[0].hb&1))!=(polyA[m]&1) && (bitB^(p[1].hb&1))==(polyB[m]&1) ) bits[n+m] = 0x39;
else bits[n+m] = 0x38;
errors = n;
break;
@ -754,11 +764,16 @@ static int print_frame(gpx_t *gpx, int crc_err, int len) {
static int print_thd_frame(gpx_t *gpx, int crc_err, int len, dsp_t *dsp) {
int ret = 0;
pthread_mutex_lock( dsp->thd.mutex );
printf("<%d> ", dsp->thd.tn);
pthread_mutex_lock( dsp->thd->mutex );
//printf("<%d> ", dsp->thd->tn);
fprintf(stdout, "<%d: ", dsp->thd->tn);
fprintf(stdout, "s=%+.4f, ", dsp->mv);
fprintf(stdout, "f=%+.4f", -dsp->thd->xlt_fq);
if (dsp->opt_dc) fprintf(stdout, "%+.6f", dsp->Df/(double)dsp->sr);
fprintf(stdout, "> ");
ret = print_frame(gpx, crc_err, len);
if (ret==0) printf("\n");
pthread_mutex_unlock( dsp->thd.mutex );
if (ret==0) fprintf(stdout, "\n");
pthread_mutex_unlock( dsp->thd->mutex );
return ret;
}
@ -803,7 +818,7 @@ static void proc_frame(gpx_t *gpx, int len, dsp_t *dsp) {
ui8_t block_bytes[FRAME_LEN+8];
ui8_t rs_cw[rs_N];
char frame_bits[BITFRAME_LEN+OVERLAP*BITS +8]; // init L-1 bits mit 0
char *rawbits = NULL;
hsbit_t *rawbits = NULL;
int i, j;
int err = 0;
int errs = 0;
@ -812,13 +827,17 @@ static void proc_frame(gpx_t *gpx, int len, dsp_t *dsp) {
if ((len % 8) > 4) {
while (len % 8) gpx->blk_rawbits[len++] = '0';
while (len % 8) {
gpx->blk_rawbits[len].hb = '0';
gpx->blk_rawbits[len].sb = -1;
len++;
}
}
gpx->blk_rawbits[len] = '\0';
gpx->blk_rawbits[len].hb = '\0';
flen = len / (2*BITS);
if (gpx->option.vit == 1) {
if (gpx->option.vit) {
viterbi(gpx->vit, gpx->blk_rawbits);
rawbits = gpx->vit->rawbits;
}
@ -948,14 +967,14 @@ void *thd_lms6X(void *targs) { // pcm_t *pcm, double xlt_fq
int k;
int bit, rbit;
hsbit_t hsbit, rhsbit;
int bitpos = 0;
int bitQ = 0;
int pos;
int header_found = 0;
float thres = 0.68;
float thres = 0.65;
float _mv = 0.0;
int symlen = 1;
@ -996,17 +1015,21 @@ void *thd_lms6X(void *targs) { // pcm_t *pcm, double xlt_fq
//gpx->option.ptu = 1;
gpx->option.jsn = tharg->option_jsn;
gpx->option.vit = 1;
gpx->option.vit = 2;
gpx->option.ecc = 1;
gpx->jsn_freq = tharg->jsn_freq;
memcpy(gpx->blk_rawbits, blk_syncbits, sizeof(blk_syncbits));
memcpy(gpx->frame, frm_sync6, sizeof(frm_sync6));
gpx->frm_pos = 0; // ecc_blk <-> frm_blk
gpx->sf6 = 0;
gpx->sfX = 0;
//memcpy(gpx->blk_rawbits, blk_syncbits, sizeof(blk_syncbits));
for (k = 0; k < strlen(blk_syncbits); k++) { // strlen(blk_syncbits)=BLOCKSTART
int hbit = blk_syncbits[k] & 1;
gpx->blk_rawbits[k].hb = hbit + 0x30;
gpx->blk_rawbits[k].sb = 2*hbit-1;
}
gpx->week = gpsweek;
@ -1023,7 +1046,7 @@ void *thd_lms6X(void *targs) { // pcm_t *pcm, double xlt_fq
dsp.dectaps = pcm->dectaps;
dsp.decM = pcm->decM;
dsp.thd = tharg->thd;
dsp.thd = &(tharg->thd);
dsp.bps = pcm->bps;
dsp.nch = pcm->nch;
@ -1041,23 +1064,26 @@ void *thd_lms6X(void *targs) { // pcm_t *pcm, double xlt_fq
dsp.opt_lp = 1;
dsp.lpIQ_bw = 8e3; // IF lowpass bandwidth
dsp.lpFM_bw = 6e3; // FM audio lowpass
dsp.opt_dc = tharg->option_dc;
dsp.opt_dc = tharg->option_dc;
dsp.opt_cnt = tharg->option_cnt;
if ( dsp.sps < 8 ) {
fprintf(stderr, "note: sample rate low (%.1f sps)\n", dsp.sps);
//fprintf(stderr, "note: sample rate low (%.1f sps)\n", dsp.sps);
}
k = init_buffers(&dsp); // baud difference not significant
if ( k < 0 ) {
fprintf(stderr, "error: init buffers\n");
return NULL;
goto exit_thread;
};
if (gpx->option.vit) {
k = vit_initCodes(gpx);
if (k < 0) return NULL;
if (k < 0) {
goto exit_thread;
}
}
if (gpx->option.ecc) {
rs_init_RS255ccsds(&gpx->RS); // bch_ecc.c
@ -1090,7 +1116,7 @@ void *thd_lms6X(void *targs) { // pcm_t *pcm, double xlt_fq
bitQ = 0;
while ( 1 && bitQ != EOF )
{
header_found = find_header(&dsp, thres, 6, bitofs, dsp.opt_dc);
header_found = find_header(&dsp, thres, 10, bitofs, dsp.opt_dc);
_mv = dsp.mv;
if (header_found == EOF) break;
@ -1115,19 +1141,26 @@ void *thd_lms6X(void *targs) { // pcm_t *pcm, double xlt_fq
while ( pos < rawbitblock_len ) {
bitQ = read_slbit(&dsp, &rbit, 0, bitofs, bitpos, -1, 0); // symlen=1
bitQ = read_softbit(&dsp, &rhsbit, 0, bitofs, bitpos, -1, 0); // symlen=1
if (bitQ == EOF) { break; }
bit = rbit ^ (bc%2); // (c0,inv(c1))
gpx->blk_rawbits[pos] = 0x30 + bit;
hsbit.hb = rhsbit.hb ^ (bc%2); // (c0,inv(c1))
int sgn = -2*(((unsigned int)bc)%2)+1;
hsbit.sb = sgn * rhsbit.sb;
if (gpx->option.vit == 1) { // hard decision
hsbit.sb = 2*hsbit.hb -1;
}
gpx->blk_rawbits[pos] = hsbit;
gpx->blk_rawbits[pos].hb += 0x30;
bc++;
pos++;
bitpos += 1;
}
gpx->blk_rawbits[pos] = '\0';
gpx->blk_rawbits[pos].hb = '\0';
time_elapsed_sec = dsp.sample_in / (double)dsp.sr;
proc_frame(gpx, pos, &dsp);
@ -1173,6 +1206,9 @@ void *thd_lms6X(void *targs) { // pcm_t *pcm, double xlt_fq
free_buffers(&dsp);
if (gpx->vit) { free(gpx->vit); gpx->vit = NULL; }
exit_thread:
reset_blockread(&dsp);
(dsp.thd)->used = 0;
return NULL;
}

Wyświetl plik

@ -1101,11 +1101,16 @@ static int print_frame(gpx_t *gpx, int pos, dsp_t *dsp) {
}
else {
int ret = 0;
pthread_mutex_lock( dsp->thd.mutex );
fprintf(stdout, "<%d> ", dsp->thd.tn);
pthread_mutex_lock( dsp->thd->mutex );
//fprintf(stdout, "<%d> ", dsp->thd->tn);
fprintf(stdout, "<%d: ", dsp->thd->tn);
fprintf(stdout, "s=%+.4f, ", dsp->mv);
fprintf(stdout, "f=%+.4f", -dsp->thd->xlt_fq);
if (dsp->opt_dc) fprintf(stdout, "%+.6f", dsp->Df/(double)dsp->sr);
fprintf(stdout, "> ");
ret = print_pos(gpx, cs1 == cs2);
if (ret==0) fprintf(stdout, "\n");
pthread_mutex_unlock( dsp->thd.mutex );
pthread_mutex_unlock( dsp->thd->mutex );
}
return (gpx->frame_bytes[0]<<8)|gpx->frame_bytes[1];
@ -1179,7 +1184,7 @@ void *thd_m10(void *targs) { // pcm_t *pcm, double xlt_fq
dsp.dectaps = pcm->dectaps;
dsp.decM = pcm->decM;
dsp.thd = tharg->thd;
dsp.thd = &(tharg->thd);
dsp.bps = pcm->bps;
dsp.nch = pcm->nch;
@ -1197,10 +1202,11 @@ void *thd_m10(void *targs) { // pcm_t *pcm, double xlt_fq
dsp.opt_lp = 1;
dsp.lpIQ_bw = 24e3; // IF lowpass bandwidth
dsp.lpFM_bw = 10e3; // FM audio lowpass
dsp.opt_dc = tharg->option_dc;
dsp.opt_dc = tharg->option_dc;
dsp.opt_cnt = tharg->option_cnt;
if ( dsp.sps < 8 ) {
fprintf(stderr, "note: sample rate low (%.1f sps)\n", dsp.sps);
//fprintf(stderr, "note: sample rate low (%.1f sps)\n", dsp.sps);
}
//headerlen = dsp.hdrlen;
@ -1208,7 +1214,7 @@ void *thd_m10(void *targs) { // pcm_t *pcm, double xlt_fq
k = init_buffers(&dsp);
if ( k < 0 ) {
fprintf(stderr, "error: init buffers\n");
return NULL;
goto exit_thread;
};
@ -1274,6 +1280,10 @@ void *thd_m10(void *targs) { // pcm_t *pcm, double xlt_fq
free_buffers(&dsp);
exit_thread:
reset_blockread(&dsp);
(dsp.thd)->used = 0;
return NULL;
}

Wyświetl plik

@ -40,7 +40,8 @@ typedef struct {
i8_t crc; // CRC check output
i8_t ecc; // Reed-Solomon ECC
i8_t sat; // GPS sat data
i8_t ptu; // PTU: temperature
i8_t ptu; // PTU: temperature humidity (pressure)
i8_t dwp; // PTU derived: dew point
i8_t inv;
i8_t aut;
i8_t jsn; // JSON output (auto_rx)
@ -77,7 +78,8 @@ typedef struct {
int std; int min; float sek;
double lat; double lon; double alt;
double vH; double vD; double vV;
float T; float RH;
float T; float RH; float TH;
float P; float RH2;
ui32_t crc;
ui8_t frame[FRAME_LEN];
ui8_t calibytes[51*16];
@ -89,6 +91,12 @@ typedef struct {
float ptu_co2[3]; // { -243.911 , 0.187654 , 8.2e-06 }
float ptu_calT2[3]; // calibration T2-Hum
float ptu_calH[2]; // calibration Hum
float ptu_mtxH[42];
float ptu_corHp[3];
float ptu_corHt[12];
float ptu_Cf1;
float ptu_Cf2;
float ptu_calP[25];
ui32_t freq; // freq/kHz (RS41)
int jsn_freq; // freq/kHz (SDR)
float batt; // battery voltage (V)
@ -203,6 +211,13 @@ static ui32_t u2(ui8_t *bytes) { // 16bit unsigned int
return bytes[0] | (bytes[1]<<8);
}
static int i2(ui8_t *bytes) { // 16bit signed int
//return (i16_t)u2(bytes);
int val = bytes[0] | (bytes[1]<<8);
if (val & 0x8000) val -= 0x10000;
return val;
}
/*
double r8(ui8_t *bytes) {
double val = 0;
@ -342,7 +357,7 @@ static int frametype(gpx_t *gpx) { // -4..+4: 0xF0 -> -4 , 0x0F -> +4
return ft;
}
static int get_FrameNb(gpx_t *gpx, int ofs) {
static int get_FrameNb(gpx_t *gpx, int crc, int ofs) {
int i;
unsigned byte;
ui8_t frnr_bytes[2];
@ -363,14 +378,14 @@ static int get_BattVolts(gpx_t *gpx, int ofs) {
int i;
unsigned byte;
ui8_t batt_bytes[2];
float batt_volts;
ui16_t batt_volts; // signed voltage?
for (i = 0; i < 2; i++) {
byte = gpx->frame[pos_BattVolts+ofs + i];
batt_bytes[i] = byte;
}
batt_volts = (float)(batt_bytes[0] + (batt_bytes[1] << 8));
// 2 bytes? V > 25.5 ?
batt_volts = batt_bytes[0]; // + (batt_bytes[1] << 8);
gpx->batt = batt_volts/10.0;
return 0;
@ -402,6 +417,8 @@ static int get_SondeID(gpx_t *gpx, int crc, int ofs) {
// don't reset gpx->frame[] !
gpx->T = -273.15;
gpx->RH = -1.0;
gpx->P = -1.0;
gpx->RH2 = -1.0;
// new ID:
memcpy(gpx->id, sondeid_bytes, 8);
gpx->id[8] = '\0';
@ -421,7 +438,7 @@ static int get_FrameConf(gpx_t *gpx, int ofs) {
err = crc;
err |= get_SondeID(gpx, crc, ofs);
err |= get_FrameNb(gpx, ofs);
err |= get_FrameNb(gpx, crc, ofs);
err |= get_BattVolts(gpx, ofs);
if (crc == 0) {
@ -440,6 +457,8 @@ static int get_FrameConf(gpx_t *gpx, int ofs) {
static int get_CalData(gpx_t *gpx) {
int j;
memcpy(&(gpx->ptu_Rf1), gpx->calibytes+61, 4); // 0x03*0x10+13
memcpy(&(gpx->ptu_Rf2), gpx->calibytes+65, 4); // 0x04*0x10+ 1
@ -450,6 +469,7 @@ static int get_CalData(gpx_t *gpx) {
memcpy(gpx->ptu_calT1+0, gpx->calibytes+89, 4); // 0x05*0x10+ 9
memcpy(gpx->ptu_calT1+1, gpx->calibytes+93, 4); // 0x05*0x10+13
memcpy(gpx->ptu_calT1+2, gpx->calibytes+97, 4); // 0x06*0x10+ 1
// ptu_calT1[3..6]
memcpy(gpx->ptu_calH+0, gpx->calibytes+117, 4); // 0x07*0x10+ 5
memcpy(gpx->ptu_calH+1, gpx->calibytes+121, 4); // 0x07*0x10+ 9
@ -461,82 +481,236 @@ static int get_CalData(gpx_t *gpx) {
memcpy(gpx->ptu_calT2+0, gpx->calibytes+305, 4); // 0x13*0x10+ 1
memcpy(gpx->ptu_calT2+1, gpx->calibytes+309, 4); // 0x13*0x10+ 5
memcpy(gpx->ptu_calT2+2, gpx->calibytes+313, 4); // 0x13*0x10+ 9
// ptu_calT2[3..6]
// cf. DF9DQ
memcpy(&(gpx->ptu_Cf1), gpx->calibytes+69, 4); // 0x04*0x10+ 5
memcpy(&(gpx->ptu_Cf2), gpx->calibytes+73, 4); // 0x04*0x10+ 9
for (j = 0; j < 42; j++) { // 0x07*0x10+13 = 0x07D = 125
memcpy(gpx->ptu_mtxH+j, gpx->calibytes+125+4*j, 4);
}
for (j = 0; j < 3; j++) { // 0x2A*0x10+6 = 0x2A6 = 678
memcpy(gpx->ptu_corHp+j, gpx->calibytes+678+4*j, 4);
}
for (j = 0; j < 12; j++) { // 0x2B*0x10+A = 0x2BA = 698
memcpy(gpx->ptu_corHt+j, gpx->calibytes+698+4*j, 4);
}
// cf. DF9DQ or stsst/RS-fork
memcpy(gpx->ptu_calP+0, gpx->calibytes+606, 4); // 0x25*0x10+14 = 0x25E
memcpy(gpx->ptu_calP+4, gpx->calibytes+610, 4); // ..
memcpy(gpx->ptu_calP+8, gpx->calibytes+614, 4);
memcpy(gpx->ptu_calP+12, gpx->calibytes+618, 4);
memcpy(gpx->ptu_calP+16, gpx->calibytes+622, 4);
memcpy(gpx->ptu_calP+20, gpx->calibytes+626, 4);
memcpy(gpx->ptu_calP+24, gpx->calibytes+630, 4);
memcpy(gpx->ptu_calP+1, gpx->calibytes+634, 4);
memcpy(gpx->ptu_calP+5, gpx->calibytes+638, 4);
memcpy(gpx->ptu_calP+9, gpx->calibytes+642, 4);
memcpy(gpx->ptu_calP+13, gpx->calibytes+646, 4);
memcpy(gpx->ptu_calP+2, gpx->calibytes+650, 4);
memcpy(gpx->ptu_calP+6, gpx->calibytes+654, 4);
memcpy(gpx->ptu_calP+10, gpx->calibytes+658, 4);
memcpy(gpx->ptu_calP+14, gpx->calibytes+662, 4);
memcpy(gpx->ptu_calP+3, gpx->calibytes+666, 4);
memcpy(gpx->ptu_calP+7, gpx->calibytes+670, 4); // ..
memcpy(gpx->ptu_calP+11, gpx->calibytes+674, 4); // 0x2A*0x10+ 2
return 0;
}
/*
static float get_Tc0(gpx_t *gpx, ui32_t f, ui32_t f1, ui32_t f2) {
// y = (f - f1) / (f2 - f1);
// y1 = (f - f1) / f2; // = (1 - f1/f2)*y
float a = 3.9083e-3, // Pt1000 platinum resistance
b = -5.775e-7,
c = -4.183e-12; // below 0C, else C=0
float *cal = gpx->ptu_calT1;
float Rb = (f1*gpx->ptu_Rf2-f2*gpx->ptu_Rf1)/(f2-f1), // ofs
Ra = f * (gpx->ptu_Rf2-gpx->ptu_Rf1)/(f2-f1) - Rb,
raw = Ra/1000.0,
g_r = 0.8024*cal[0] + 0.0176, // empirisch
r_o = 0.0705*cal[1] + 0.0011, // empirisch
r = raw * g_r + r_o,
t = (-a + sqrt(a*a + 4*b*(r-1)))/(2*b); // t>0: c=0
// R/R0 = 1 + at + bt^2 + c(t-100)t^3 , R0 = 1000 Ohm, t/Celsius
return t;
}
*/
// T_RH-sensor
static float get_TH(gpx_t *gpx, ui32_t f, ui32_t f1, ui32_t f2) {
float *p = gpx->ptu_co2;
float *c = gpx->ptu_calT2;
// temperature, platinum resistor
// T-sensor: gpx->ptu_co1 , gpx->ptu_calT1
// T_RH-sensor: gpx->ptu_co2 , gpx->ptu_calT2
static float get_T(gpx_t *gpx, ui32_t f, ui32_t f1, ui32_t f2, float *ptu_co, float *ptu_calT) {
float *p = ptu_co;
float *c = ptu_calT;
float g = (float)(f2-f1)/(gpx->ptu_Rf2-gpx->ptu_Rf1), // gain
Rb = (f1*gpx->ptu_Rf2-f2*gpx->ptu_Rf1)/(float)(f2-f1), // ofs
Rc = f/g - Rb,
//R = (Rc + c[1]) * c[0],
//T = p[0] + p[1]*R + p[2]*R*R;
R = Rc * c[0],
T = (p[0] + p[1]*R + p[2]*R*R + c[1])*(1.0 + c[2]);
return T;
}
// T-sensor, platinum resistor
static float get_Tc(gpx_t *gpx, ui32_t f, ui32_t f1, ui32_t f2) {
float *p = gpx->ptu_co1;
float *c = gpx->ptu_calT1;
float g = (float)(f2-f1)/(gpx->ptu_Rf2-gpx->ptu_Rf1), // gain
Rb = (f1*gpx->ptu_Rf2-f2*gpx->ptu_Rf1)/(float)(f2-f1), // ofs
Rc = f/g - Rb,
//R = (Rc + c[1]) * c[0],
//T = p[0] + p[1]*R + p[2]*R*R;
R = Rc * c[0],
T = (p[0] + p[1]*R + p[2]*R*R + c[1])*(1.0 + c[2]);
return T;
return T; // [Celsius]
}
// rel.hum., capacitor
// (data:) ftp://ftp-cdc.dwd.de/climate_environment/CDC/observations_germany/radiosondes/
// (diffAlt: Ellipsoid-Geoid)
// (note: humidity sensor has significant time-lag at low temperatures)
static float get_RH(gpx_t *gpx, ui32_t f, ui32_t f1, ui32_t f2, float T) {
float a0 = 7.5; // empirical
float a1 = 350.0/gpx->ptu_calH[0]; // empirical
float fh = (f-f1)/(float)(f2-f1);
float rh = 100.0 * ( a1*fh - a0 );
float T0 = 0.0, T1 = -25.0; // T/C
rh += T0 - T/5.5; // empir. temperature compensation
if (T < T1) rh *= 1.0 + (T1-T)/90.0; // empir. temperature compensation
float T0 = 0.0, T1 = -20.0, T2 = -40.0; // T/C v0.4
rh += T0 - T/5.5; // empir. temperature compensation
if (T < T1) rh *= 1.0 + (T1-T)/100.0; // empir. temperature compensation
if (T < T2) rh *= 1.0 + (T2-T)/120.0; // empir. temperature compensation
if (rh < 0.0) rh = 0.0;
if (rh > 100.0) rh = 100.0;
if (T < -273.0) rh = -1.0;
return rh;
}
static int get_PTU(gpx_t *gpx, int ofs, int pck) {
// ---------------------------------------------------------------------------------------
//
// cf. github DF9DQ
// water vapor saturation pressure (Hyland and Wexler)
static float vaporSatP(float Tc) {
double T = Tc + 273.15;
// Apply some correction magic
// T = -0.4931358 + (1.0 + 4.61e-3) * T - 1.3746454e-5 * T*T + 1.2743214e-8 * T*T*T;
// H+W equation
double p = expf(-5800.2206 / T
+1.3914993
+6.5459673 * log(T)
-4.8640239e-2 * T
+4.1764768e-5 * T*T
-1.4452093e-8 * T*T*T
);
return (float)p; // [Pa]
}
// cf. github DF9DQ
static float get_RH2adv(gpx_t *gpx, ui32_t f, ui32_t f1, ui32_t f2, float T, float TH, float P) {
float rh = 0.0;
float cfh = (f-f1)/(float)(f2-f1);
float cap = gpx->ptu_Cf1+(gpx->ptu_Cf2-gpx->ptu_Cf1)*cfh;
double Cp = (cap / gpx->ptu_calH[0] - 1.0) * gpx->ptu_calH[1];
double Trh_20_180 = (TH - 20.0) / 180.0;
double _rh = 0.0;
double aj = 1.0;
double bk = 1.0, b[6];
int j, k;
bk = 1.0;
for (k = 0; k < 6; k++) {
b[k] = bk; // Tp^k
bk *= Trh_20_180;
}
if (P > 0.0) // in particular if P<200hPa , T<-40
{
double _p = P / 1000.0; // bar
double _cpj = 1.0;
double corrCp = 0.0;
double bt, bp[3];
for (j = 0; j < 3; j++) {
bp[j] = gpx->ptu_corHp[j] * (_p/(1.0 + gpx->ptu_corHp[j]*_p) - _cpj/(1.0 + gpx->ptu_corHp[j]));
_cpj *= Cp; // Cp^j
}
corrCp = 0.0;
for (j = 0; j < 3; j++) {
bt = 0.0;
for (k = 0; k < 4; k++) {
bt += gpx->ptu_corHt[4*j+k] * b[k];
}
corrCp += bp[j] * bt;
}
Cp -= corrCp;
}
aj = 1.0;
for (j = 0; j < 7; j++) {
for (k = 0; k < 6; k++) {
_rh += aj * b[k] * gpx->ptu_mtxH[6*j+k];
}
aj *= Cp;
}
if ( P <= 0.0 ) { // empirical correction
float T2 = -40;
if (T < T2) { _rh += (T-T2)/12.0; }
}
rh = _rh * vaporSatP(TH)/vaporSatP(T);
if (rh < 0.0) rh = 0.0;
if (rh > 100.0) rh = 100.0;
return rh;
}
//
// cf. github DF9DQ or stsst/RS-fork
static float get_P(gpx_t *gpx, ui32_t f, ui32_t f1, ui32_t f2, int fx)
{
double p = 0.0;
double a0, a1, a0j, a1k;
int j, k;
if (f1 == f2 || f1 == f) return 0.0;
a0 = gpx->ptu_calP[24] / ((float)(f - f1) / (float)(f2 - f1));
a1 = fx * 0.01;
a0j = 1.0;
for (j = 0; j < 6; j++) {
a1k = 1.0;
for (k = 0; k < 4; k++) {
p += a0j * a1k * gpx->ptu_calP[j*4+k];
a1k *= a1;
}
a0j *= a0;
}
return (float)p;
}
// ---------------------------------------------------------------------------------------
//
// barometric formula https://en.wikipedia.org/wiki/Barometric_formula
static float Ph(float h) {
double Pb, Tb, Lb, hb;
//double RgM = 8.31446/(9.80665*0.0289644);
double gMR = 9.80665*0.0289644/8.31446;
float P = 0.0;
if (h > 32000.0) { //P < 8.6802
Pb = 8.6802;
Tb = 228.65;
Lb = 0.0028;
hb = 32000.0;
}
else if (h > 20000.0) { // P < 54.7489 (&& P >= 8.6802)
Pb = 54.7489;
Tb = 216.65;
Lb = 0.001;
hb = 20000.0;
}
else if (h > 11000.0) { // P < 226.321 (&& P >= 54.7489)
Pb = 226.321;
Tb = 216.65;
Lb = 0.0;
hb = 11000.0;
}
else { // P >= 226.321
Pb = 1013.25;
Tb = 288.15;
Lb = -0.0065;
hb = 0.0;
}
//if (Lb == 0.0) altP = -RgM*Tb * log(P/Pb) + hb;
//else altP = Tb/Lb * (pow(P/Pb, -RgM*Lb)-1.0) + hb;
if (Lb == 0.0) P = Pb * exp( -gMR*(h-hb)/Tb );
else P = Pb * pow( 1.0+Lb*(h-hb)/Tb , -gMR/Lb);
return P;
}
static int get_PTU(gpx_t *gpx, int ofs, int pck, int valid_alt) {
int err=0, i;
int bR, bc1, bT1,
bc2, bT2;
int bH;
int bH2;
int bP;
ui32_t meas[12];
float Tc = -273.15;
float TH = -273.15;
float RH = -1.0;
float RH2 = -1.0;
float P = -1.0;
get_CalData(gpx);
@ -558,21 +732,53 @@ static int get_PTU(gpx_t *gpx, int ofs, int pck) {
bT2 = gpx->calfrchk[0x13];
bH = gpx->calfrchk[0x07];
bH2 = gpx->calfrchk[0x07] && gpx->calfrchk[0x08]
&& gpx->calfrchk[0x09] && gpx->calfrchk[0x0A]
&& gpx->calfrchk[0x0B] && gpx->calfrchk[0x0C]
&& gpx->calfrchk[0x0D] && gpx->calfrchk[0x0E]
&& gpx->calfrchk[0x0F] && gpx->calfrchk[0x10]
&& gpx->calfrchk[0x11] && gpx->calfrchk[0x12]
&& gpx->calfrchk[0x2A] && gpx->calfrchk[0x2B]
&& gpx->calfrchk[0x2C] && gpx->calfrchk[0x2D]
&& gpx->calfrchk[0x2E];
bP = gpx->calfrchk[0x21] && gpx->calibytes[0x21F] == 'P'
&& gpx->calfrchk[0x25] && gpx->calfrchk[0x26]
&& gpx->calfrchk[0x27] && gpx->calfrchk[0x28]
&& gpx->calfrchk[0x29] && gpx->calfrchk[0x2A];
if (bR && bc1 && bT1) {
Tc = get_Tc(gpx, meas[0], meas[1], meas[2]);
//Tc0 = get_Tc0(gpx, meas[0], meas[1], meas[2]);
Tc = get_T(gpx, meas[0], meas[1], meas[2], gpx->ptu_co1, gpx->ptu_calT1);
}
gpx->T = Tc;
if (bR && bc2 && bT2) {
TH = get_TH(gpx, meas[6], meas[7], meas[8]);
TH = get_T(gpx, meas[6], meas[7], meas[8], gpx->ptu_co2, gpx->ptu_calT2);
}
gpx->TH = TH;
if (bH) {
if (bH && Tc > -273.0) {
RH = get_RH(gpx, meas[3], meas[4], meas[5], Tc); // TH, TH-Tc (sensorT - T)
}
gpx->RH = RH;
// cf. DF9DQ, stsst/RS-fork
if (bP) {
P = get_P(gpx, meas[9], meas[10], meas[11], i2(gpx->frame+pos_PTU+ofs+2+38));
}
gpx->P = P;
if (gpx->option.ptu == 2) {
float _P = -1.0;
if (bP) _P = P;
else { // approx
if (valid_alt > 0) _P = Ph(gpx->alt);
}
if (bH && bH2 && Tc > -273.0 && TH > -273.0) {
RH2 = get_RH2adv(gpx, meas[3], meas[4], meas[5], Tc, TH, _P);
}
}
gpx->RH2 = RH2;
if (gpx->option.vbs == 4 && (gpx->crc & (crc_PTU | crc_GPS3))==0)
{
@ -585,7 +791,7 @@ static int get_PTU(gpx_t *gpx, int ofs, int pck) {
printf("3: %8d %8d %8d", meas[6], meas[7], meas[8]);
printf(" # ");
//if (Tc > -273.0 && RH > -0.5)
if (0 && Tc > -273.0 && RH > -0.5)
{
printf(" ");
printf(" Tc:%.2f ", Tc);
@ -696,6 +902,7 @@ static int get_GPS2(gpx_t *gpx, int ofs) {
return err;
}
// WGS84/GRS80 Ellipsoid
#define EARTH_a 6378137.0
#define EARTH_b 6356752.31424518
#define EARTH_a2_b2 (EARTH_a*EARTH_a - EARTH_b*EARTH_b)
@ -889,10 +1096,10 @@ static int get_Calconf(gpx_t *gpx, int out, int ofs) {
byte = gpx->frame[pos_CalData+ofs+1+i];
fprintf(stdout, "%02x ", byte);
}
/*
/*
if (err == 0) fprintf(stdout, "[OK]");
else fprintf(stdout, "[NO]");
*/
*/
fprintf(stdout, " ");
}
@ -1001,7 +1208,7 @@ static int rs41_ecc(gpx_t *gpx, int frmlen) {
if (gpx->option.ecc == 2 && (errors1 < 0 || errors2 < 0))
{ // 2nd pass
{ // 2nd pass: set packet-IDs
gpx->frame[pos_FRAME] = (pck_FRAME>>8)&0xFF; gpx->frame[pos_FRAME+1] = pck_FRAME&0xFF;
gpx->frame[pos_PTU] = (pck_PTU >>8)&0xFF; gpx->frame[pos_PTU +1] = pck_PTU &0xFF;
gpx->frame[pos_GPS1] = (pck_GPS1 >>8)&0xFF; gpx->frame[pos_GPS1 +1] = pck_GPS1 &0xFF;
@ -1073,12 +1280,32 @@ static int prn_frm(gpx_t *gpx) {
static int prn_ptu(gpx_t *gpx) {
fprintf(stdout, " ");
if (gpx->T > -273.0) fprintf(stdout, " T=%.1fC ", gpx->T);
if (gpx->RH > -0.5) fprintf(stdout, " RH=%.0f%% ", gpx->RH);
if (gpx->RH > -0.5 && gpx->option.ptu != 2) fprintf(stdout, " _RH=%.0f%% ", gpx->RH);
if (gpx->P > 0.0) {
if (gpx->P < 100.0) fprintf(stdout, " P=%.2fhPa ", gpx->P);
else fprintf(stdout, " P=%.1fhPa ", gpx->P);
}
if (gpx->option.ptu == 2) {
if (gpx->RH2 > -0.5) fprintf(stdout, " RH2=%.0f%% ", gpx->RH2);
}
// dew point
if (gpx->option.dwp)
{
float rh = gpx->RH;
float Td = -273.15f; // dew point Td
if (gpx->option.ptu == 2) rh = gpx->RH2;
if (rh > 0.0f && gpx->T > -273.0f) {
float gamma = logf(rh / 100.0f) + (17.625f * gpx->T / (243.04f + gpx->T));
Td = 243.04f * gamma / (17.625f - gamma);
fprintf(stdout, " Td=%.1fC ", Td);
}
}
return 0;
}
static int prn_gpstime(gpx_t *gpx) {
Gps2Date(gpx);
//Gps2Date(gpx);
fprintf(stdout, "%s ", weekday[gpx->wday]);
fprintf(stdout, "%04d-%02d-%02d %02d:%02d:%06.3f",
gpx->jahr, gpx->monat, gpx->tag, gpx->std, gpx->min, gpx->sek);
@ -1154,7 +1381,7 @@ static int prn_sat3(gpx_t *gpx, int ofs) {
pDOP = gpx->frame[pos_pDOP+ofs]/10.0; if (gpx->frame[pos_pDOP+ofs] == 0xFF) pDOP = -1.0;
fprintf(stdout, "numSatsFix: %2d sAcc: %.1f pDOP: %.1f\n", numSV, sAcc, pDOP);
/*
/*
fprintf(stdout, "CRC: ");
fprintf(stdout, " %04X", pck_GPS1);
if (check_CRC(gpx, pos_GPS1+ofs, pck_GPS1)==0) fprintf(stdout, "[OK]"); else fprintf(stdout, "[NO]");
@ -1167,12 +1394,12 @@ static int prn_sat3(gpx_t *gpx, int ofs) {
//fprintf(stdout, "[%+d]", check_CRC(gpx, pos_GPS3, pck_GPS3));
fprintf(stdout, "\n");
*/
*/
return 0;
}
static int print_position(gpx_t *gpx, int ec) {
int i, j;
int i;
int err, err0, err1, err2, err3;
//int output, out_mask;
int encrypted = 0;
@ -1180,6 +1407,7 @@ static int print_position(gpx_t *gpx, int ec) {
int out = 1;
int sat = 0;
int pos_aux = 0, cnt_aux = 0;
int ofs_ptu = 0, pck_ptu = 0;
int ret = 0;
//gpx->out = 0;
@ -1234,10 +1462,11 @@ static int print_position(gpx_t *gpx, int ec) {
break;
case pck_PTU: // 0x7A2A
ofs = pos - pos_PTU;
err0 = get_PTU(gpx, ofs, pck_PTU);
if ( 0 && !err0 && gpx->option.ptu ) {
prn_ptu(gpx);
ofs_ptu = pos - pos_PTU;
pck_ptu = pck_PTU;
if ( 0 && gpx->option.ptu ) {
//err0 = get_PTU(gpx, ofs_ptu, pck_ptu);
// if (!err0) prn_ptu(gpx);
}
break;
@ -1245,6 +1474,7 @@ static int print_position(gpx_t *gpx, int ec) {
ofs = pos - pos_GPS1;
err1 = get_GPS1(gpx, ofs);
if ( !err1 ) {
Gps2Date(gpx);
if (out) prn_gpstime(gpx);
if (sat) prn_sat1(gpx, ofs);
}
@ -1268,8 +1498,11 @@ static int print_position(gpx_t *gpx, int ec) {
break;
case pck_SGM_xTU: // 0x7F1B
ofs = pos - pos_PTU;
err0 = get_PTU(gpx, ofs, pck);
ofs_ptu = pos - pos_PTU;
pck_ptu = pck;
if ( 0 ) {
//err0 = get_PTU(gpx, ofs_ptu, pck_ptu);
}
break;
case pck_SGM_CRYPT: // 0x80A7
@ -1301,9 +1534,11 @@ static int print_position(gpx_t *gpx, int ec) {
if ( pos > frm_end ) // end of (sub)frame
{
if (gpx->option.ptu && out && !sat && !err0 && !encrypted) {
prn_ptu(gpx);
if (gpx->option.ptu && !sat && !encrypted && pck_ptu > 0) {
err0 = get_PTU(gpx, ofs_ptu, pck_ptu, !err3);
if (!err0 && out) prn_ptu(gpx);
}
pck_ptu = 0;
get_Calconf(gpx, out, ofs_cal);
@ -1325,11 +1560,18 @@ static int print_position(gpx_t *gpx, int ec) {
fprintf(stdout, "{ \"type\": \"%s\"", "RS41");
fprintf(stdout, ", \"frame\": %d, \"id\": \"%s\", \"datetime\": \"%04d-%02d-%02dT%02d:%02d:%06.3fZ\", \"lat\": %.5f, \"lon\": %.5f, \"alt\": %.5f, \"vel_h\": %.5f, \"heading\": %.5f, \"vel_v\": %.5f, \"sats\": %d, \"bt\": %d, \"batt\": %.2f",
gpx->frnr, gpx->id, gpx->jahr, gpx->monat, gpx->tag, gpx->std, gpx->min, gpx->sek, gpx->lat, gpx->lon, gpx->alt, gpx->vH, gpx->vD, gpx->vV, gpx->numSV, gpx->conf_cd, gpx->batt );
if (gpx->option.ptu && !err0 && gpx->T > -273.0) {
fprintf(stdout, ", \"temp\": %.1f", gpx->T );
}
if (gpx->option.ptu && !err0 && gpx->RH > -0.5) {
fprintf(stdout, ", \"humidity\": %.1f", gpx->RH );
if (gpx->option.ptu && !err0) {
float _RH = gpx->RH;
if (gpx->option.ptu == 2) _RH = gpx->RH2;
if (gpx->T > -273.0) {
fprintf(stdout, ", \"temp\": %.1f", gpx->T );
}
if (_RH > -0.5) {
fprintf(stdout, ", \"humidity\": %.1f", _RH );
}
if (gpx->P > 0.0) {
fprintf(stdout, ", \"pressure\": %.2f", gpx->P );
}
}
if (gpx->aux) { // <=> gpx->xdata[0]!='\0'
fprintf(stdout, ", \"aux\": \"%s\"", gpx->xdata );
@ -1374,13 +1616,16 @@ static int print_position(gpx_t *gpx, int ec) {
ofs = 0;
if (pck < 0x8000) {
err0 = get_PTU(gpx, 0, pck);
//err0 = get_PTU(gpx, 0, pck, 0);
if (pck == pck_PTU) ofs = 0;
else if (pck == pck_SGM_xTU) ofs = 0x1B-0x2A;
err1 = get_GPS1(gpx, ofs);
err2 = get_GPS2(gpx, ofs);
err3 = get_GPS3(gpx, ofs);
if (!err1) Gps2Date(gpx);
err0 = get_PTU(gpx, 0, pck, !err3);
if (out) {
@ -1464,11 +1709,16 @@ static void print_frame(gpx_t *gpx, int len, dsp_t *dsp) {
fprintf(stdout, "\n");
}
else {
pthread_mutex_lock( dsp->thd.mutex );
fprintf(stdout, "<%d> ", dsp->thd.tn);
pthread_mutex_lock( dsp->thd->mutex );
//fprintf(stdout, "<%d> ", dsp->thd->tn);
fprintf(stdout, "<%d: ", dsp->thd->tn);
fprintf(stdout, "s=%+.4f, ", dsp->mv);
fprintf(stdout, "f=%+.4f", -dsp->thd->xlt_fq);
if (dsp->opt_dc) fprintf(stdout, "%+.6f", dsp->Df/(double)dsp->sr);
fprintf(stdout, "> ");
ret = print_position(gpx, ec);
if (ret==0) fprintf(stdout, "\n");
pthread_mutex_unlock( dsp->thd.mutex );
pthread_mutex_unlock( dsp->thd->mutex );
}
}
@ -1518,7 +1768,7 @@ void *thd_rs41(void *targs) { // pcm_t *pcm, double xlt_fq
// init gpx
gpx.option.vbs = 1;
gpx.option.ptu = 1;
gpx.option.ptu = 2;
gpx.option.aut = 1;
gpx.option.jsn = tharg->option_jsn;
@ -1546,7 +1796,7 @@ void *thd_rs41(void *targs) { // pcm_t *pcm, double xlt_fq
dsp.dectaps = pcm->dectaps;
dsp.decM = pcm->decM;
dsp.thd = tharg->thd;
dsp.thd = &(tharg->thd);
dsp.bps = pcm->bps;
dsp.nch = pcm->nch;
@ -1564,17 +1814,18 @@ void *thd_rs41(void *targs) { // pcm_t *pcm, double xlt_fq
dsp.opt_lp = 1;
dsp.lpIQ_bw = 8e3; // IF lowpass bandwidth
dsp.lpFM_bw = 6e3; // FM audio lowpass
dsp.opt_dc = tharg->option_dc;
dsp.opt_dc = tharg->option_dc;
dsp.opt_cnt = tharg->option_cnt;
if ( dsp.sps < 8 ) {
fprintf(stderr, "note: sample rate low (%.1f sps)\n", dsp.sps);
//fprintf(stderr, "note: sample rate low (%.1f sps)\n", dsp.sps);
}
k = init_buffers(&dsp); // BT=0.5 (IQ-Int: BT > 0.5 ?)
if ( k < 0 ) {
fprintf(stderr, "error: init buffers\n");
return NULL;
goto exit_thread;
};
//if (option_iq: 2,3) bitofs += 1; // FM: +1 , IQ: +2, IQ5: +1
@ -1634,6 +1885,9 @@ void *thd_rs41(void *targs) { // pcm_t *pcm, double xlt_fq
free_buffers(&dsp);
exit_thread:
reset_blockread(&dsp);
(dsp.thd)->used = 0;
return NULL;
}

Wyświetl plik

@ -1,5 +1,6 @@
/*
gcc -O2 -c demod_base.c
gcc -O2 -c bch_ecc_mod.c
gcc -O2 -c rs41base.c
@ -10,13 +11,26 @@ gcc -O2 rs_multi.c demod_base.o bch_ecc_mod.o rs41base.o dfm09base.o m10base.o l
./a.out --rs41 <fq0> --dfm <fq1> --m10 <fq2> baseband_IQ.wav
-0.5 < fq < 0.5 , fq=freq/sr
e.g.
[terminal 1]
./a.out --dc --fifo rsfifo --m10 -0.28034 --lms 0.12636 baseband_IQ.wav
[terminal 2]
echo "m10 -0.36" > rsfifo
echo "-1" > rsfifo
echo "lms 0.02428" > rsfifo
*/
#include <stdio.h>
#include <sys/stat.h> // mkfifo()
#include <unistd.h> // open(),close(),unlink()
#include <fcntl.h> // O_RDONLY //... setmode()/cygwin
#ifdef CYGWIN
#include <fcntl.h> // cygwin: _setmode()
//#include <fcntl.h> // cygwin: _setmode()
#include <io.h>
#endif
@ -29,8 +43,8 @@ static pthread_cond_t cond = PTHREAD_COND_INITIALIZER;
static float complex *block_decMB;
int rbf1; // extern in demod_base.c
extern int rbf1; // demod_base.c
extern int bufeof; // demod_base.c
void *thd_rs41(void *);
void *thd_dfm09(void *);
@ -85,6 +99,7 @@ static int pcm_dec_init(pcm_t *p) {
return 0;
}
#define FIFOBUF_LEN 20
int main(int argc, char **argv) {
@ -98,7 +113,15 @@ int main(int argc, char **argv) {
int option_pcmraw = 0,
option_jsn = 0,
option_dc = 0,
option_min = 0;
option_min = 0,
option_cont = 0;
// FIFO
int option_fifo = 0;
char *rs_fifo = NULL;
int fd = -1;
char fifo_buf[FIFOBUF_LEN];
int th_used = 0;
#ifdef CYGWIN
_setmode(fileno(stdin), _O_BINARY); // _fileno(stdin)
@ -115,8 +138,7 @@ int main(int argc, char **argv) {
if (strcmp(*argv, "--rs41") == 0) {
double fq = 0.0;
++argv;
if (*argv) fq = atof(*argv);
else return -1;
if (*argv) fq = atof(*argv); else return -1;
if (fq < -0.5) fq = -0.5;
if (fq > 0.5) fq = 0.5;
if (xlt_cnt < MAX_FQ) {
@ -128,8 +150,7 @@ int main(int argc, char **argv) {
else if (strcmp(*argv, "--dfm") == 0) {
double fq = 0.0;
++argv;
if (*argv) fq = atof(*argv);
else return -1;
if (*argv) fq = atof(*argv); else return -1;
if (fq < -0.5) fq = -0.5;
if (fq > 0.5) fq = 0.5;
if (xlt_cnt < MAX_FQ) {
@ -141,8 +162,7 @@ int main(int argc, char **argv) {
else if (strcmp(*argv, "--m10") == 0) {
double fq = 0.0;
++argv;
if (*argv) fq = atof(*argv);
else return -1;
if (*argv) fq = atof(*argv); else return -1;
if (fq < -0.5) fq = -0.5;
if (fq > 0.5) fq = 0.5;
if (xlt_cnt < MAX_FQ) {
@ -154,8 +174,7 @@ int main(int argc, char **argv) {
else if (strcmp(*argv, "--lms") == 0) {
double fq = 0.0;
++argv;
if (*argv) fq = atof(*argv);
else return -1;
if (*argv) fq = atof(*argv); else return -1;
if (fq < -0.5) fq = -0.5;
if (fq > 0.5) fq = 0.5;
if (xlt_cnt < MAX_FQ) {
@ -180,6 +199,18 @@ int main(int argc, char **argv) {
else if (strcmp(*argv, "--min") == 0) {
option_min = 1;
}
else if ( (strcmp(*argv, "-c") == 0) || (strcmp(*argv, "--cnt") == 0) ) {
option_cont = 1;
}
else if (strcmp(*argv, "--fifo") == 0) {
++argv;
if (*argv) rs_fifo = *argv; else return -1;
unlink(rs_fifo);
if (mkfifo(rs_fifo, 0666) < 0) { // evtl. mode : S_IWUSR | S_IRUSR
//fprintf(stderr, "error open %s\n", rs_fifo); // exit, if exists ?
}
option_fifo = 1;
}
else if (strcmp(*argv, "-") == 0) {
int sample_rate = 0, bits_sample = 0, channels = 0;
++argv;
@ -199,7 +230,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;
@ -229,7 +260,8 @@ int main(int argc, char **argv) {
block_decMB = calloc(pcm.decM*blk_sz+1, sizeof(float complex)); if (block_decMB == NULL) return -1;
thargs_t tharg[xlt_cnt];
thargs_t tharg[MAX_FQ]; // xlt_cnt<=MAX_FQ
for (k = 0; k < MAX_FQ; k++) tharg[k].thd.used = 0;
for (k = 0; k < xlt_cnt; k++) {
tharg[k].thd.tn = k;
@ -249,8 +281,10 @@ int main(int argc, char **argv) {
tharg[k].option_jsn = option_jsn;
tharg[k].option_dc = option_dc;
tharg[k].option_cnt = option_cont;
rbf1 |= tharg[k].thd.tn_bit;
tharg[k].thd.used = 1;
}
for (k = 0; k < xlt_cnt; k++) {
@ -258,16 +292,129 @@ int main(int argc, char **argv) {
}
// FIFO
//
if (option_fifo)
{
fd = open(rs_fifo, O_RDONLY | O_NONBLOCK); //fcntl.h
if (fd < 0) {
fprintf(stderr, "error open %s\n", rs_fifo);
return -1;
}
while ( !bufeof ) {
int l = 0;
memset(fifo_buf, 0, FIFOBUF_LEN);
th_used = 0;
for (k = 0; k < MAX_FQ; k++) th_used += tharg[k].thd.used;
if (th_used == 0) break;
l = read(fd, fifo_buf, FIFOBUF_LEN);
if ( l > 1 ) {
void *rstype = NULL;
char *fifo_fq = fifo_buf;
while (l > 1 && fifo_buf[l-1] < 0x20) l--;
fifo_buf[l] = '\0'; // remove \n, terminate string
if (strncmp(fifo_buf, "rs41", 4) == 0) {
fifo_fq = fifo_buf + 4;
rstype = thd_rs41;
}
else if (strncmp(fifo_buf, "dfm", 3) == 0) {
fifo_fq = fifo_buf + 3;
rstype = thd_dfm09;
}
else if (strncmp(fifo_buf, "m10", 3) == 0) {
fifo_fq = fifo_buf + 3;
rstype = thd_m10;
}
else if (strncmp(fifo_buf, "lms", 3) == 0) {
fifo_fq = fifo_buf + 3;
rstype = thd_lms6X;
}
else {
if (fifo_buf[0] == '-') { // -<n> : close <n>
int num = atoi(fifo_buf+1);
if (num >= 0 && num < MAX_FQ) {
if (tharg[num].thd.used) {
tharg[num].thd.used = 0;
}
}
}
continue;
}
double fq = 0.0;
if (fifo_fq) fq = atof(fifo_fq);
else return -1;
if (fq < -0.5) fq = -0.5;
if (fq > 0.5) fq = 0.5;
// find slot
for (k = 0; k < MAX_FQ; k++) {
if (tharg[k].thd.used == 0) break;
}
if (k < MAX_FQ) {
double base_fq = fq;
tharg[k].thd.tn = k;
tharg[k].thd.tn_bit = (1<<k);
tharg[k].thd.mutex = &mutex;
tharg[k].thd.cond = &cond;
//tharg[k].thd.lock = &lock;
tharg[k].thd.blk = block_decMB;
tharg[k].thd.xlt_fq = -base_fq;
if (cfreq > 0) {
int fq_kHz = (cfreq - tharg[k].thd.xlt_fq*pcm.sr_base + 500)/1e3;
tharg[k].jsn_freq = fq_kHz;
}
tharg[k].pcm = pcm;
tharg[k].option_jsn = option_jsn;
tharg[k].option_dc = option_dc;
rbf1 |= tharg[k].thd.tn_bit;
tharg[k].thd.used = 1;
pthread_create(&tharg[k].thd.tid, NULL, rstype, &tharg[k]);
pthread_mutex_lock( &mutex );
fprintf(stdout, "<%d: add f=%+.4f>\n", k, base_fq);
pthread_mutex_unlock( &mutex );
k++;
if (k > xlt_cnt) xlt_cnt = k;
for (k = 0; k < xlt_cnt; k++) {
tharg[k].thd.max_fq = xlt_cnt;
}
}
}
sleep(1);
}
}
for (k = 0; k < xlt_cnt; k++) {
pthread_join(tharg[k].thd.tid, NULL);
}
th_used = 1;
while (th_used) {
th_used = 0;
for (k = 0; k < MAX_FQ; k++) th_used += tharg[k].thd.used;
}
if (block_decMB) { free(block_decMB); block_decMB = NULL; }
decimate_free();
fclose(fp);
if (fd >= 0) { // FIFO
close(fd);
unlink(rs_fifo);
}
return 0;
}