From 4230198d91211716e448eeaf2e1cc178e88aee86 Mon Sep 17 00:00:00 2001 From: ha7ilm Date: Sun, 20 Mar 2016 16:41:37 +0100 Subject: [PATCH] Added a squelch --- README.md | 54 +++++----- csdr.c | 301 ++++++++++++++++++++++++++++++++---------------------- libcsdr.c | 210 ++++++++++++++++++++----------------- libcsdr.h | 26 ++--- 4 files changed, 335 insertions(+), 256 deletions(-) diff --git a/README.md b/README.md index b01e807..4f32598 100644 --- a/README.md +++ b/README.md @@ -36,7 +36,7 @@ Usage by example - Baseband I/Q signal is coming from an RTL-SDR USB dongle, with a center frequency of `-f 104300000` Hz, a sampling rate of `-s 240000` samples per second. - The `rtl_sdr` tool outputs an unsigned 8-bit I/Q signal (one byte of I sample and one byte of Q coming after each other), but `libcsdr` DSP routines internally use floating point data type, so we convert the data stream of `unsigned char` to `float` by `csdr convert_u8_f`. - We want to listen one radio station at the frequency `-f 89500000` Hz (89.5 MHz). -- No other radio station is within the sampled bandwidth, so we send the signal directly to the demodulator. (This is an easy, but not perfect solution as the anti-aliasing filter at RTL-SDR DDC is too short.) +- No other radio station is within the sampled bandwidth, so we send the signal directly to the demodulator. (This is an easy, but not perfect solution as the anti-aliasing filter at RTL-SDR DDC is too short.) - After FM demodulation we decimate the signal by a factor of 5 to match the rate of the audio card (240000 / 5 = 48000). - A de-emphasis filter is used, because pre-emphasis is applied at the transmitter to compensate noise at higher frequencies. The time constant for de-emphasis for FM broadcasting in Europe is 50 microseconds (hence the `50e-6`). - Also, `mplayer` cannot play floating point audio, so we convert our signal to a stream of 16-bit integers. @@ -58,12 +58,12 @@ Sample rates look like this: *Note:* there is an example shell script that does this for you (without the unnecessary shift operation). If you just want to listen to FM radio, type: - csdr-fm 89.5 20 + csdr-fm 89.5 20 The first parameter is the frequency in MHz, and the second optional parameter is the RTL-SDR tuner gain in dB. ### Demodulate NFM - + rtl_sdr -s 2400000 -f 145000000 -g 20 - | csdr convert_u8_f | csdr shift_addition_cc `python -c "print float(145000000-145350000)/2400000"` | csdr fir_decimate_cc 50 0.005 HAMMING | csdr fmdemod_quadri_cf | csdr limit_ff | csdr deemphasis_nfm_ff 48000 | csdr fastagc_ff | csdr convert_f_i16 | mplayer -cache 1024 -quiet -rawaudio samplesize=2:channels=1:rate=48000 -demuxer rawaudio - - Note that the decimation factor is higher (we want to select a ~25 kHz channel). @@ -74,7 +74,7 @@ The first parameter is the frequency in MHz, and the second optional parameter i rtl_sdr -s 2400000 -f 145000000 -g 20 - | csdr convert_u8_f | csdr shift_addition_cc `python -c "print float(145000000-144400000)/2400000"` | csdr fir_decimate_cc 50 0.005 HAMMING | csdr amdemod_cf | csdr fastdcblock_ff | csdr agc_ff | csdr limit_ff | csdr convert_f_i16 | mplayer -cache 1024 -quiet -rawaudio samplesize=2:channels=1:rate=48000 -demuxer rawaudio - -- `amdemod_cf` is used as demodulator. +- `amdemod_cf` is used as demodulator. - `agc_ff` should be used for AM and SSB. ### Design FIR band-pass filter (with complex taps) @@ -89,14 +89,14 @@ The first parameter is the frequency in MHz, and the second optional parameter i rtl_sdr -s 2400000 -f 145000000 -g 20 - | csdr convert_u8_f | csdr shift_addition_cc `python -c "print float(145000000-144400000)/2400000"` | csdr fir_decimate_cc 50 0.005 HAMMING | csdr bandpass_fir_fft_cc 0 0.1 0.05 | csdr realpart_cf | csdr agc_ff | csdr limit_ff | csdr convert_f_i16 | mplayer -cache 1024 -quiet -rawaudio samplesize=2:channels=1:rate=48000 -demuxer rawaudio - -- It is a modified Weaver-demodulator. The complex FIR filter removes the lower sideband and lets only the upper pass (USB). If you want to demodulate LSB, change `bandpass_fir_fft_cc 0 0.05` to `bandpass_fir_fft_cc -0.05 0`. +- It is a modified Weaver-demodulator. The complex FIR filter removes the lower sideband and lets only the upper pass (USB). If you want to demodulate LSB, change `bandpass_fir_fft_cc 0 0.05` to `bandpass_fir_fft_cc -0.05 0`. ### Draw FFT rtl_sdr -s 2400000 -f 104300000 -g 20 - | csdr convert_u8_f | csdr fft_cc 1024 1200000 HAMMING --octave | octave -i > /dev/null - We calculate the Fast Fourier Transform by `csdr fft_cc` on the first 1024 samples of every block of 1200000 complex samples coming after each other. (We calculate FFT from 1024 samples and then skip 1200000-1024=1198976 samples. This way we will calculate FFT two times every second.) -- The window used for FFT is the Hamming window, and the output consists of commands that can be directly interpreted by GNU Octave which plots us the spectrum. +- The window used for FFT is the Hamming window, and the output consists of commands that can be directly interpreted by GNU Octave which plots us the spectrum. Usage ----- @@ -107,7 +107,7 @@ Function name endings found in *libcsdr* mean the input and output data types of Data types are noted as it follows: - `f` is `float` (single percision) -- `c` is `complexf` (two single precision floating point values in a struct) +- `c` is `complexf` (two single precision floating point values in a struct) - `u8` is `unsigned char` of 1 byte/8 bits (e. g. the output of `rtl_sdr` is of `u8`) - `s16` is `signed short` of 2 bytes/16 bits (e. g. sound card input is usually `s16`) @@ -117,14 +117,14 @@ Functions usually end as: - `_cf` complex input, float output - `_cc` complex input, complex output -Regarding *csdr*, it can convert a real/complex stream from one data format to another, to interface it with other SDR tools and the sound card. -The following commands are available: +Regarding *csdr*, it can convert a real/complex stream from one data format to another, to interface it with other SDR tools and the sound card. +The following commands are available: -- `csdr convert_u8_f` -- `csdr convert_f_u8` +- `csdr convert_u8_f` +- `csdr convert_f_u8` - `csdr convert_s8_f` - `csdr convert_f_s8` -- `csdr convert_s16_f` +- `csdr convert_s16_f` - `csdr convert_f_s16` How to interpret: `csdr convert__` @@ -169,7 +169,7 @@ The `csdr` process just exits with 0. yes_f [buf_times] -It outputs continously the `to_repeat` float number. +It outputs continously the `to_repeat` float number. If `buf_times` is not given, it never stops. Else, after outputing `buf_times` number of buffers (the size of which is stated in the `BUFSIZE` macro), it exits. @@ -179,7 +179,7 @@ Along with copying its input samples to the output, it prints a warning message floatdump_f -It prints any floating point input samples. +It prints any floating point input samples. The format string used is `"%g "`. flowcontrol @@ -190,8 +190,8 @@ It copies `data_rate / reads_per_second` bytes from the input to the output, doi shift_math_cc It shifts the signal in the frequency domain by `rate`. -`rate` is a floating point number between -0.5 and 0.5. -`rate` is relative to the sampling rate. +`rate` is a floating point number between -0.5 and 0.5. +`rate` is relative to the sampling rate. Internally, a sine and cosine wave is generated to perform this function, and this function uses `math.h` for this purpose, which is quite accurate, but not always very fast. @@ -249,7 +249,7 @@ It uses fixed filters so it works only on predefined sample rates, for the actua amdemod_cf -It is an AM demodulator that uses `sqrt`. On some architectures `sqrt` can be directly calculated by dedicated CPU instructions, but on others it may be slower. +It is an AM demodulator that uses `sqrt`. On some architectures `sqrt` can be directly calculated by dedicated CPU instructions, but on others it may be slower. amdemod_estimator_cf @@ -278,7 +278,7 @@ Other parameters were explained above at `firdes_lowpass_f`. fir_decimate_cc [transition_bw [window]] -It is a decimator that keeps one sample out of `decimation_factor` samples. +It is a decimator that keeps one sample out of `decimation_factor` samples. To avoid aliasing, it runs a filter on the signal and removes spectral components above `0.5 × nyquist_frequency × decimation_factor`. `transition_bw` and `window` are the parameters of the filter. @@ -304,7 +304,7 @@ Parameters are described under `firdes_bandpass_c` and `firdes_lowpass_f`. agc_ff [hang_time [reference [attack_rate [decay_rate [max_gain [attack_wait [filter_alpha]]]]]]] -It is an automatic gain control function. +It is an automatic gain control function. - `hang_time` is the number of samples to wait before strating to increase the gain after a peak. - `reference` is the reference level for the AGC. It tries to keep the amplitude of the output signal close to that. @@ -322,7 +322,7 @@ It is a faster AGC that linearly changes the gain, taking the highest amplitude fft_cc [window [--octave] [--benchmark]] -It performs an FFT on the first `fft_size` samples out of `out_of_every_n_samples`, thus skipping `out_of_every_n_samples - fft_size` samples in the input. +It performs an FFT on the first `fft_size` samples out of `out_of_every_n_samples`, thus skipping `out_of_every_n_samples - fft_size` samples in the input. It can draw the spectrum by using `--octave`, for more information, look at the [Usage by example] section. @@ -357,7 +357,7 @@ It exchanges the first and second part of the FFT vector, to prepare it for the dsb_fc [q_value] -It converts a real signal to a double sideband complex signal centered around DC. +It converts a real signal to a double sideband complex signal centered around DC. It does so by generating a complex signal: * the real part of which is the input real signal, * the imaginary part of which is `q_value` (0 by default). @@ -387,6 +387,10 @@ It doubles every input sample. See the [buffer sizes](#buffer_sizes) section. + squelch_and_smeter_cc --fifo --outfifo + +This is a controllable squelch, which reads the squelch level input from `` and writes the power level output to ``. Both input and output are in the format of `%g\n`. While calculating the power level, it takes only every `` sample into consideration. It writes the S-meter value for every `` buffer to ``. If the squelch level is set to 0, it it forces the squelch to be open. If the squelch is closed, it fills the output with zero. + #### Control via pipes Some parameters can be changed while the `csdr` process is running. To achieve this, some `csdr` functions have special parameters. You have to supply a fifo previously created by the `mkfifo` command. Processing will only start after the first control command has been received by `csdr` over the FIFO. @@ -406,7 +410,7 @@ Processing will only start after the first control command has been received by By writing to the given FIFO file with the syntax below, you can control the shift rate: \n - + E.g. you can send `-0.05 0.02\n` #### Buffer sizes @@ -416,12 +420,12 @@ E.g. you can send `-0.05 0.02\n` * *dynamic buffer size determination:* input buffer size is recommended by the previous process, output buffer size is determined by the process, * *fixed buffer sizes*. -*csdr* can choose from two different buffer sizes by **default**. +*csdr* can choose from two different buffer sizes by **default**. * For operations handling the full-bandwidth I/Q data from the receiver, a buffer size of 16384 samples is used (see `env_csdr_fixed_big_bufsize` in the code). * For operations handling only a selected channel, a buffer size of 1024 samples is used (see `env_csdr_fixed_bufsize` in the code). *csdr* now has an experimental feature called **dynamic buffer size determination**, which is switched on by issuing `export CSDR_DYNAMIC_BUFSIZE_ON=1` in the shell before running `csdr`. If it is enabled: -* All `csdr` processes in a DSP chain acquire their recommended input buffer size from the previous `csdr` process. This information is in the first 8 bytes of the input stream. +* All `csdr` processes in a DSP chain acquire their recommended input buffer size from the previous `csdr` process. This information is in the first 8 bytes of the input stream. * Each process can decide whether to use this or choose another input buffer size (if that's more practical). * Every process sends out its output buffer size to the next process. Then it startss processing data. * The DSP chain should start with a `csdr setbuf ` process, which only copies data from the input to the output, but also sends out the given buffer size information to the next process. @@ -451,7 +455,7 @@ Example of initalization if the process generates N/D output samples for N input Example of initialization if the process allocates memory for itself, and it doesn't want to use the global buffers: - getbufsize(); //dummy + getbufsize(); //dummy sendbufsize(my_own_bufsize); Example of initialization if the process always works with a fixed output size, regardless of the input: diff --git a/csdr.c b/csdr.c index affcc4c..4ea544f 100644 --- a/csdr.c +++ b/csdr.c @@ -1,5 +1,5 @@ /* -This software is part of libcsdr, a set of simple DSP routines for +This software is part of libcsdr, a set of simple DSP routines for Software Defined Radio. Copyright (c) 2014, Andras Retzler @@ -35,14 +35,14 @@ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. #include #include #include -#include +#include #include #include #include #include #include #include -#include +#include #include "libcsdr.h" #include "libcsdr_gpl.h" #include "ima_adpcm.h" @@ -71,7 +71,9 @@ char usage[]= " floatdump_f\n" " flowcontrol [prebuffer_sec] [thrust]\n" " shift_math_cc \n" +" shift_math_cc --fifo \n" " shift_addition_cc \n" +" shift_addition_cc --fifo \n" " shift_addition_cc_test\n" " shift_table_cc [table_size]\n" " decimating_shift_addition_cc [decimation]\n" @@ -95,6 +97,7 @@ char usage[]= " logpower_cf [add_db]\n" " fft_benchmark [--benchmark]\n" " bandpass_fir_fft_cc [window]\n" +" bandpass_fir_fft_cc --fifo [window]\n" " encode_ima_adpcm_s16_u8\n" " decode_ima_adpcm_u8_s16\n" " compress_fft_adpcm_f_u8 \n" @@ -107,6 +110,7 @@ char usage[]= " monos2stereo_s16\n" " setbuf \n" " fft_exchange_sides_ff \n" +" squelch_and_smeter_cc --fifo --outfifo \n" " \n" ; @@ -192,12 +196,12 @@ int read_fifo_ctl(int fd, char* format, ...) static int buffer_index=0; int bytes_read=read(fd,buffer+buffer_index,(RFCTL_BUFSIZE-buffer_index)*sizeof(char)); if(bytes_read<=0) return 0; - + int prev_newline_at=0; int last_newline_at=0; - for(int i=0;i=4) sscanf(argv[3],"%d",&buf_times); if(!sendbufsize(initialize_buffers())) return -2; for(int i=0;i/dev/null //csdr yes_f 1 1000000 | time csdr shift_addition_cc 0.2 >/dev/null //csdr yes_f 1 1000000 | time csdr shift_table_cc 0.2 >/dev/null @@ -497,7 +501,7 @@ int main(int argc, char *argv[]) if(!strcmp(argv[1],"shift_table_cc")) { bigbufs=1; - if(argc<=2) return badsyntax("need required parameter (rate)"); + if(argc<=2) return badsyntax("need required parameter (rate)"); float starting_phase=0; float rate; int table_size=65536; @@ -521,7 +525,7 @@ int main(int argc, char *argv[]) if(!strcmp(argv[1],"decimating_shift_addition_cc")) { bigbufs=1; - if(argc<=2) return badsyntax("need required parameter (rate)"); + if(argc<=2) return badsyntax("need required parameter (rate)"); float starting_phase=0; float rate; int decimation=1; @@ -558,7 +562,7 @@ int main(int argc, char *argv[]) } else { - if(argc<=2) return badsyntax("need required parameter (rate)"); + if(argc<=2) return badsyntax("need required parameter (rate)"); sscanf(argv[2],"%g",&rate); } @@ -595,7 +599,7 @@ int main(int argc, char *argv[]) if(!strcmp(argv[1],"shift_addition_cc_test")) { - if(argc<=2) return badsyntax("need required parameter (rate)"); + if(argc<=2) return badsyntax("need required parameter (rate)"); float rate; sscanf(argv[2],"%g",&rate); //if(initialize_buffers()) return -2; //most likely we don't need this here @@ -607,7 +611,7 @@ int main(int argc, char *argv[]) if(!strcmp(argv[1],"dcblock_ff")) { static dcblock_preserve_t dcp; //will be 0 as .bss is set to 0 - if(!sendbufsize(initialize_buffers())) return -2; + if(!sendbufsize(initialize_buffers())) return -2; for(;;) { FEOF_CHECK; @@ -682,7 +686,7 @@ int main(int argc, char *argv[]) } if(!strcmp(argv[1],"deemphasis_wfm_ff")) { - if(argc<=3) return badsyntax("need required parameters (sample rate, tau)"); + if(argc<=3) return badsyntax("need required parameters (sample rate, tau)"); if(!sendbufsize(initialize_buffers())) return -2; int sample_rate; sscanf(argv[2],"%d",&sample_rate); @@ -707,19 +711,19 @@ int main(int argc, char *argv[]) { FEOF_CHECK; FREAD_R; - int nan_detect=0; - for(int i=0; i=3) sscanf(argv[2],"%d",&input.input_size); - getbufsize(); //dummy + getbufsize(); //dummy sendbufsize(input.input_size); input.reference=1.0; @@ -996,7 +1000,7 @@ int main(int argc, char *argv[]) //input.max_peak_ratio=12.0; //if(argc>=5) sscanf(argv[3],"%g",&input.max_peak_ratio); - + input.buffer_1=(float*)calloc(input.input_size,sizeof(float)); input.buffer_2=(float*)calloc(input.input_size,sizeof(float)); input.buffer_input=(float*)malloc(sizeof(float)*input.input_size); @@ -1005,7 +1009,7 @@ int main(int argc, char *argv[]) { FEOF_CHECK; fread(input.buffer_input, sizeof(float), input.input_size, stdin); - fastagc_ff(&input, agc_output_buffer); + fastagc_ff(&input, agc_output_buffer); fwrite(agc_output_buffer, sizeof(float), input.input_size, stdout); TRY_YIELD; } @@ -1014,11 +1018,11 @@ int main(int argc, char *argv[]) int suboptimal; if( (suboptimal=!strcmp(argv[1],"suboptimal_rational_resampler_ff"))||(!strcmp(argv[1],"rational_resampler_ff")) ) { - + //last@2014-11-06: ./docompile; ./csdr yes_f 1.0 | ./csdr suboptimal_rational_resampler_ff 5 2 //Process the params - if(argc<=3) return badsyntax("need required parameters (interpolation, decimation)"); + if(argc<=3) return badsyntax("need required parameters (interpolation, decimation)"); int interpolation; sscanf(argv[2],"%d",&interpolation); int decimation; @@ -1050,7 +1054,7 @@ int main(int argc, char *argv[]) int taps_length = firdes_filter_len(transition_bw); float* taps = (float*)malloc(sizeof(float)*taps_length); rational_resampler_get_lowpass_f(taps, taps_length, interpolation, decimation, window); - + static rational_resampler_ff_t d; //in .bss => initialized to zero for(;;) @@ -1071,7 +1075,7 @@ int main(int argc, char *argv[]) if(!strcmp(argv[1],"fractional_decimator_ff")) { //Process the params - if(argc<=2) return badsyntax("need required parameters (rate)"); + if(argc<=2) return badsyntax("need required parameters (rate)"); float rate; sscanf(argv[2],"%g",&rate); @@ -1086,7 +1090,7 @@ int main(int argc, char *argv[]) else fprintf(stderr,"fractional_decimator_ff: window = %s\n",firdes_get_string_from_window(window)); if(!initialize_buffers()) return -2; - sendbufsize(the_bufsize / rate); + sendbufsize(the_bufsize / rate); if(rate==1) clone_(the_bufsize); //copy input to output in this special case (and stick in this function). @@ -1094,7 +1098,7 @@ int main(int argc, char *argv[]) int taps_length = firdes_filter_len(transition_bw); fprintf(stderr,"fractional_decimator_ff: taps_length = %d\n",taps_length); float* taps = (float*)malloc(sizeof(float)*taps_length); - firdes_lowpass_f(taps, taps_length, 0.59*0.5/(rate-transition_bw), window); //0.6 const to compensate rolloff + firdes_lowpass_f(taps, taps_length, 0.59*0.5/(rate-transition_bw), window); //0.6 const to compensate rolloff //for(int=0;i initialized to zero @@ -1112,10 +1116,10 @@ int main(int argc, char *argv[]) if(!strcmp(argv[1],"fft_cc")) { - if(argc<=3) return badsyntax("need required parameters (fft_size, out_of_every_n_samples)"); + if(argc<=3) return badsyntax("need required parameters (fft_size, out_of_every_n_samples)"); int fft_size; sscanf(argv[2],"%d",&fft_size); - if(log2n(fft_size)==-1) return badsyntax("fft_size should be power of 2"); + if(log2n(fft_size)==-1) return badsyntax("fft_size should be power of 2"); int every_n_samples; sscanf(argv[3],"%d",&every_n_samples); int benchmark=0; @@ -1125,12 +1129,12 @@ int main(int argc, char *argv[]) { window=firdes_get_window_from_string(argv[4]); } - if(argc>=6) + if(argc>=6) { benchmark|=!strcmp("--benchmark",argv[5]); octave|=!strcmp("--octave",argv[5]); } - if(argc>=7) + if(argc>=7) { benchmark|=!strcmp("--benchmark",argv[6]); octave|=!strcmp("--octave",argv[6]); @@ -1172,7 +1176,7 @@ int main(int argc, char *argv[]) printf("fftdata=["); //we have to swap the two parts of the array to get a valid spectrum for(int i=fft_size/2;i=3) sscanf(argv[2],"%g",&add_db); - if(!sendbufsize(initialize_buffers())) return -2; + if(!sendbufsize(initialize_buffers())) return -2; for(;;) { @@ -1203,7 +1207,7 @@ int main(int argc, char *argv[]) if(!strcmp(argv[1],"fft_exchange_sides_ff")) { - if(argc<=2) return badsyntax("need required parameters (fft_size)"); + if(argc<=2) return badsyntax("need required parameters (fft_size)"); int fft_size; sscanf(argv[2],"%d",&fft_size); if(!getbufsize()) return -2; //dummy @@ -1226,13 +1230,13 @@ int main(int argc, char *argv[]) #define COMPRESS_FFT_PAD_N 10 //We will pad the FFT at the beginning, with the first value of the input data, COMPRESS_FFT_PAD_N times. -//No, this is not advanced DSP, just the ADPCM codec produces some gabarge samples at the beginning, -//so we just add data to become garbage and get skipped. +//No, this is not advanced DSP, just the ADPCM codec produces some gabarge samples at the beginning, +//so we just add data to become garbage and get skipped. //COMPRESS_FFT_PAD_N should be even. if(!strcmp(argv[1],"compress_fft_adpcm_f_u8")) { - if(argc<=2) return badsyntax("need required parameters (fft_size)"); + if(argc<=2) return badsyntax("need required parameters (fft_size)"); int fft_size; sscanf(argv[2],"%d",&fft_size); int real_data_size=fft_size+COMPRESS_FFT_PAD_N; @@ -1260,12 +1264,12 @@ int main(int argc, char *argv[]) if(!strcmp(argv[1],"fft_benchmark")) { - if(argc<=3) return badsyntax("need required parameters (fft_size, fft_cycles)"); + if(argc<=3) return badsyntax("need required parameters (fft_size, fft_cycles)"); int fft_size; sscanf(argv[2],"%d",&fft_size); int fft_cycles; sscanf(argv[3],"%d",&fft_cycles); - + int benchmark=(argc>=5)&&!strcmp(argv[4],"--benchmark"); fprintf(stderr,"fft_benchmark: FFT library used: %s\n",FFT_LIBRARY_USED); @@ -1274,20 +1278,20 @@ int main(int argc, char *argv[]) //fill input with random data srand(time(NULL)); - for(int i=0;ioutput,i)=qof(plan_inverse_2->output,i)=0; - + FFT_PLAN_T* plan_inverse_2 = make_fft_c2c(fft_size, output_fourier, output_2, 0, 1); + //we initialize this buffer to 0 as it will be taken as the overlap source for the first time: + for(int i=0;ioutput,i)=qof(plan_inverse_2->output,i)=0; + for(int i=input_size;i0) flowcontrol_bufindex+=read_return; + if(read_return>0) flowcontrol_bufindex+=read_return; if(flowcontrol_is_buffering) @@ -1537,7 +1541,7 @@ int main(int argc, char *argv[]) if(!strcmp(argv[1],"through")) { - struct timespec start_time, end_time; + struct timespec start_time, end_time; if(!sendbufsize(initialize_buffers())) return -2; int time_now_sec=0; @@ -1545,21 +1549,21 @@ int main(int argc, char *argv[]) unsigned char* through_buffer; through_buffer = (unsigned char*)malloc(the_bufsize*sizeof(float)); - - + + for(;;) { FEOF_CHECK; fread(through_buffer, sizeof(float), the_bufsize, stdin); - if(!time_now_sec) - { - time_now_sec=1; - clock_gettime(CLOCK_MONOTONIC_RAW, &start_time); - } - else + if(!time_now_sec) { - clock_gettime(CLOCK_MONOTONIC_RAW, &end_time); + time_now_sec=1; + clock_gettime(CLOCK_MONOTONIC_RAW, &start_time); + } + else + { + clock_gettime(CLOCK_MONOTONIC_RAW, &end_time); float timetaken; if(time_now_sec<(timetaken=TIME_TAKEN(start_time,end_time))) { @@ -1583,20 +1587,20 @@ int main(int argc, char *argv[]) { FEOF_CHECK; FREAD_R; - for(int i=0;i)"); + fprintf(stderr, "squelch_and_power_cc: initial squelch level is %g\n", squelch_level); + if((argc<=5)||((argc>5)&&(strcmp(argv[4],"--outfifo")))) return badsyntax("need required parameter (--outfifo )"); + int fd2 = open(argv[5], O_WRONLY); + if(fd2==-1) return badsyntax("error while opening --outfifo"); + int flags = fcntl(fd2, F_GETFL, 0); + fcntl(fd2, F_SETFL, flags | O_NONBLOCK); + if(argc<=6) return badsyntax("need required parameter (use_every_nth)"); + sscanf(argv[6],"%d",&decimation); + if(decimation<=0) return badsyntax("use_every_nth <= 0 is invalid"); + if(argc<=7) return badsyntax("need required parameter (report_every_nth)"); + sscanf(argv[7],"%d",&report_every_nth); + if(report_every_nth<=0) return badsyntax("report_every_nth <= 0 is invalid"); + for(;;) + { + FEOF_CHECK; + FREAD_C; //read input data + power = get_power_c((complexf*)input_buffer, the_bufsize, decimation); + if(report_cntr++>report_every_nth) + { + report_cntr=0; + power_value_buf_size=snprintf(power_value_buf,100,"%g\n",power); + write(fd2,power_value_buf,power_value_buf_size*sizeof(char)); + } + if(squelch_level==0||power>=squelch_level) + { + //fprintf(stderr,"P"); + fwrite(input_buffer, sizeof(complexf), the_bufsize, stdout); + } + else + { + //fprintf(stderr,"S"); + fwrite(zerobuf, sizeof(complexf), the_bufsize, stdout); + } + if(read_fifo_ctl(fd,"%g\n",&squelch_level)) fprintf(stderr, "squelch_and_power_cc: new squelch level is %g\n", squelch_level); + TRY_YIELD; + } } if(!strcmp(argv[1],"none")) { return 0; } - + return badsyntax("function name given in argument 1 does not exist. Possible causes:\n- You mistyped the commandline.\n- You need to update csdr to a newer version (if available)."); } - - diff --git a/libcsdr.c b/libcsdr.c index 81474db..1b4125b 100644 --- a/libcsdr.c +++ b/libcsdr.c @@ -1,5 +1,5 @@ /* -This software is part of libcsdr, a set of simple DSP routines for +This software is part of libcsdr, a set of simple DSP routines for Software Defined Radio. Copyright (c) 2014, Andras Retzler @@ -40,14 +40,14 @@ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. #include /* - _ _ __ _ _ - (_) | | / _| | | (_) - __ ___ _ __ __| | _____ __ | |_ _ _ _ __ ___| |_ _ ___ _ __ ___ + _ _ __ _ _ + (_) | | / _| | | (_) + __ ___ _ __ __| | _____ __ | |_ _ _ _ __ ___| |_ _ ___ _ __ ___ \ \ /\ / / | '_ \ / _` |/ _ \ \ /\ / / | _| | | | '_ \ / __| __| |/ _ \| '_ \/ __| \ V V /| | | | | (_| | (_) \ V V / | | | |_| | | | | (__| |_| | (_) | | | \__ \ \_/\_/ |_|_| |_|\__,_|\___/ \_/\_/ |_| \__,_|_| |_|\___|\__|_|\___/|_| |_|___/ - - + + */ #define MFIRDES_GWS(NAME) \ @@ -95,7 +95,7 @@ float firdes_wkernel_boxcar(float rate) } float (*firdes_get_window_kernel(window_t window))(float) -{ +{ if(window==WINDOW_HAMMING) return firdes_wkernel_hamming; else if(window==WINDOW_BLACKMAN) return firdes_wkernel_blackman; else if(window==WINDOW_BOXCAR) return firdes_wkernel_boxcar; @@ -103,14 +103,14 @@ float (*firdes_get_window_kernel(window_t window))(float) } /* - ______ _____ _____ __ _ _ _ _ _ - | ____|_ _| __ \ / _(_) | | | | (_) - | |__ | | | |__) | | |_ _| | |_ ___ _ __ __| | ___ ___ _ __ _ _ __ - | __| | | | _ / | _| | | __/ _ \ '__| / _` |/ _ \/ __| |/ _` | '_ \ + ______ _____ _____ __ _ _ _ _ _ + | ____|_ _| __ \ / _(_) | | | | (_) + | |__ | | | |__) | | |_ _| | |_ ___ _ __ __| | ___ ___ _ __ _ _ __ + | __| | | | _ / | _| | | __/ _ \ '__| / _` |/ _ \/ __| |/ _` | '_ \ | | _| |_| | \ \ | | | | | || __/ | | (_| | __/\__ \ | (_| | | | | |_| |_____|_| \_\ |_| |_|_|\__\___|_| \__,_|\___||___/_|\__, |_| |_| - __/ | - |___/ + __/ | + |___/ */ void firdes_lowpass_f(float *output, int length, float cutoff_rate, window_t window) @@ -127,7 +127,7 @@ void firdes_lowpass_f(float *output, int length, float cutoff_rate, window_t win output[middle-i]=output[middle+i]=(sin(2*PI*cutoff_rate*i)/i)*window_function((float)i/middle); //printf("%g %d %d %d %d | %g\n",output[middle-i],i,middle,middle+i,middle-i,sin(2*PI*cutoff_rate*i)); } - + //Normalize filter kernel float sum=0; for(int i=0;i2*PI) phase-=2*PI; //@@firdes_bandpass_c - while(phase<0) phase+=2*PI; + while(phase<0) phase+=2*PI; iof(output,i)=cosval*realtaps[i]; qof(output,i)=sinval*realtaps[i]; //output[i] := realtaps[i] * e^j*w @@ -173,14 +173,14 @@ int firdes_filter_len(float transition_bw) } /* - _____ _____ _____ __ _ _ - | __ \ / ____| __ \ / _| | | (_) - | | | | (___ | |__) | | |_ _ _ _ __ ___| |_ _ ___ _ __ ___ + _____ _____ _____ __ _ _ + | __ \ / ____| __ \ / _| | | (_) + | | | | (___ | |__) | | |_ _ _ _ __ ___| |_ _ ___ _ __ ___ | | | |\___ \| ___/ | _| | | | '_ \ / __| __| |/ _ \| '_ \/ __| | |__| |____) | | | | | |_| | | | | (__| |_| | (_) | | | \__ \ |_____/|_____/|_| |_| \__,_|_| |_|\___|\__|_|\___/|_| |_|___/ -*/ +*/ float shift_math_cc(complexf *input, complexf* output, int input_size, float rate, float starting_phase) { @@ -239,7 +239,7 @@ float shift_table_cc(complexf* input, complexf* output, int input_size, float ra //float vphase=fmodf(phase,PI/2); //between 0 and 90deg int quadrant=phase/(PI/2); //between 0 and 3 float vphase=phase-quadrant*(PI/2); - sin_index=(vphase/(PI/2))*table_data.table_size; + sin_index=(vphase/(PI/2))*table_data.table_size; cos_index=table_data.table_size-1-sin_index; if(quadrant&1) //in quadrant 1 and 3 { @@ -280,8 +280,8 @@ int fir_decimate_cc(complexf *input, complexf *output, int input_size, int decim for(int i=0; iinput_size) break; - register float acci=0; - register float accq=0; + register float acci=0; + register float accq=0; register int ti=0; register float* pinput=(float*)&(input[i+ti]); @@ -289,10 +289,10 @@ int fir_decimate_cc(complexf *input, complexf *output, int input_size, int decim register float* ptaps_end=taps+taps_length; float quad_acciq [8]; - + /* q0, q1: input signal I sample and Q sample -q2: taps +q2: taps q4, q5: accumulator for I branch and Q branch (will be the output) */ @@ -300,7 +300,7 @@ q4, q5: accumulator for I branch and Q branch (will be the output) " vmov.f32 q4, #0.0\n\t" //another way to null the accumulators " vmov.f32 q5, #0.0\n\t" "for_fdccasm: vld2.32 {q0-q1}, [%[pinput]]!\n\t" //load q0 and q1 directly from the memory address stored in pinput, with interleaving (so that we get the I samples in q0 and the Q samples in q1), also increment the memory address in pinput (hence the "!" mark) //http://community.arm.com/groups/processors/blog/2010/03/17/coding-for-neon--part-1-load-and-stores - " vld1.32 {q2}, [%[ptaps]]!\n\t" + " vld1.32 {q2}, [%[ptaps]]!\n\t" " vmla.f32 q4, q0, q2\n\t" //quad_acc_i += quad_input_i * quad_taps_1 //http://stackoverflow.com/questions/3240440/how-to-use-the-multiply-and-accumulate-intrinsics-in-arm-cortex-a8 //http://infocenter.arm.com/help/index.jsp?topic=/com.arm.doc.dui0489e/CIHEJBIE.html " vmla.f32 q5, q1, q2\n\t" //quad_acc_q += quad_input_q * quad_taps_1 " cmp %[ptaps], %[ptaps_end]\n\t" //if(ptaps == ptaps_end) @@ -311,13 +311,13 @@ q4, q5: accumulator for I branch and Q branch (will be the output) [pinput]"+r"(pinput), [ptaps]"+r"(ptaps) //output operand list : [ptaps_end]"r"(ptaps_end), [quad_acci]"r"(quad_acciq), [quad_accq]"r"(quad_acciq+4) //input operand list - : + : "memory", "q0", "q1", "q2", "q4", "q5", "cc" //clobber list ); //original for loops for reference: - //for(int ti=0; ti> [%d] %g \n", n, quad_acciq[n]); iof(output,oi)=quad_acciq[0]+quad_acciq[1]+quad_acciq[2]+quad_acciq[3]; //we're still not ready, as we have to add up the contents of a quad accumulator register to get a single accumulated value qof(output,oi)=quad_acciq[4]+quad_acciq[5]+quad_acciq[6]+quad_acciq[7]; @@ -381,7 +381,7 @@ int fir_decimate_cc(complexf *input, complexf *output, int input_size, int decim rational_resampler_ff_t rational_resampler_ff(float *input, float *output, int input_size, int interpolation, int decimation, float *taps, int taps_length, int last_taps_delay) { - + //Theory: http://www.dspguru.com/dsp/faqs/multirate/resampling //oi: output index, i: tap index int output_size=input_size*interpolation/decimation; @@ -412,7 +412,7 @@ rational_resampler_ff_t rational_resampler_ff(float *input, float *output, int i /* -The greatest challenge in resampling is figuring out which tap should be applied to which sample. +The greatest challenge in resampling is figuring out which tap should be applied to which sample. Typical test patterns for rational_resampler_ff: @@ -457,13 +457,13 @@ float inline fir_one_pass_ff(float* input, float* taps, int taps_length) fractional_decimator_ff_t fractional_decimator_ff(float* input, float* output, int input_size, float rate, float *taps, int taps_length, fractional_decimator_ff_t d) { if(rate<=1.0) return d; //sanity check, can't decimate <=1.0 - //This routine can handle floating point decimation rates. - //It linearly interpolates between two samples that are taken into consideration from the filtered input. + //This routine can handle floating point decimation rates. + //It linearly interpolates between two samples that are taken into consideration from the filtered input. int oi=0; int index_high; float where=d.remain; float result_high, result_low; - if(where==0.0) //in the first iteration index_high may be zero (so using the item index_high-1 would lead to invalid memory access). + if(where==0.0) //in the first iteration index_high may be zero (so using the item index_high-1 would lead to invalid memory access). { output[oi++]=fir_one_pass_ff(input,taps,taps_length); where+=rate; @@ -474,7 +474,7 @@ fractional_decimator_ff_t fractional_decimator_ff(float* input, float* output, i for(;(index_high=ceilf(where))+taps_lengthoutput; complexf* out = plan_inverse->input; - + for(int i=0;isize;i++) //@apply_fir_fft_cc: multiplication { iof(out,i)=iof(in,i)*iof(taps_fft,i)-qof(in,i)*qof(taps_fft,i); qof(out,i)=iof(in,i)*qof(taps_fft,i)+qof(in,i)*iof(taps_fft,i); } - + //calculate inverse FFT on multiplied buffer fft_execute(plan_inverse); - + //add the overlap of the previous segment complexf* result = plan_inverse->output; @@ -516,35 +516,35 @@ void apply_fir_fft_cc(FFT_PLAN_T* plan, FFT_PLAN_T* plan_inverse, complexf* taps iof(result,i)/=plan->size; qof(result,i)/=plan->size; } - + for(int i=0;ipeak_2) target_peak=input->peak_2; if(target_peakpeak_1) target_peak=input->peak_1; - + //we change the gain linearly on the apply_block from the last_gain to target_gain. float target_gain=input->reference/target_peak; if(target_gain>FASTAGC_MAX_GAIN) target_gain=FASTAGC_MAX_GAIN; @@ -668,9 +668,9 @@ void fastagc_ff(fastagc_ff_t* input, float* output) } /* - ______ __ __ _ _ _ _ - | ____| \/ | | | | | | | | | - | |__ | \ / | __| | ___ _ __ ___ ___ __| |_ _| | __ _| |_ ___ _ __ ___ + ______ __ __ _ _ _ _ + | ____| \/ | | | | | | | | | + | |__ | \ / | __| | ___ _ __ ___ ___ __| |_ _| | __ _| |_ ___ _ __ ___ | __| | |\/| | / _` |/ _ \ '_ ` _ \ / _ \ / _` | | | | |/ _` | __/ _ \| '__/ __| | | | | | | | (_| | __/ | | | | | (_) | (_| | |_| | | (_| | || (_) | | \__ \ |_| |_| |_| \__,_|\___|_| |_| |_|\___/ \__,_|\__,_|_|\__,_|\__\___/|_| |___/ @@ -721,33 +721,33 @@ complexf fmdemod_quadri_cf(complexf* input, float* output, int input_size, float temp_dq[0]=qof(input,0)-last_sample.q; for (int i=1; i tau = 75e-6 WFM transmission in EU: 50 us -> tau = 50e-6 @@ -767,7 +767,7 @@ float deemphasis_wfm_ff (float* input, float* output, int input_size, float tau, float dt = 1.0/sample_rate; float alpha = dt/(tau+dt); if(is_nan(last_output)) last_output=0.0; //if last_output is NaN - output[0]=alpha*input[0]+(1-alpha)*last_output; + output[0]=alpha*input[0]+(1-alpha)*last_output; for (int i=1;ioutput[i])?-max_amplitude:output[i]; - } + } } void gain_ff(float* input, float* output, int input_size, float gain) @@ -818,20 +818,40 @@ void gain_ff(float* input, float* output, int input_size, float gain) for(int i=0;i 0) ? new_amplitude / amplitude_now : 0; iof(output,i)=iof(input,i)*gain; qof(output,i)=qof(input,i)*gain; - } + } } /* - ______ _ ______ _ _______ __ - | ____| | | | ____| (_) |__ __| / _| - | |__ __ _ ___| |_ | |__ ___ _ _ _ __ _ ___ _ __ | |_ __ __ _ _ __ ___| |_ ___ _ __ _ __ ___ - | __/ _` / __| __| | __/ _ \| | | | '__| |/ _ \ '__| | | '__/ _` | '_ \/ __| _/ _ \| '__| '_ ` _ \ + ______ _ ______ _ _______ __ + | ____| | | | ____| (_) |__ __| / _| + | |__ __ _ ___| |_ | |__ ___ _ _ _ __ _ ___ _ __ | |_ __ __ _ _ __ ___| |_ ___ _ __ _ __ ___ + | __/ _` / __| __| | __/ _ \| | | | '__| |/ _ \ '__| | | '__/ _` | '_ \/ __| _/ _ \| '__| '_ ` _ \ | | | (_| \__ \ |_ | | | (_) | |_| | | | | __/ | | | | | (_| | | | \__ \ || (_) | | | | | | | | - |_| \__,_|___/\__| |_| \___/ \__,_|_| |_|\___|_| |_|_| \__,_|_| |_|___/_| \___/|_| |_| |_| |_| + |_| \__,_|___/\__| |_| \___/ \__,_|_| |_|\___|_| |_|_| \__,_|_| |_|___/_| \___/|_| |_| |_| |_| */ int log2n(int x) { int result=-1; - for(int i=0;i<31;i++) + for(int i=0;i<31;i++) { if((x>>i)&1) //@@log2n { @@ -892,7 +912,7 @@ int next_pow2(int x) { int pow2; //portability? (31 is the problem) - for(int i=0;i<31;i++) + for(int i=0;i<31;i++) { if(x<(pow2=1< convert_u8_f } @@ -968,7 +988,7 @@ void convert_f_s8(float* input, signed char* output, int input_size) void convert_f_s16(float* input, short* output, int input_size) { - /*for(int i=0;i1.0) input[i]=1.0; if(input[i]<-1.0) input[i]=-1.0; @@ -989,5 +1009,3 @@ int trivial_vectorize() } return c[0]; } - - diff --git a/libcsdr.h b/libcsdr.h index 79696b3..8bcbc2a 100644 --- a/libcsdr.h +++ b/libcsdr.h @@ -1,5 +1,5 @@ /* -This software is part of libcsdr, a set of simple DSP routines for +This software is part of libcsdr, a set of simple DSP routines for Software Defined Radio. Copyright (c) 2014, Andras Retzler @@ -32,14 +32,14 @@ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. #define MIN_M(x,y) (((x)>(y))?(y):(x)) /* - _____ _ - / ____| | | - | | ___ _ __ ___ _ __ | | _____ __ - | | / _ \| '_ ` _ \| '_ \| |/ _ \ \/ / - | |___| (_) | | | | | | |_) | | __/> < - \_____\___/|_| |_| |_| .__/|_|\___/_/\_\ - | | - |_| + _____ _ + / ____| | | + | | ___ _ __ ___ _ __ | | _____ __ + | | / _ \| '_ ` _ \| '_ \| |/ _ \ \/ / + | |___| (_) | | | | | | |_) | | __/> < + \_____\___/|_| |_| |_| .__/|_|\___/_/\_\ + | | + |_| */ typedef struct complexf_s { float i; float q; } complexf; @@ -64,7 +64,7 @@ typedef struct complexf_s { float i; float q; } complexf; #define PI ((float)3.14159265358979323846) //window -typedef enum window_s +typedef enum window_s { WINDOW_BOXCAR, WINDOW_BLACKMAN, WINDOW_HAMMING } window_t; @@ -113,7 +113,7 @@ float fastdcblock_ff(float* input, float* output, int input_size, float last_dc_ typedef struct fastagc_ff_s { - float* buffer_1; + float* buffer_1; float* buffer_2; float* buffer_input; //it is the actual input buffer to fill float peak_1; @@ -150,7 +150,7 @@ fractional_decimator_ff_t fractional_decimator_ff(float* input, float* output, i typedef struct shift_table_data_s { float* table; - int table_size; + int table_size; } shift_table_data_t; void shift_table_deinit(shift_table_data_t table_data); shift_table_data_t shift_table_init(int table_size); @@ -161,6 +161,8 @@ int log2n(int x); int next_pow2(int x); void apply_fir_fft_cc(FFT_PLAN_T* plan, FFT_PLAN_T* plan_inverse, complexf* taps_fft, complexf* last_overlap, int overlap_size); void gain_ff(float* input, float* output, int input_size, float gain); +float get_power_f(float* input, int input_size, int decimation); +float get_power_c(complexf* input, int input_size, int decimation); void add_dcoffset_cc(complexf* input, complexf* output, int input_size); float fmmod_fc(float* input, complexf* output, int input_size, float last_phase);