diff --git a/Makefile b/Makefile index 06f0709..920bb46 100644 --- a/Makefile +++ b/Makefile @@ -29,7 +29,7 @@ LIBSOURCES = fft_fftw.c libcsdr_wrapper.c #SOURCES = csdr.c $(LIBSOURCES) cpufeature = $(if $(findstring $(1),$(shell cat /proc/cpuinfo)),$(2)) -PARAMS_SSE = $(call cpufeature,sse,-msse) $(call cpufeature,sse2,-msse2) $(call cpufeature,sse3,-msse3) $(call cpufeature,sse4,-msse4) $(call cpufeature,sse4_1,-msse4.1) $(call cpufeature,sse4_2,-msse4.2) -mfpmath=sse +PARAMS_SSE = $(call cpufeature,sse,-msse) $(call cpufeature,sse2,-msse2) $(call cpufeature,sse3,-msse3) $(call cpufeature,sse4a,-msse4a) $(call cpufeature,sse4_1,-msse4.1) $(call cpufeature,sse4_2,-msse4.2 -msse4) -mfpmath=sse PARAMS_NEON = -mfloat-abi=hard -march=armv7-a -mtune=cortex-a8 -mfpu=neon -mvectorize-with-neon-quad -funsafe-math-optimizations -Wformat=0 -DNEON_OPTS #tnx Jan Szumiec for the Raspberry Pi support PARAMS_RASPI = -mfloat-abi=hard -mcpu=arm1176jzf-s -mfpu=vfp -funsafe-math-optimizations -Wformat=0 @@ -92,7 +92,7 @@ emcc-get-deps: emmake make; \ emmake make install emcc: - emcc -O3 -Isdr.js/$(FFTW_PACKAGE)/api -Lsdr.js/$(FFTW_PACKAGE)/emscripten-lib -o sdr.js/sdrjs-compiled.js fft_fftw.c libcsdr_wrapper.c -DLIBCSDR_GPL -DUSE_IMA_ADPCM -DUSE_FFTW -lfftw3f -s EXPORTED_FUNCTIONS="`python sdr.js/exported_functions.py`" + emcc -O3 -Isdr.js/$(FFTW_PACKAGE)/api -Lsdr.js/$(FFTW_PACKAGE)/emscripten-lib -o sdr.js/sdrjs-compiled.js fft_fftw.c libcsdr_wrapper.c -s TOTAL_MEMORY=67108864 -DLIBCSDR_GPL -DUSE_IMA_ADPCM -DUSE_FFTW -lfftw3f -s EXPORTED_FUNCTIONS="`python sdr.js/exported_functions.py`" cat sdr.js/sdrjs-header.js sdr.js/sdrjs-compiled.js sdr.js/sdrjs-footer.js > sdr.js/sdr.js emcc-beautify: bash -c 'type js-beautify >/dev/null 2>&1; if [ $$? -eq 0 ]; then js-beautify sdr.js/sdr.js >sdr.js/sdr.js.beautiful; mv sdr.js/sdr.js.beautiful sdr.js/sdr.js; fi' diff --git a/README.md b/README.md index 0aedfcd..db2816d 100644 --- a/README.md +++ b/README.md @@ -12,14 +12,22 @@ Most of the code is available under the permissive BSD license, with some option How to compile -------------- -The project was only tested on Linux. It has the following dependencies: `libfftw3-dev` make sudo make install +The project was only tested on Linux. It has the following dependencies: `libfftw3-dev` + If you compile on ARM, please edit the Makefile and tailor `PARAMS_NEON` for your CPU. + To run the examples, you will also need rtl_sdr from Osmocom, and the following packages (at least on Debian): `mplayer octave gnuplot gnuplot-x11` +If you compile *fftw3* from sources for use with *libcsdr*, you need to configure it with 32-bit float support enabled: + + ./configure --enable-float + +(This is for *fftw3*, not *libcsdr*. You do not need to run the configure script before compiling *libcsdr*.) + Credits ------- The library was written by Andras Retzler, HA7ILM <>. @@ -31,19 +39,19 @@ Usage by example ### Demodulate WFM - rtl_sdr -s 240000 -f 89500000 -g 20 - | csdr convert_u8_f | csdr fmdemod_quadri_cf | csdr fractional_decimator_ff 5 | csdr deemphasis_wfm_ff 48000 50e-6 | csdr convert_f_i16 | mplayer -cache 1024 -quiet -rawaudio samplesize=2:channels=1:rate=48000 -demuxer rawaudio - + rtl_sdr -s 240000 -f 89500000 -g 20 - | csdr convert_u8_f | csdr fmdemod_quadri_cf | csdr fractional_decimator_ff 5 | csdr deemphasis_wfm_ff 48000 50e-6 | csdr convert_f_s16 | mplayer -cache 1024 -quiet -rawaudio samplesize=2:channels=1:rate=48000 -demuxer rawaudio - - 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. ### Demodulate WFM: advanced - rtl_sdr -s 2400000 -f 89300000 -g 20 - | csdr convert_u8_f | csdr shift_addition_cc -0.085 | csdr fir_decimate_cc 10 0.05 HAMMING | csdr fmdemod_quadri_cf | csdr fractional_decimator_ff 5 | csdr deemphasis_wfm_ff 48000 50e-6 | csdr convert_f_i16 | mplayer -cache 1024 -quiet -rawaudio samplesize=2:channels=1:rate=48000 -demuxer rawaudio - + rtl_sdr -s 2400000 -f 89300000 -g 20 - | csdr convert_u8_f | csdr shift_addition_cc -0.085 | csdr fir_decimate_cc 10 0.05 HAMMING | csdr fmdemod_quadri_cf | csdr fractional_decimator_ff 5 | csdr deemphasis_wfm_ff 48000 50e-6 | csdr convert_f_s16 | mplayer -cache 1024 -quiet -rawaudio samplesize=2:channels=1:rate=48000 -demuxer rawaudio - - We want to listen to one radio station, but input signal contains multiple stations, and its bandwidth is too large for sending it directly to the FM demodulator. - We shift the signal to the center frequency of the station we want to receive: `-0.085*2400000 = -204000`, so basically we will listen to the radio station centered at 89504000 Hz. @@ -58,13 +66,13 @@ 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 - + + 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_s16 | 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). - Also there is a python hack to calculate the relative shift offset. The real receiver frequency is `145350000` Hz. @@ -72,9 +80,9 @@ The first parameter is the frequency in MHz, and the second optional parameter i ### Demodulate AM - 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 - + 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_s16 | 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) @@ -87,16 +95,16 @@ The first parameter is the frequency in MHz, and the second optional parameter i ### Demodulate SSB - 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 - + 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_s16 | 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,9 +115,9 @@ 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`) -- `i16` is `signed short` of 2 bytes/16 bits (e. g. sound card input is usually `i16`) +- `s16` is `signed short` of 2 bytes/16 bits (e. g. sound card input is usually `s16`) Functions usually end as: @@ -117,19 +125,21 @@ 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_i16_f` -- `csdr convert_f_i16` +- `csdr convert_s16_f` +- `csdr convert_f_s16` How to interpret: `csdr convert__` You can use these commands on complex streams, too, as they are only interleaved values (I,Q,I,Q,I,Q... coming after each other). +> Note: The the functions with `i16` in their names have been renamed, but still work (e.g. `csdr convert_f_i16`). + #### csdr commands `csdr` should be considered as a reference implementation on using `libcsdr`. For additional details on how to use the library, check `csdr.c` and `libcsdr.c`. @@ -157,13 +167,17 @@ It multiplies all samples by `gain`. It copies the input to the output. + through + +It copies the input to the output, while also displaying the speed of the data going through it. + none 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. @@ -173,7 +187,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 @@ -184,8 +198,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. @@ -257,7 +271,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 @@ -286,7 +300,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. @@ -312,7 +326,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. @@ -330,7 +344,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. @@ -363,11 +377,45 @@ The actual number of padding samples can be determined by running `cat csdr.c | It exchanges the first and second part of the FFT vector, to prepare it for the waterfall/spectrum display. It should operate on the data output from `logpower_cf`. + dsb_fc [q_value] + +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). +With `q_value = 0` it is an AM-DSB/SC modulator. If you want to get an AM-DSB signal, you will have to add a carrier to it. + + add_dcoffset_cc + +It adds a DC offset to the complex signal: `i_output = 0.5 + i_input / 2, q_output = q_input / 2` + + convert_f_samplerf + +It converts a real signal to the `-mRF` input format of [https://github.com/F5OEO/rpitx](rpitx), so it allows you to generate frequency modulation. The input signal will be the modulating signal. The `` parameter is the value for `rpitx` indicating the time to wait between samples. For a sampling rate of 48 ksps, this is 20833. + + fmmod_fc + +It generates a complex FM modulated output from a real input signal. + + fixed_amplitude_cc + +It changes the amplitude of every complex input sample to a fixed value. It does not change the phase information of the samples. + + mono2stereo_s16 + +It doubles every input sample. + setbuf -If the environment variable `CSDR_DYNAMIC_BUFSIZE_ON` is set to 1, then you can use this command to set the input buffer size for the next `csdr` process in the chain. +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. + + fifo + +It is similar to `clone`, but internally it uses a circular buffer. It reads as much as possible from the input. It discards input samples if the input buffer is full. #### Control via pipes @@ -388,9 +436,59 @@ 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 + +*csdr* has three modes of determining the buffer sizes, which can be chosen by the appropriate environment variables: +* *default:* 16k or 1k buffer is chosen based on function, +* *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**. +* 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. +* 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. +* The 8 bytes of information included in the beginning of the stream is: + * a preamble of the bytes 'c','s','d','r' (4 bytes), + * the buffer size stored as `int` (4 bytes). +* This size always counts as samples, as we expect that the user takes care of connecting the functions with right data types to each other. + +> I added this feature while researching how to decrease the latency of a DSP chain consisting of several multirate algorithms.
+> For example, a `csdr fir_decimate_cc 10` would use an input buffer of 10240, and an output buffer of 1024. The next process in the chain, `csdr bandpass_fir_fft_cc` would automatically adjust to it, using a buffer of 1024 for both input and output.
+> In contrast to original expectations, using dynamic buffer sizes didn't decrease the latency much. + +If dynamic buffer size determination is disabled, you can still set a **fixed buffer size** with `export CSDR_FIXED_BUFSIZE=`. + +For debug purposes, buffer sizes of all processes can be printed using `export CSDR_PRINT_BUFSIZES=1`. + +If you add your own functions to `csdr`, you have to initialize the buffers before doing the processing. Buffer size will be stored in the global variable `the_bufsize`. + +Example of initialization if the process generates N output samples for N input samples: + + if(!sendbufsize(initialize_buffers())) return -2; + +Example of initalization if the process generates N/D output samples for N input samples: + + if(!initialize_buffers()) return -2; + sendbufsize(the_bufsize/D); + +Example of initialization if the process allocates memory for itself, and it doesn't want to use the global buffers: + + getbufsize(); //dummy + sendbufsize(my_own_bufsize); + +Example of initialization if the process always works with a fixed output size, regardless of the input: + + if(!initialize_buffers()) return -2; + sendbufsize(fft_size); + #### Testbench `csdr` was tested with GNU Radio Companion flowgraphs. These flowgraphs are available under the directory `grc_tests`, and they require the gr-ha5kfu set of blocks for GNU Radio. diff --git a/csdr.c b/csdr.c index ee32400..5506ab8 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,20 +35,21 @@ 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" #include #include #include +#include #include "fastddc.h" char usage[]= @@ -60,8 +61,10 @@ char usage[]= " convert_f_u8\n" " convert_s8_f\n" " convert_f_s8\n" -" convert_f_i16\n" -" convert_i16_f\n" +" convert_f_s16\n" +" convert_s16_f\n" +" convert_f_s24 [--bigendian]\n" +" convert_s24_f [--bigendian]\n" " realpart_cf\n" " clipdetect_ff\n" " limit_ff [max_amplitude]\n" @@ -73,7 +76,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" @@ -97,10 +102,21 @@ char usage[]= " logpower_cf [add_db]\n" " fft_benchmark [--benchmark]\n" " bandpass_fir_fft_cc [window]\n" -" encode_ima_adpcm_i16_u8\n" -" decode_ima_adpcm_u8_i16\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" +" flowcontrol \n" +" through\n" +" dsb_fc [q_value]\n" +" convert_f_samperf \n" +" fmmod_fc\n" +" fixed_amplitude_cc \n" +" monos2stereo_s16\n" +" setbuf \n" " fft_exchange_sides_ff \n" +" squelch_and_smeter_cc --fifo --outfifo \n" +" fifo \n" " \n" ; @@ -202,12 +218,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 STDIN_FILENO) ? STDOUT_FILENO : STDIN_FILENO) + 1; + + int fifo_actual_buffer_wr = fifo_num_buffers - 1; + int fifo_actual_buffer_rd = 0; + int fifo_actual_buffer_wr_pos = 0; + int fifo_actual_buffer_rd_pos = 0; + int fifo_error = 0; + int fifo_overrun_shown = 0; + + for(;;) + { + select(highfd, &read_fds, NULL, NULL, NULL); + + //try to read until buffer is full + if(FD_ISSET(STDIN_FILENO, &read_fds)) for(;;) + { + int read_bytes=read(STDIN_FILENO, fifo_buffers[fifo_actual_buffer_rd]+fifo_actual_buffer_rd_pos, fifo_buffer_size-fifo_actual_buffer_rd_pos); + //fprintf(stderr, "r %d %d | %d %d\n", read_bytes, fifo_buffer_size-fifo_actual_buffer_rd_pos, fifo_actual_buffer_rd, fifo_actual_buffer_rd_pos); + if(!read_bytes || ((read_bytes<0)&&(fifo_error=read_bytes)) ) break; + fifo_actual_buffer_rd_pos+=read_bytes; + if(!((fifo_actual_buffer_rd==fifo_actual_buffer_wr-1)||(fifo_actual_buffer_wr==0&&fifo_actual_buffer_rd==fifo_num_buffers-1))) + { + if(fifo_actual_buffer_rd_pos==fifo_buffer_size) + { + fifo_overrun_shown = 0; + fifo_actual_buffer_rd++; + fifo_actual_buffer_rd_pos = 0; + if(fifo_actual_buffer_rd>=fifo_num_buffers) fifo_actual_buffer_rd=0; + } + } + else + { + if(fifo_actual_buffer_rd_pos==fifo_buffer_size) + { + fifo_actual_buffer_rd_pos = 0; //rewrite same buffer + if(!fifo_overrun_shown) { fifo_overrun_shown=1; fprintf(stderr, "fifo: circular buffer full, dropping samples\n"); } + } + } + } + //try to write until buffer is empty + if(FD_ISSET(STDOUT_FILENO, &write_fds)) for(;;) + { + if(fifo_actual_buffer_wr == fifo_actual_buffer_rd) break; + int written_bytes=write(STDOUT_FILENO, fifo_buffers[fifo_actual_buffer_wr]+fifo_actual_buffer_wr_pos, fifo_buffer_size-fifo_actual_buffer_wr_pos); + //fprintf(stderr, "w %d %d | %d %d\n", written_bytes, fifo_buffer_size-fifo_actual_buffer_wr_pos, fifo_actual_buffer_wr, fifo_actual_buffer_wr_pos); + if(!written_bytes || ((written_bytes<0)&&(fifo_error=written_bytes)) ) break; + fifo_actual_buffer_wr_pos+=written_bytes; + if(fifo_actual_buffer_wr_pos==fifo_buffer_size) + { + fifo_actual_buffer_wr++; + fifo_actual_buffer_wr_pos = 0; + if(fifo_actual_buffer_wr>=fifo_num_buffers) fifo_actual_buffer_wr=0; + } + + } + if(fifo_error&&errno!=11) { fprintf(stderr,"fifo: fifo_error (%d)", errno); return -1; } + } + + return -1; + + } + if(!strcmp(argv[1],"convert_u8_f")) { @@ -368,27 +472,29 @@ int main(int argc, char *argv[]) } if(!strcmp(argv[1],"convert_s8_f")) { + if(!sendbufsize(initialize_buffers())) return -2; for(;;) { FEOF_CHECK; - fread(buffer_s8, sizeof(signed char), BUFSIZE, stdin); - convert_s8_f(buffer_s8, output_buffer, BUFSIZE); + fread((signed char*)buffer_u8, sizeof(signed char), the_bufsize, stdin); + convert_s8_f((signed char*)buffer_u8, output_buffer, the_bufsize); FWRITE_R; TRY_YIELD; } } if(!strcmp(argv[1],"convert_f_s8")) //not tested { + if(!sendbufsize(initialize_buffers())) return -2; for(;;) { FEOF_CHECK; FREAD_R; - convert_f_s8(input_buffer, buffer_s8, BUFSIZE); - fwrite(buffer_s8, sizeof(signed char), BUFSIZE, stdout); + convert_f_s8(input_buffer, (signed char*)buffer_u8, the_bufsize); + fwrite((signed char*)buffer_u8, sizeof(signed char), the_bufsize, stdout); TRY_YIELD; } } - if(!strcmp(argv[1],"convert_f_i16")) + if((!strcmp(argv[1],"convert_f_i16")) || (!strcmp(argv[1],"convert_f_s16"))) { if(!sendbufsize(initialize_buffers())) return -2; for(;;) @@ -400,7 +506,7 @@ int main(int argc, char *argv[]) TRY_YIELD; } } - if(!strcmp(argv[1],"convert_i16_f")) //not tested + if((!strcmp(argv[1],"convert_i16_f")) || (!strcmp(argv[1],"convert_s16_f"))) { if(!sendbufsize(initialize_buffers())) return -2; for(;;) @@ -412,6 +518,34 @@ int main(int argc, char *argv[]) TRY_YIELD; } } + if(!strcmp(argv[1],"convert_f_s24")) + { + int bigendian = (argc>2) && (!strcmp(argv[2],"--bigendian")); + unsigned char* s24buffer = (unsigned char*)malloc(sizeof(unsigned char)*the_bufsize*3); + if(!sendbufsize(initialize_buffers())) return -2; + for(;;) + { + FEOF_CHECK; + FREAD_R; + convert_f_s24(input_buffer, s24buffer, the_bufsize, bigendian); + fwrite(s24buffer, sizeof(unsigned char)*3, the_bufsize, stdout); + TRY_YIELD; + } + } + if(!strcmp(argv[1],"convert_s24_f")) + { + int bigendian = (argc>2) && (!strcmp(argv[2],"--bigendian")); + unsigned char* s24buffer = (unsigned char*)malloc(sizeof(unsigned char)*the_bufsize*3); + if(!sendbufsize(initialize_buffers())) return -2; + for(;;) + { + FEOF_CHECK; + fread(s24buffer, sizeof(unsigned char)*3, the_bufsize, stdin); + convert_s24_f(s24buffer, output_buffer, the_bufsize, bigendian); + FWRITE_R; + TRY_YIELD; + } + } if(!strcmp(argv[1],"realpart_cf")) { if(!sendbufsize(initialize_buffers())) return -2; @@ -438,9 +572,9 @@ int main(int argc, char *argv[]) } if(!strcmp(argv[1],"gain_ff")) { - if(argc<=2) return badsyntax("need required parameter (gain)"); + if(argc<=2) return badsyntax("need required parameter (gain)"); float gain; - sscanf(argv[2],"%g",&gain); + sscanf(argv[2],"%g",&gain); if(!sendbufsize(initialize_buffers())) return -2; for(;;) { @@ -467,23 +601,23 @@ int main(int argc, char *argv[]) } if(!strcmp(argv[1],"yes_f")) { - if(argc<=2) return badsyntax("need required parameter (to_repeat)"); + if(argc<=2) return badsyntax("need required parameter (to_repeat)"); float to_repeat; sscanf(argv[2],"%g",&to_repeat); int buf_times = 0; if(argc>=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 @@ -506,7 +640,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; @@ -629,7 +763,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; @@ -666,7 +800,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); } @@ -703,7 +837,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 @@ -715,7 +849,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; @@ -790,7 +924,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); @@ -815,19 +949,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; @@ -1104,7 +1242,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); @@ -1113,7 +1251,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; } @@ -1122,11 +1260,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; @@ -1158,7 +1296,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(;;) @@ -1179,7 +1317,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); @@ -1194,7 +1332,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). @@ -1202,7 +1340,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 @@ -1220,10 +1358,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; @@ -1233,12 +1371,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]); @@ -1255,6 +1393,8 @@ int main(int argc, char *argv[]) FFT_PLAN_T* plan=make_fft_c2c(fft_size, windowed, output, 1, benchmark); if(benchmark) fprintf(stderr," done\n"); if(octave) printf("setenv(\"GNUTERM\",\"X11 noraise\");y=zeros(1,%d);semilogy(y,\"ydatasource\",\"y\");\n",fft_size); + float *windowt; + windowt = precalculate_window(fft_size, window); for(;;) { FEOF_CHECK; @@ -1273,14 +1413,15 @@ int main(int argc, char *argv[]) for(int i=0;i=3) sscanf(argv[2],"%g",&add_db); - if(!sendbufsize(initialize_buffers())) return -2; + if(!sendbufsize(initialize_buffers())) return -2; for(;;) { @@ -1309,9 +1450,43 @@ int main(int argc, char *argv[]) } } + if(!strcmp(argv[1],"logaveragepower_cf")) + { + bigbufs=1; + if(argc<=4) return badsyntax("need required parameters (add_db, table_size, avgnumber)"); + float add_db=0; + int avgnumber=0; + int fft_size=0; + + sscanf(argv[2],"%g",&add_db); + sscanf(argv[3],"%d",&fft_size); + sscanf(argv[4],"%d",&avgnumber); + + float *input = malloc(sizeof(float)*2 * fft_size); + float *output = malloc(sizeof(float) * fft_size); + + add_db -= 10.0*log10(avgnumber); + for(;;) + { + int i,n; + for(i = 0; i < fft_size; i++) { + output[i] = 0; + } + FEOF_CHECK; + for(n = 0; n < avgnumber; n++) { + fread (input, sizeof(float)*2, fft_size, stdin); + accumulate_power_cf((complexf*)input, output, fft_size); + } + log_ff(output, output, fft_size, add_db); + fwrite (output, sizeof(float), fft_size, stdout); + TRY_YIELD; + } + return 0; + } + 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 @@ -1334,13 +1509,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; @@ -1366,12 +1541,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); @@ -1380,20 +1555,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) @@ -1638,10 +1812,11 @@ int main(int argc, char *argv[]) TRY_YIELD; } } +#endif 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; @@ -1649,21 +1824,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))) { @@ -1687,20 +1862,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],"fastddc_fwd_cc") ) // [transition_bw [window]] @@ -1942,9 +2170,5 @@ int main(int argc, char *argv[]) return 0; } - fprintf(stderr,"csdr: function name given in argument 1 (%s) does not exist. Possible causes:\n- You mistyped the commandline.\n- You need to update csdr to a newer version (if available).", argv[1]); - return -1; - + fprintf(stderr,"csdr: function name given in argument 1 (%s) does not exist. Possible causes:\n- You mistyped the commandline.\n- You need to update csdr to a newer version (if available).", argv[1]); return -1; } - - diff --git a/libcsdr.c b/libcsdr.c index 13203b8..972d216 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 { @@ -480,8 +480,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]); @@ -489,10 +489,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) */ @@ -500,7 +500,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) @@ -511,13 +511,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]; @@ -581,7 +581,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; @@ -612,7 +612,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: @@ -657,13 +657,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; @@ -674,7 +674,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; @@ -716,35 +716,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; @@ -868,9 +868,9 @@ void fastagc_ff(fastagc_ff_t* input, float* output) } /* - ______ __ __ _ _ _ _ - | ____| \/ | | | | | | | | | - | |__ | \ / | __| | ___ _ __ ___ ___ __| |_ _| | __ _| |_ ___ _ __ ___ + ______ __ __ _ _ _ _ + | ____| \/ | | | | | | | | | + | |__ | \ / | __| | ___ _ __ ___ ___ __| |_ _| | __ _| |_ ___ _ __ ___ | __| | |\/| | / _` |/ _ \ '_ ` _ \ / _ \ / _` | | | | |/ _` | __/ _ \| '__/ __| | | | | | | | (_| | __/ | | | | | (_) | (_| | |_| | | (_| | || (_) | | \__ \ |_| |_| |_| \__,_|\___|_| |_| |_|\___/ \__,_|\__,_|_|\__,_|\__\___/|_| |___/ @@ -921,33 +921,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 @@ -967,7 +967,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) @@ -1018,20 +1018,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 { @@ -1092,7 +1112,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 } @@ -1166,16 +1222,56 @@ void convert_f_s8(float* input, signed char* output, int input_size) for(int i=0;i1.0) input[i]=1.0; if(input[i]<-1.0) input[i]=-1.0; }*/ - for(int i=0;i>8); + unsigned char* ptemp=(unsigned char*)&temp; + output[k++]=*ptemp; + output[k++]=*(ptemp+1); + output[k++]=*(ptemp+2); + } + else for(int i=0;i>8); + unsigned char* ptemp=(unsigned char*)&temp; + output[k++]=*(ptemp+2); + output[k++]=*(ptemp+1); + output[k++]=*ptemp; + } +} + +void convert_s24_f(unsigned char* input, float* output, int input_size, int bigendian) +{ + int k=0; + if(bigendian) for(int i=0;i @@ -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; @@ -66,7 +66,7 @@ typedef struct complexf_s { float i; float q; } complexf; #define TIME_TAKEN(start,end) ((end.tv_sec-start.tv_sec)+(end.tv_nsec-start.tv_nsec)/1e9) //window -typedef enum window_s +typedef enum window_s { WINDOW_BOXCAR, WINDOW_BLACKMAN, WINDOW_HAMMING } window_t; @@ -115,7 +115,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; @@ -137,9 +137,13 @@ typedef struct rational_resampler_ff_s 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); void rational_resampler_get_lowpass_f(float* output, int output_size, int interpolation, int decimation, window_t window); +float *precalculate_window(int size, window_t window); void apply_window_c(complexf* input, complexf* output, int size, window_t window); +void apply_precalculated_window_c(complexf* input, complexf* output, int size, float *windowt); void apply_window_f(float* input, float* output, int size, window_t window); void logpower_cf(complexf* input, float* output, int size, float add_db); +void accumulate_power_cf(complexf* input, float* output, int size); +void log_ff(float* input, float* output, int size, float add_db); typedef struct fractional_decimator_ff_s { @@ -152,7 +156,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); @@ -182,6 +186,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); @@ -191,7 +197,12 @@ void convert_u8_f(unsigned char* input, float* output, int input_size); void convert_f_u8(float* input, unsigned char* output, int input_size); void convert_s8_f(signed char* input, float* output, int input_size); void convert_f_s8(float* input, signed char* output, int input_size); +void convert_f_s16(float* input, short* output, int input_size); +void convert_s16_f(short* input, float* output, int input_size); void convert_f_i16(float* input, short* output, int input_size); void convert_i16_f(short* input, float* output, int input_size); +void convert_f_s24(float* input, unsigned char* output, int input_size, int bigendian); +void convert_s24_f(unsigned char* input, float* output, int input_size, int bigendian); + int is_nan(float f); diff --git a/libcsdr_gpl.c b/libcsdr_gpl.c index 738a3b5..b895e00 100644 --- a/libcsdr_gpl.c +++ b/libcsdr_gpl.c @@ -47,7 +47,7 @@ float shift_addition_cc(complexf *input, complexf* output, int input_size, shift } starting_phase+=d.rate*PI*input_size; while(starting_phase>PI) starting_phase-=2*PI; //@shift_addition_cc: normalize starting_phase - while(starting_phase<-PI) starting_phase+=2*PI; + while(starting_phase<-PI) starting_phase+=2*PI; return starting_phase; } @@ -82,12 +82,12 @@ void shift_addition_cc_test(shift_addition_data_t d) sinphi=sinphi_last*d.cosdelta+cosphi_last*d.sindelta; phi+=d.rate*PI; while(phi>2*PI) phi-=2*PI; //@shift_addition_cc: normalize phase - if(i%SACCTEST_STEP==0) + if(i%SACCTEST_STEP==0) { avg_counter=avg_size; avg=0; } - if(avg_counter) + if(avg_counter) { avg+=fabs(cosphi-cos(phi)); if(!--avg_counter) printf("%g ", avg/avg_size); @@ -128,7 +128,7 @@ decimating_shift_addition_status_t decimating_shift_addition_cc(complexf *input, s.starting_phase+=d.rate*PI*k; s.output_size=k; while(s.starting_phase>PI) s.starting_phase-=2*PI; //@shift_addition_cc: normalize starting_phase - while(s.starting_phase<-PI) s.starting_phase+=2*PI; + while(s.starting_phase<-PI) s.starting_phase+=2*PI; return s; } @@ -142,10 +142,10 @@ float agc_ff(float* input, float* output, int input_size, float reference, float hang_time = (hang_time_ms / 1000) * sample_rate hang_time is given in samples, and should be about 4ms. hang_time can be switched off by setting it to zero (not recommended). - + max_gain = pow(2, adc_bits) - max_gain should be no more than the dynamic range of your A/D converter. - gain_filter_alpha = 1 / ((fs/(2*PI*fc))+1) + max_gain should be no more than the dynamic range of your A/D converter. + gain_filter_alpha = 1 / ((fs/(2*PI*fc))+1) >>> 1 / ((48000./(2*3.141592654*100))+1) 0.012920836043344543 @@ -153,7 +153,7 @@ float agc_ff(float* input, float* output, int input_size, float reference, float 0.0013072857061786625 - Literature: + Literature: ww.qsl.net/va3iul/Files/Automatic_Gain_Control.pdf page 7 of http://www.arrl.org/files/file/Technology/tis/info/pdf/021112qex027.pdf @@ -170,10 +170,10 @@ float agc_ff(float* input, float* output, int input_size, float reference, float output[0]=last_gain*input[0]; //we skip this one sample, because it is easier this way for(int i=1;imax_gain) gain=max_gain; //We also have to limit our gain, it can't be infinity. - if(gain<0) gain=0; } + if(gain>max_gain) gain=max_gain; //We also have to limit our gain, it can't be infinity. + if(gain<0) gain=0; //output[i]=gain*input[i]; //Here we do the actual scaling of the samples. //Here we do the actual scaling of the samples, but we run an IIR filter on the gain values: output[i]=(gain=gain+last_gain-gain_filter_alpha*last_gain)*input[i]; //dc-pass-filter: freqz([1 -1],[1 -0.99]) y[i]=x[i]+y[i-1]-alpha*x[i-1] @@ -230,8 +229,7 @@ float agc_ff(float* input, float* output, int input_size, float reference, float last_gain=gain; } - return gain; //this will be the last_gain next time + return gain; //this will be the last_gain next time } #endif -