sdr (empty) → 0.1.0.0
raw patch · 23 files changed
+2804/−0 lines, 23 filesdep +Chartdep +Chart-cairodep +Decimalsetup-changed
Dependencies added: Chart, Chart-cairo, Decimal, GLFW-b, OpenGL, QuickCheck, array, base, bytestring, cairo, cereal, colour, containers, criterion, dynamic-graph, either, fftwRaw, optparse-applicative, pango, pipes, pipes-bytestring, pipes-concurrency, primitive, pulse-simple, rtlsdr, sdr, storable-complex, test-framework, test-framework-quickcheck2, time, tuple, vector
Files
- LICENSE +30/−0
- Setup.hs +2/−0
- benchmarks/Benchmarks.hs +127/−0
- c_sources/convert.c +50/−0
- c_sources/decimate.c +147/−0
- c_sources/filter.c +162/−0
- c_sources/resample.c +87/−0
- c_sources/scale.c +38/−0
- hs_sources/SDR/ArgUtils.hs +19/−0
- hs_sources/SDR/Demod.hs +43/−0
- hs_sources/SDR/FFT.hs +162/−0
- hs_sources/SDR/Filter.hs +496/−0
- hs_sources/SDR/FilterDesign.hs +76/−0
- hs_sources/SDR/FilterInternal.hs +401/−0
- hs_sources/SDR/PipeUtils.hs +55/−0
- hs_sources/SDR/Plot.hs +172/−0
- hs_sources/SDR/Pulse.hs +33/−0
- hs_sources/SDR/RTLSDRStream.hs +53/−0
- hs_sources/SDR/Serialize.hs +84/−0
- hs_sources/SDR/Util.hs +154/−0
- hs_sources/SDR/VectorUtils.hs +64/−0
- sdr.cabal +136/−0
- tests/TestSuite.hs +213/−0
+ LICENSE view
@@ -0,0 +1,30 @@+Copyright (c) 2014, Adam Walker++All rights reserved.++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:++ * Redistributions of source code must retain the above copyright+ notice, this list of conditions and the following disclaimer.++ * Redistributions in binary form must reproduce the above+ copyright notice, this list of conditions and the following+ disclaimer in the documentation and/or other materials provided+ with the distribution.++ * Neither the name of Adam Walker nor the names of other+ contributors may be used to endorse or promote products derived+ from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ Setup.hs view
@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain
+ benchmarks/Benchmarks.hs view
@@ -0,0 +1,127 @@+{-# LANGUAGE ScopedTypeVariables, BangPatterns #-}++import Control.Monad.Primitive +import Control.Monad+import Foreign.C.Types+import Foreign.Ptr+import Unsafe.Coerce+import Data.Complex++import qualified Data.Vector.Generic as VG+import qualified Data.Vector.Generic.Mutable as VGM+import qualified Data.Vector.Storable as VS+import qualified Data.Vector.Storable.Mutable as VSM+import qualified Data.Vector.Fusion.Stream as VFS+import qualified Data.Vector.Fusion.Stream.Monadic as VFSM++import Foreign.Storable.Complex+import Criterion.Main++import SDR.FilterInternal+import SDR.Util++theBench :: IO ()+theBench = do+ --Setup+ let size = 16384+ numCoeffs = 128+ num = size - numCoeffs + 1+ decimation = 4+ interpolation = 3+ numCoeffsDiv2 = 64++ coeffsList :: [Float]+ coeffsList = take numCoeffs [0 ..]+ coeffs :: VS.Vector Float+ coeffs = VG.fromList $ take numCoeffs [0 ..]+ coeffsSym :: VS.Vector Float+ coeffsSym = VG.fromList $ take numCoeffsDiv2 [0 ..]+ inBuf :: VS.Vector Float+ inBuf = VG.fromList $ take size [0 ..]+ inBufComplex :: VS.Vector (Complex Float)+ inBufComplex = VG.fromList $ take size $ do+ i <- [0..]+ return $ i :+ i++ numConv = 16386+ inBufConv :: VS.Vector CUChar+ inBufConv = VG.fromList $ take size $ concat $ repeat [0 .. 255]++ duplicate :: [a] -> [a]+ duplicate = concat . map func + where func x = [x, x]++ coeffs2 :: VS.Vector Float+ coeffs2 = VG.fromList $ duplicate $ take numCoeffs [0 ..]++ outBuf :: VS.MVector RealWorld Float <- VGM.new size+ outBufComplex :: VS.MVector RealWorld (Complex Float) <- VGM.new size++ resampler3 <- resampleCRR2 interpolation decimation coeffsList+ resampler4 <- resampleCSSERR interpolation decimation coeffsList+ resampler5 <- resampleCAVXRR interpolation decimation coeffsList++ --Benchmarks+ defaultMain [+ bgroup "filter" [+ bgroup "real" [+ bench "highLevel" $ nfIO $ filterHighLevel coeffs num inBuf outBuf,+ bench "imperative1" $ nfIO $ filterImperative1 coeffs num inBuf outBuf,+ bench "imperative2" $ nfIO $ filterImperative2 coeffs num inBuf outBuf,+ bench "c" $ nfIO $ filterCRR coeffs num inBuf outBuf,+ bench "cSSE" $ nfIO $ filterCSSERR coeffs num inBuf outBuf,+ bench "cSSESym" $ nfIO $ filterCSSESymmetricRR coeffsSym num inBuf outBuf,+ bench "cAVX" $ nfIO $ filterCAVXRR coeffs num inBuf outBuf,+ bench "cAVXSym" $ nfIO $ filterCAVXSymmetricRR coeffsSym num inBuf outBuf+ ],+ bgroup "complex" [+ bench "highLevel" $ nfIO $ filterHighLevel coeffs num inBufComplex outBufComplex,+ bench "c" $ nfIO $ filterCRC coeffs num inBufComplex outBufComplex,+ bench "cSSE" $ nfIO $ filterCSSERC coeffs2 num inBufComplex outBufComplex,+ bench "cSSE2" $ nfIO $ filterCSSERC2 coeffs num inBufComplex outBufComplex,+ bench "cSSESym" $ nfIO $ filterCSSESymmetricRC coeffsSym num inBufComplex outBufComplex,+ bench "cAVX" $ nfIO $ filterCAVXRC coeffs2 num inBufComplex outBufComplex,+ bench "cAVX2" $ nfIO $ filterCAVXRC2 coeffs num inBufComplex outBufComplex,+ bench "cAVXSym" $ nfIO $ filterCAVXSymmetricRC coeffsSym num inBufComplex outBufComplex+ ]+ ],+ bgroup "decimate" [+ bgroup "real" [+ bench "highLevel" $ nfIO $ decimateHighLevel decimation coeffs (num `quot` decimation) inBuf outBuf,+ bench "c" $ nfIO $ decimateCRR decimation coeffs (num `quot` decimation) inBuf outBuf,+ bench "cSSE" $ nfIO $ decimateCSSERR decimation coeffs (num `quot` decimation) inBuf outBuf,+ bench "cSSESym" $ nfIO $ decimateCSSESymmetricRR decimation coeffsSym (num `quot` decimation) inBuf outBuf,+ bench "cAVX" $ nfIO $ decimateCAVXRR decimation coeffs (num `quot` decimation) inBuf outBuf,+ bench "cAVXSym" $ nfIO $ decimateCAVXSymmetricRR decimation coeffsSym (num `quot` decimation) inBuf outBuf+ ],+ bgroup "complex" [+ bench "highLevel" $ nfIO $ decimateHighLevel decimation coeffs (num `quot` decimation) inBufComplex outBufComplex,+ bench "c" $ nfIO $ decimateCRC decimation coeffs (num `quot` decimation) inBufComplex outBufComplex,+ bench "cSSE" $ nfIO $ decimateCSSERC decimation coeffs2 (num `quot` decimation) inBufComplex outBufComplex,+ bench "cSSE2" $ nfIO $ decimateCSSERC2 decimation coeffs (num `quot` decimation) inBufComplex outBufComplex,+ bench "cSSESym" $ nfIO $ decimateCSSESymmetricRC decimation coeffsSym (num `quot` decimation) inBufComplex outBufComplex,+ bench "cAVX" $ nfIO $ decimateCAVXRC decimation coeffs2 (num `quot` decimation) inBufComplex outBufComplex,+ bench "cAVX2" $ nfIO $ decimateCAVXRC2 decimation coeffs (num `quot` decimation) inBufComplex outBufComplex,+ bench "cAVXSym" $ nfIO $ decimateCAVXSymmetricRC decimation coeffsSym (num `quot` decimation) inBufComplex outBufComplex+ ]+ ],+ bgroup "resample" [+ bgroup "real" [+ bench "highLevel" $ nfIO $ resampleHighLevel interpolation decimation coeffs 0 (num `quot` decimation) inBuf outBuf,+ bench "c" $ nfIO $ resampleCRR (num `quot` decimation) interpolation decimation 0 coeffs inBuf outBuf,+ bench "c2" $ nfIO $ resampler3 (num `quot` decimation) 0 inBuf outBuf,+ bench "cSSE" $ nfIO $ resampler4 (num `quot` decimation) 0 inBuf outBuf,+ bench "cAVX" $ nfIO $ resampler5 (num `quot` decimation) 0 inBuf outBuf+ ],+ bgroup "complex" [+ bench "highLevel" $ nfIO $ resampleHighLevel interpolation decimation coeffs 0 (num `quot` decimation) inBufComplex outBufComplex+ ]+ ],+ bgroup "scaling" [+ bench "c" $ nfIO $ scaleC size 0.3 inBuf outBuf,+ bench "cSSE" $ nfIO $ scaleCSSE size 0.3 inBuf outBuf,+ bench "cAVX" $ nfIO $ scaleCAVX size 0.3 inBuf outBuf+ ]+ ]++main = theBench
+ c_sources/convert.c view
@@ -0,0 +1,50 @@+/*+ * Conversion of byte IQ data from RTLSDR device to complex floats.+ * These exist because the pure Haskell implementations are slow.+ * Uses SIMD instructions for performance.+ */++#include <stdio.h>+#include <stdint.h>+#include <x86intrin.h>++/*+ * Conversion+ */++void convertC(int num, uint8_t *in, float *out){+ int i;+ for(i=0; i<num; i++){+ out[i] = ((float) in[i] - 128.0) * (1/128.0);+ }+}++void convertCSSE(int num, uint8_t *in, float *out){+ int i;+ __m128 sub = _mm_set1_ps(128.0);+ __m128 mul = _mm_set1_ps(1/128.0);+ for(i=0; i<num; i+=4){+ __m128i val = _mm_loadu_si128((__m128i *)(in + i));+ __m128i ints = _mm_cvtepu8_epi32(val);+ __m128 cvtd = _mm_cvtepi32_ps(ints);++ __m128 res = _mm_mul_ps(_mm_sub_ps(cvtd, sub), mul);++ _mm_storeu_ps(out + i, res);+ }+}++void convertCAVX(int num, uint8_t *in, float *out){+ int i;+ __m256 sub = _mm256_set1_ps(128.0);+ __m256 mul = _mm256_set1_ps(1/128.0);+ for(i=0; i<num; i+=8){+ __m128i val = _mm_loadu_si128((__m128i *)(in + i));+ __m256i ints = _mm256_cvtepu8_epi32(val);+ __m256 cvtd = _mm256_cvtepi32_ps(ints);++ __m256 res = _mm256_mul_ps(_mm256_sub_ps(cvtd, sub), mul);++ _mm256_storeu_ps(out + i, res);+ }+}
+ c_sources/decimate.c view
@@ -0,0 +1,147 @@+/*+ * FIR decimation of complex and real data with real coefficients.+ * These exist because the pure Haskell implementations are slow.+ * Uses SIMD instructions for performance.+ */++#include <stdio.h>+#include <stdint.h>+#include <x86intrin.h>++#include "common.h"++/*+ * Real coefficients, real inputs+ */+void decimateRR(int num, int factor, int numCoeffs, float *coeffs, float *inBuf, float *outBuf){+ int i, k;+ for(i=0, k=0; i<num; i++, k+=factor){+ float *startPtr = inBuf + k;+ outBuf[i] = dotprod_R(numCoeffs, coeffs, startPtr);+ }+}++/*+ * SIMD versions+ */+void decimateSSERR(int num, int factor, int numCoeffs, float *coeffs, float *inBuf, float *outBuf){+ int i, k;+ for(i=0, k=0; i<num; i++, k+=factor){+ float *startPtr = inBuf + k;+ __m128 accum = sse_dotprod_R(numCoeffs, coeffs, startPtr);+ accum = sse_hadd_R(accum);+ _mm_store_ss(outBuf + i, accum);+ }+}++void decimateAVXRR(int num, int factor, int numCoeffs, float *coeffs, float *inBuf, float *outBuf){+ int i, k;+ for(i=0, k=0; i<num; i++, k+=factor){+ float *startPtr = inBuf + k;+ __m256 accum = avx_dotprod_R(numCoeffs, coeffs, startPtr);+ __m128 accum2 = avx_hadd_R(accum);+ _mm_store_ss(outBuf + i, accum2);+ }+}++/*+ * Symmetric versions+ */+void decimateSSESymmetricRR(int num, int factor, int numCoeffs, float *coeffs, float *inBuf, float *outBuf){+ int i, k;+ for(i=0, k=0; i<num; i++, k+=factor){+ float *startPtr = inBuf + k;+ __m128 accum = sse_sym_dotprod_R(numCoeffs, coeffs, startPtr);+ accum = sse_hadd_R(accum);+ _mm_store_ss(outBuf + i, accum);+ }+}++void decimateAVXSymmetricRR(int num, int factor, int numCoeffs, float *coeffs, float *inBuf, float *outBuf){+ int i, j, k;+ for(i=0, k=0; i<num; i++, k+=factor){+ float *startPtr = inBuf + k;+ __m256 accum = avx_sym_dotprod_R(numCoeffs, coeffs, startPtr);+ __m128 accum2 = avx_hadd_R(accum);+ _mm_store_ss(outBuf + i, accum2);+ }+}++/*+ * Real coefficients, complex inputs+ */+void decimateRC(int num, int factor, int numCoeffs, float *coeffs, float *inBuf, float *outBuf){+ int i, k;+ for(i=0, k=0; i<num*2; i+=2, k+=factor*2){+ float *startPtr = inBuf + k;+ dotprod_C(numCoeffs, coeffs, startPtr, outBuf + i);+ }+}++/*+ * SIMD versions+ */+void decimateSSERC(int num, int factor, int numCoeffs, float *coeffs, float *inBuf, float *outBuf){+ int i, k;+ for(i=0, k=0; i<num*2; i+=2, k+=factor*2){+ float *startPtr = inBuf + k;+ __m128 accum = sse_dotprod_R(numCoeffs, coeffs, startPtr);+ accum = sse_hadd_C(accum);+ store_complex(outBuf + i, accum);+ }+}++void decimateSSERC2(int num, int factor, int numCoeffs, float *coeffs, float *inBuf, float *outBuf){+ int i, j, k;+ for(i=0, k=0; i<num*2; i+=2, k+=factor*2){+ float *startPtr = inBuf + k;+ __m128 accum = sse_dotprod_C(numCoeffs, coeffs, startPtr);+ accum = sse_hadd_C(accum);+ store_complex(outBuf + i, accum);+ }+}+++void decimateAVXRC(int num, int factor, int numCoeffs, float *coeffs, float *inBuf, float *outBuf){+ int i, j, k;+ for(i=0, k=0; i<num*2; i+=2, k+=factor*2){+ float *startPtr = inBuf + k;+ __m256 accum = avx_dotprod_R(numCoeffs, coeffs, startPtr);+ __m128 accum2 = avx_hadd_C(accum);+ store_complex(outBuf + i, accum2);+ }+}++void decimateAVXRC2(int num, int factor, int numCoeffs, float *coeffs, float *inBuf, float *outBuf){+ int i, j, k;+ for(i=0, k=0; i<num*2; i+=2, k+=factor*2){+ float *startPtr = inBuf + k;+ __m256 accum = avx_dotprod_C(numCoeffs, coeffs, startPtr);+ __m128 accum2 = avx_hadd_C(accum);+ store_complex(outBuf + i, accum2);+ }+}++/*+ * Symmetric versions+ */+void decimateSSESymmetricRC(int num, int factor, int numCoeffs, float *coeffs, float *inBuf, float *outBuf){+ int i, k;+ for(i=0, k=0; i<num*2; i+=2, k+=factor*2){+ float *startPtr = inBuf + k;+ __m128 accum = sse_sym_dotprod_C(numCoeffs, coeffs, startPtr);+ accum = sse_hadd_C(accum);+ store_complex(outBuf + i, accum);+ }+}++void decimateAVXSymmetricRC(int num, int factor, int numCoeffs, float *coeffs, float *inBuf, float *outBuf){+ int i, k;+ for(i=0, k=0; i<num*2; i+=2, k+=factor*2){+ float *startPtr = inBuf + k;+ __m256 accum = avx_sym_dotprod_C(numCoeffs, coeffs, startPtr);+ __m128 accum1 = avx_hadd_C(accum);+ store_complex(outBuf + i, accum1);+ }+}+
+ c_sources/filter.c view
@@ -0,0 +1,162 @@+/*+ * FIR filtering of complex and real data with real coefficients.+ * These exist because the pure Haskell implementations are slow.+ * Uses SIMD instructions for performance.+ */++#include <stdio.h>+#include <stdint.h>+#include <x86intrin.h>++#include "common.h"++/*+ * Real coefficients, real inputs+ */+void filterRR(int num, int numCoeffs, float *coeffs, float *inBuf, float *outBuf){+ int i;+ for(i=0; i<num; i++){+ float *startPtr = inBuf + i;+ outBuf[i] = dotprod_R(numCoeffs, coeffs, startPtr);+ }+}++/*+ * SIMD versions+ */+void filterSSERR(int num, int numCoeffs, float *coeffs, float *inBuf, float *outBuf){+ int i;+ for(i=0; i<num; i++){+ float *startPtr = inBuf + i;+ __m128 accum = sse_dotprod_R(numCoeffs, coeffs, startPtr);+ accum = sse_hadd_R(accum);+ _mm_store_ss(outBuf + i, accum);+ }+}++void filterAVXRR(int num, int numCoeffs, float *coeffs, float *inBuf, float *outBuf){+ int i;+ for(i=0; i<num; i++){+ float *startPtr = inBuf + i;+ __m256 accum = avx_dotprod_R(numCoeffs, coeffs, startPtr);+ __m128 accum2 = avx_hadd_R(accum);+ _mm_store_ss(outBuf + i, accum2);+ }+}++/*+ * Symmetric versions+ */+void filterSSESymmetricRR(int num, int numCoeffs, float *coeffs, float *inBuf, float *outBuf){+ int i;+ for(i=0; i<num; i++){+ float *startPtr = inBuf + i;+ __m128 accum = sse_sym_dotprod_R(numCoeffs, coeffs, startPtr);+ accum = sse_hadd_R(accum);+ _mm_store_ss(outBuf + i, accum);+ }+}++void filterAVXSymmetricRR(int num, int numCoeffs, float *coeffs, float *inBuf, float *outBuf){+ int i;+ for(i=0; i<num; i++){+ float *startPtr = inBuf + i;+ __m256 accum = avx_sym_dotprod_R(numCoeffs, coeffs, startPtr);+ __m128 accum2 = avx_hadd_R(accum);+ _mm_store_ss(outBuf + i, accum2);+ }+}++/*+ * Real coefficients, complex input+ */++void filterRC(int num, int numCoeffs, float *coeffs, float *inBuf, float *outBuf){+ int i, j;+ for(i=0; i<num*2; i+=2){+ float *startPtr = inBuf + i;+ dotprod_C(numCoeffs, coeffs, startPtr, outBuf + i);+ }+}++/*+ * SIMD versions+ */++void filterSSERC(int num, int numCoeffs, float *coeffs, float *inBuf, float *outBuf){+ int i;+ for(i=0; i<num*2; i+=2){+ float *startPtr = inBuf + i;+ __m128 accum = sse_dotprod_R(numCoeffs, coeffs, startPtr);+ accum = sse_hadd_C(accum);+ store_complex(outBuf + i, accum);+ }+}++void filterSSERC2(int num, int numCoeffs, float *coeffs, float *inBuf, float *outBuf){+ int i, j;+ for(i=0; i<num*2; i+=2){+ float *startPtr = inBuf + i;+ __m128 accum = sse_dotprod_C(numCoeffs, coeffs, startPtr);+ accum = sse_hadd_C(accum);+ store_complex(outBuf + i, accum);+ }+}++void filterAVXRC(int num, int numCoeffs, float *coeffs, float *inBuf, float *outBuf){+ int i;+ for(i=0; i<num*2; i+=2){+ float *startPtr = inBuf + i;+ __m256 accum = avx_dotprod_R(numCoeffs, coeffs, startPtr);+ __m128 accum2 = avx_hadd_C(accum);+ store_complex(outBuf + i, accum2);+ }+}++void filterAVXRC2(int num, int numCoeffs, float *coeffs, float *inBuf, float *outBuf){+ int i;+ for(i=0; i<num*2; i+=2){+ float *startPtr = inBuf + i;+ __m256 accum = avx_dotprod_C(numCoeffs, coeffs, startPtr);+ __m128 accum2 = avx_hadd_C(accum);+ store_complex(outBuf + i, accum2);+ }+}++/*+ * Symmetric versions+ */+void filterSSESymmetricRC(int num, int numCoeffs, float *coeffs, float *inBuf, float *outBuf){+ int i;+ for(i=0; i<num*2; i+=2){+ float *startPtr = inBuf + i;+ __m128 accum = sse_sym_dotprod_C(numCoeffs, coeffs, startPtr);+ accum = sse_hadd_C(accum);+ store_complex(outBuf + i, accum);+ }+}++void filterAVXSymmetricRC(int num, int numCoeffs, float *coeffs, float *inBuf, float *outBuf){+ int i;+ for(i=0; i<num*2; i+=2){+ float *startPtr = inBuf + i;+ __m256 accum = avx_sym_dotprod_C(numCoeffs, coeffs, startPtr);+ __m128 accum1 = avx_hadd_C(accum);+ store_complex(outBuf + i, accum1);+ }+}++/*+ * DC blocker+ */+void dcBlocker(int num, float lastSample, float lastOutput, float *finalSample, float *finalOutput, float *inBuf, float *outBuf){+ int i;+ for(i=0; i<num; i++){+ lastOutput = inBuf[i] - lastSample + 0.997 * lastOutput;+ outBuf[i] = lastOutput;+ lastSample = inBuf[i];+ }+ *finalSample = lastSample;+ *finalOutput = lastOutput;+}+
+ c_sources/resample.c view
@@ -0,0 +1,87 @@+/*+ * FIR rational resampling of complex and real data with real coefficients.+ * These exist because the pure Haskell implementations are slow.+ * Uses SIMD instructions for performance.+ */++#include <stdio.h>+#include <stdint.h>+#include <x86intrin.h>++#include "common.h"++/*+ * Rational downsampling+ */+void resample(int buf_size, int coeff_size, int interpolation, int decimation, int filter_offset, float *coeffs, float *in_buf, float *out_buf){+ int j, k, l;+ int input_offset = 0;+ for(k=0; k<buf_size; k++) {+ float accum = 0;++ for(l=0, j=filter_offset; j<coeff_size; l++, j+=interpolation) {+ accum += in_buf[input_offset + l] * coeffs[j];+ }++ int filter_offset_new = interpolation - 1 - (decimation - filter_offset - 1) % interpolation;+ input_offset += (decimation - filter_offset - 1) / interpolation + 1;+ filter_offset = filter_offset_new;++ out_buf[k] = accum;+ }+}++int resample2(int buf_size, int num_coeffs, int starting_group, int num_groups, int *increments, float **coeffs, float *in_buf, float *out_buf){+ int i, j;+ int group = starting_group;+ float *start_ptr = in_buf;++ for(i=0; i<buf_size; i++){+ out_buf[i] = dotprod_R(num_coeffs, coeffs[group], start_ptr);++ start_ptr += increments[group];+ group++;+ if(group == num_groups) + group = 0;+ }++ return group;+}++int resampleSSERR(int buf_size, int num_coeffs, int starting_group, int num_groups, int *increments, float **coeffs, float *in_buf, float *out_buf){+ int i, j;+ int group = starting_group;+ float *startPtr = in_buf;++ for(i=0; i<buf_size; i++){+ __m128 accum = sse_dotprod_R(num_coeffs, coeffs[group], startPtr);+ accum = sse_hadd_R(accum);+ _mm_store_ss(out_buf + i, accum);++ startPtr += increments[group];+ group++;+ if(group == num_groups) + group = 0;+ }++ return group;+}++int resampleAVXRR(int buf_size, int num_coeffs, int starting_group, int num_groups, int *increments, float **coeffs, float *in_buf, float *out_buf){+ int i, j;+ int group = starting_group;+ float *startPtr = in_buf;++ for(i=0; i<buf_size; i++){+ __m256 accum = avx_dotprod_R(num_coeffs, coeffs[group], startPtr);+ __m128 res = avx_hadd_R(accum);+ _mm_store_ss(out_buf + i, res);++ startPtr += increments[group];+ group++;+ if(group == num_groups) + group = 0;+ }++ return group;+}
+ c_sources/scale.c view
@@ -0,0 +1,38 @@+/*+ * Scaling of Vectors.+ * These exist because the pure Haskell implementations are slow.+ * Uses SIMD instructions for performance.+ */++#include <stdio.h>+#include <stdint.h>+#include <x86intrin.h>++/*+ * Scaling+ */++void scale(int num, float factor, float *in_buf, float *out_buf){+ int i;+ for(i=0; i<num; i++){+ out_buf[i] = in_buf[i] * factor;+ }+}++void scaleSSE(int num, float factor, float *in_buf, float *out_buf){+ int i;+ __m128 fac = _mm_set1_ps(factor);+ for(i=0; i<num; i+=4){+ _mm_storeu_ps(out_buf + i, _mm_mul_ps(fac, _mm_loadu_ps(in_buf + i)));+ }+}++void scaleAVX(int num, float factor, float *in_buf, float *out_buf){+ int i;+ __m256 fac = _mm256_set1_ps(factor);+ for(i=0; i<num; i+=8){+ _mm256_storeu_ps(out_buf + i, _mm256_mul_ps(fac, _mm256_loadu_ps(in_buf + i)));+ }+}++
+ hs_sources/SDR/ArgUtils.hs view
@@ -0,0 +1,19 @@+{-| Utilities for parsing command line arguments that might be useful when writing a SDR application. Uses the optparse-applicative library. -}+module SDR.ArgUtils (+ parseSize+ ) where++import Options.Applicative+import Data.Decimal++{-| Parse a number that may have a decimal point and a suffix, e.g. 2.56M -}+parseSize :: ReadM Integer+parseSize = eitherReader $ \arg -> case reads arg of+ [(r, suffix)] -> case suffix of + [] -> return $ round (r :: Decimal)+ "K" -> return $ round $ r * 1000 + "M" -> return $ round $ r * 1000000+ "G" -> return $ round $ r * 1000000000+ x -> Left $ "Cannot parse suffix: `" ++ x ++ "'"+ _ -> Left $ "Cannot parse value: `" ++ arg ++ "'"+
+ hs_sources/SDR/Demod.hs view
@@ -0,0 +1,43 @@+{-# LANGUAGE FlexibleContexts #-}++{-| FM demodulation pipes -}+module SDR.Demod (+ fmDemodStr,+ fmDemodVec,+ fmDemod+ ) where++import Data.Complex+import Data.Vector.Generic as VG+import Data.Vector.Fusion.Stream++import SDR.VectorUtils+import Pipes++-- | FM demodulate a stream of complex samples+{-# INLINE fmDemodStr #-}+fmDemodStr :: (RealFloat a) + => Complex a -- ^ The starting sample - i.e. the last sample in the last buffer+ -> Stream (Complex a) -- ^ The input stream+ -> Stream a -- ^ The output stream+fmDemodStr = mapAccumMV func + where+ {-# INLINE func #-}+ func last sample = return (sample, phase (sample * conjugate last))++-- | FM demodulate a vector of complex samples+{-# INLINE fmDemodVec #-}+fmDemodVec :: (RealFloat a, Vector v (Complex a), Vector v a) + => Complex a -- ^ The starting sample - i.e. the last sample in the last buffer+ -> v (Complex a) -- ^ The input Vector+ -> v a -- ^ The output Vector+fmDemodVec init = unstream . fmDemodStr init . stream++{-# INLINE fmDemod #-}+fmDemod :: (RealFloat a, Vector v (Complex a), Vector v a) => Pipe (v (Complex a)) (v a) IO ()+fmDemod = func 0+ where+ func lastSample = do+ dat <- await+ yield $ fmDemodVec lastSample dat+ func $ VG.unsafeIndex dat (VG.length dat - 1)
+ hs_sources/SDR/FFT.hs view
@@ -0,0 +1,162 @@+{-# LANGUAGE FlexibleContexts, ScopedTypeVariables #-}++{-| Fast FFTs using FFTW -}+module SDR.FFT (+ -- * Windows+ fftFixup,+ hamming,+ hanning,+ blackman,++ -- * FFTs+ fftw,+ fftwReal,+ fftwParallel+ ) where++import Control.Monad as CM+import Foreign.Storable+import Foreign.Storable.Complex+import Foreign.C.Types+import Data.Complex+import Foreign.ForeignPtr+import Control.Concurrent hiding (yield)+import qualified Data.Map as Map++import Data.Vector.Generic as VG+import Data.Vector.Generic.Mutable as VGM+import Data.Vector.Storable as VS+import Data.Vector.Storable.Mutable as VSM+import Data.Vector.Fusion.Stream as VFS++import Pipes+import Numeric.FFTW++import SDR.FilterDesign++mallocForeignBufferAligned :: forall a. Storable a => Int -> IO (ForeignPtr a)+mallocForeignBufferAligned elems = do+ ptr <- fftwMalloc $ fromIntegral $ elems * sizeOf (undefined :: a)+ newForeignPtr fftwFreePtr ptr++-- | Compute a vector of alternating 1s and 0s of the given size.+fftFixup :: (VG.Vector v n, Num n) + => Int -- ^ The length of the Vector+ -> v n+fftFixup size = VG.generate size func+ where+ func idx + | even idx = 1+ | otherwise = -1++-- | Creates a Pipe that performs a complex to complex DFT.+fftw :: (VG.Vector v (Complex CDouble)) + => Int -- ^ The size of the input and output buffers+ -> IO (Pipe (v (Complex CDouble)) (VS.Vector (Complex CDouble)) IO ())+fftw samples = do+ ina <- mallocForeignBufferAligned samples+ out <- mallocForeignBufferAligned samples++ plan <- withForeignPtr ina $ \ip -> + withForeignPtr out $ \op -> + planDFT1d samples ip op Forward fftwEstimate+ + return $ for cat $ \inv' -> do+ out <- lift $ mallocForeignBufferAligned samples+ ina <- lift $ mallocForeignBufferAligned samples+ let inv = VSM.unsafeFromForeignPtr0 ina samples++ lift $ VGM.fill inv $ VFS.mapM return $ VG.stream inv'++ let (fp, offset, length) = VSM.unsafeToForeignPtr inv++ lift $ withForeignPtr fp $ \fpp -> + withForeignPtr out $ \op -> + executeDFT plan fpp op++ yield $ VS.unsafeFromForeignPtr0 out samples++-- | Creates a pipe that performs a real to complex DFT.+fftwReal :: (VG.Vector v CDouble) + => Int -- ^ The size of the input Vector+ -> IO (Pipe (v CDouble) (VS.Vector (Complex CDouble)) IO ())+fftwReal samples = do+ --Allocate in and out buffers that wont be used because there doesnt seem to be a way to create a plan without them+ ina <- mallocForeignBufferAligned samples+ out <- mallocForeignBufferAligned samples++ plan <- withForeignPtr ina $ \ip -> + withForeignPtr out $ \op -> + planDFTR2C1d samples ip op fftwEstimate++ return $ for cat $ \inv' -> do+ out <- lift $ mallocForeignBufferAligned ((samples `quot` 2) + 1)+ ina <- lift $ mallocForeignBufferAligned samples+ let inv = VSM.unsafeFromForeignPtr0 ina samples++ lift $ VGM.fill inv $ VFS.mapM return $ VG.stream inv'+ let (fp, offset, length) = VSM.unsafeToForeignPtr inv++ lift $ withForeignPtr fp $ \fpp -> + withForeignPtr out $ \op -> + executeDFTR2C plan fpp op++ yield $ VS.unsafeFromForeignPtr0 out samples++{-| Creates a pipe that uses multiple threads to perform complex to complex DFTs in+ a pipelined fashion. Each time a buffer is consumed, it is given to+ a pool of threads to perform the DFT. Then, if a thread has finished+ performing a previous DFT, the result is yielded.+-}+fftwParallel :: (VG.Vector v (Complex CDouble)) + => Int -- ^ The number of threads to use+ -> Int -- ^ The size of the input Vector+ -> IO (Pipe (v (Complex CDouble)) (VS.Vector (Complex CDouble)) IO ())+fftwParallel threads samples = do+ --plan the DFT+ ina <- mallocForeignBufferAligned samples+ out <- mallocForeignBufferAligned samples++ plan <- withForeignPtr ina $ \ip -> + withForeignPtr out $ \op -> + planDFT1d samples ip op Forward fftwEstimate++ --setup the channels and worker threads+ inChan <- newChan + outMap <- newMVar Map.empty++ CM.replicateM threads $ forkIO $ forever $ do+ (idx, res) <- readChan inChan+ + out <- mallocForeignBufferAligned samples+ ina <- mallocForeignBufferAligned samples+ let inv = VSM.unsafeFromForeignPtr0 ina samples++ VGM.fill inv $ VFS.mapM return $ VG.stream res++ let (fp, offset, length) = VSM.unsafeToForeignPtr inv++ withForeignPtr fp $ \fpp -> + withForeignPtr out $ \op -> + executeDFT plan fpp op++ theMap <- takeMVar outMap+ putMVar outMap $ Map.insert idx (VS.unsafeFromForeignPtr0 out samples) theMap+ + --build the pipe+ let pipe nextIn nextOut = do+ dat <- await+ lift $ writeChan inChan (nextIn, dat)++ theMap <- lift $ takeMVar outMap+ case Map.lookup nextOut theMap of+ Nothing -> do+ lift $ putMVar outMap theMap+ pipe (nextIn + 1) nextOut+ Just dat -> do+ lift $ putMVar outMap $ Map.delete nextOut theMap+ yield dat+ pipe (nextIn + 1) (nextOut + 1)++ return $ pipe 0 0+
+ hs_sources/SDR/Filter.hs view
@@ -0,0 +1,496 @@+{-# LANGUAGE RecordWildCards, FlexibleContexts, GADTs #-}++{-| FIR filtering, decimation and resampling.++ FIR filters (and decimators, resamplers) work by taking successive dot products between the filter coefficients and the input data at increasing offsets. Sometimes the dot product fits entirely within one input buffer and other times it spans two input buffers (but never more because we assume that the filter length is less than the buffer size).++ We divide the filtering code by these two cases. Each filter (or decimator, resampler) is described by a data structure such as `Filter` with two functions, one for filtering within a single buffer and one that crosses buffers. ++ The user must first create one of these data structures using the helper functions and pass this data structure to one of `firFilter`, `firDecimator`, or `firResampler` to create the `Pipe` that does the filtering. For example:++ > decimatorStruct <- fastDecimatorC decimation coeffs+ > let decimatorPipe :: Pipe (Vector (Complex Float)) (Vector (Complex Float)) IO ()+ > decimatorPipe = firDecimator decimatorStruct outputSize++ There are polymorphic Haskell only implementations of filtering, decimation and resampling, for example, `haskellFilter`. In addition, there are optimised C implementations that use SIMD instructions on x86 machines, such as `fastFilterR`. These are always specialized to either real or complex numbers. There are also even faster implementations specialized for the case where the filter coefficients are symmetric as in a linear phase filter such as `fastSymmetricFilterR`.++ The Haskell implementations are reasonably fast due to the Vector library and GHC's LLVM backend, however, if speed is important you are much better off with the C implementations.++ In the future we may avoid the cross buffer filtering function by mapping the buffers consecutively in memory as (I believe) GNU Radio does.++ An extensive benchmark suite exists in the /benchmarks subdirectory of this package.+-}+module SDR.Filter (+ -- * Types+ Filter(..),+ Decimator(..),+ Resampler(..),++ -- * Helper Functions+ -- ** Filters+ haskellFilter,+ fastFilterR,+ fastFilterC,+ fastSymmetricFilterR,++ -- ** Decimators+ haskellDecimator,+ fastDecimatorR,+ fastDecimatorC,+ fastSymmetricDecimatorR,++ -- ** Resamplers+ haskellResampler,+ --fastResampler,++ -- * Filter+ firFilter,++ -- * Decimate+ firDecimator,++ -- * Resample+ firResampler,++ -- * DC Blocking Filter+ dcBlockingFilter+ ) where++import Data.Complex+import Control.Exception hiding (assert)+import qualified Data.Vector.Generic as VG+import qualified Data.Vector.Generic.Mutable as VGM+import qualified Data.Vector.Storable as VS+import Control.Monad.Primitive++import Pipes++import SDR.Util+import SDR.FilterInternal++{- | A `Filter` contains all of the information needed by the `filterr` + function to perform filtering. i.e. it contains the filter coefficients + and pointers to the functions to do the actual filtering.+-}+data Filter m v vm a = Filter {+ numCoeffsF :: Int,+ filterOne :: Int -> v a -> vm (PrimState m) a -> m (),+ filterCross :: Int -> v a -> v a -> vm (PrimState m) a -> m ()+}++{- | A `Decimator` contains all of the information needed by the `decimate`+ function to perform decimation i.e. it contains the filter coefficients + and pointers to the functions to do the actual decimation.+-}+data Decimator m v vm a = Decimator {+ numCoeffsD :: Int,+ decimationD :: Int,+ decimateOne :: Int -> v a -> vm (PrimState m) a -> m (),+ decimateCross :: Int -> v a -> v a -> vm (PrimState m) a -> m ()+}++{- | A `Resampler` contains all of the information needed by the `resample` + function to perform resampling i.e. it contains the filter coefficients + and pointers to the functions to do the actual resampling.+-}+data Resampler m v vm a = Resampler {+ numCoeffsR :: Int,+ decimationR :: Int,+ interpolationR :: Int,+ resampleOne :: Int -> Int -> v a -> vm (PrimState m) a -> m Int,+ resampleCross :: Int -> Int -> v a -> v a -> vm (PrimState m) a -> m Int+}++duplicate :: [a] -> [a]+duplicate = concat . map func + where func x = [x, x]++{-# INLINE haskellFilter #-}+-- | Returns a slow Filter data structure entirely implemented in Haskell+haskellFilter :: (PrimMonad m, Functor m, Num a, Mult a b, VG.Vector v a, VG.Vector v b, VGM.MVector vm a) + => [b] -- ^ The filter coefficients+ -> IO (Filter m v vm a) -- ^ The `Filter` data structure+haskellFilter coeffs = do+ let vCoeffs = VG.fromList coeffs+ evaluate vCoeffs+ let filterOne = filterHighLevel vCoeffs+ filterCross = filterCrossHighLevel vCoeffs+ numCoeffsF = length coeffs+ return $ Filter {..}++{-# INLINE fastFilterR #-}+-- | Returns a fast Filter data structure implemented in C using AVX instructions. For filtering real data with real coefficients.+fastFilterR :: [Float] -- ^ The filter coefficients+ -> IO (Filter IO VS.Vector VS.MVector Float) -- ^ The `Filter` data structure+fastFilterR coeffs = do+ let l = length coeffs+ ru = (l + 8 - 1) `quot` 8+ numCoeffsF = ru * 8 + diff = numCoeffsF - l+ vCoeffs = VG.fromList $ coeffs ++ replicate diff 0+ evaluate vCoeffs+ let filterOne = filterCAVXRR vCoeffs+ filterCross = filterCrossHighLevel vCoeffs+ return $ Filter {..}++{-# INLINE fastFilterC #-}+-- | Returns a fast Filter data structure implemented in C using AVX instructions. For filtering complex data with real coefficients.+fastFilterC :: [Float] -- ^ The filter coefficients+ -> IO (Filter IO VS.Vector VS.MVector (Complex Float)) -- ^ The `Filter` data structure+fastFilterC coeffs = do+ let l = length coeffs+ ru = (l + 4 - 1) `quot` 4+ numCoeffsF = ru * 4 + diff = numCoeffsF - l+ vCoeffs = VG.fromList $ duplicate $ coeffs ++ replicate diff 0+ vCoeffs2 = VG.fromList $ coeffs ++ replicate diff 0+ evaluate vCoeffs+ let filterOne = filterCAVXRC vCoeffs+ filterCross = filterCrossHighLevel vCoeffs2+ return $ Filter {..}++{-# INLINE fastSymmetricFilterR #-}+-- | Returns a fast Filter data structure implemented in C using AVX instructions. For filtering real data with real coefficients. For filters with symmetric coefficients, i.e. 'linear phase'. Coefficient length must be a multiple of 4.+fastSymmetricFilterR :: [Float] -- ^ The first half of the filter coefficients+ -> IO (Filter IO VS.Vector VS.MVector Float) -- ^ The `Filter` data structure+fastSymmetricFilterR coeffs = do+ let vCoeffs = VG.fromList coeffs + let vCoeffs2 = VG.fromList $ coeffs ++ reverse coeffs+ evaluate vCoeffs+ evaluate vCoeffs2+ let filterOne = filterCAVXSymmetricRR vCoeffs+ filterCross = filterCrossHighLevel vCoeffs2+ numCoeffsF = length coeffs * 2+ return $ Filter {..}++{-# INLINE haskellDecimator #-}+-- | Returns a slow Decimator data structure entirely implemented in Haskell+haskellDecimator :: (PrimMonad m, Functor m, Num a, Mult a b, VG.Vector v a, VG.Vector v b, VGM.MVector vm a) + => Int -- ^ The decimation factor+ -> [b] -- ^ The filter coefficients+ -> IO (Decimator m v vm a) -- ^ The `Decimator` data structure+haskellDecimator decimationD coeffs = do+ let vCoeffs = VG.fromList coeffs+ evaluate vCoeffs+ let decimateOne = decimateHighLevel decimationD vCoeffs+ decimateCross = decimateCrossHighLevel decimationD vCoeffs+ numCoeffsD = length coeffs+ return $ Decimator {..}++{-# INLINE fastDecimatorR #-}+-- | Returns a fast Decimator data structure implemented in C using AVX instructions. For decimating real data with real coefficients.+fastDecimatorR :: Int -- ^ The decimation factor+ -> [Float] -- ^ The filter coefficients+ -> IO (Decimator IO VS.Vector VS.MVector Float) -- ^ The `Decimator` data structure+fastDecimatorR decimationD coeffs = do+ let l = length coeffs+ ru = (l + 8 - 1) `quot` 8+ numCoeffsD = ru * 8 + diff = numCoeffsD - l+ vCoeffs = VG.fromList $ coeffs ++ replicate diff 0+ evaluate vCoeffs+ let decimateOne = decimateCAVXRR decimationD vCoeffs+ decimateCross = decimateCrossHighLevel decimationD vCoeffs+ return $ Decimator {..}++{-# INLINE fastDecimatorC #-}+-- | Returns a fast Decimator data structure implemented in C using AVX instructions. For decimating complex data with real coefficients.+fastDecimatorC :: Int -- ^ The decimation factor+ -> [Float] -- ^ The filter coefficients+ -> IO (Decimator IO VS.Vector VS.MVector (Complex Float)) -- ^ The `Decimator` data structure+fastDecimatorC decimationD coeffs = do+ let l = length coeffs+ ru = (l + 4 - 1) `quot` 4+ numCoeffsD = ru * 4 + diff = numCoeffsD - l+ vCoeffs = VG.fromList $ duplicate $ coeffs ++ replicate diff 0+ vCoeffs2 = VG.fromList $ coeffs ++ replicate diff 0+ evaluate vCoeffs+ let decimateOne = decimateCAVXRC decimationD vCoeffs+ decimateCross = decimateCrossHighLevel decimationD vCoeffs2+ return $ Decimator {..}++{-# INLINE fastSymmetricDecimatorR #-}+-- | Returns a fast Decimator data structure implemented in C using AVX instructions. For decimating real data with real coefficients. For decimators with symmetric coefficients, i.e. 'linear phase'. Coefficient length must be a multiple of 4.+fastSymmetricDecimatorR :: Int -- ^ The decimation factor+ -> [Float] -- ^ The first half of the filter coefficients+ -> IO (Decimator IO VS.Vector VS.MVector Float) -- ^ The `Decimator` data structure+fastSymmetricDecimatorR decimationD coeffs = do+ let vCoeffs = VG.fromList coeffs+ let vCoeffs2 = VG.fromList $ coeffs ++ reverse coeffs+ evaluate vCoeffs+ evaluate vCoeffs2+ let decimateOne = decimateCAVXSymmetricRR decimationD vCoeffs+ decimateCross = decimateCrossHighLevel decimationD vCoeffs2+ numCoeffsD = length coeffs * 2+ return $ Decimator {..}++{-# INLINE haskellResampler #-}+-- | Returns a slow Resampler data structure entirely implemented in Haskell+haskellResampler :: (PrimMonad m, Functor m, Num a, Mult a b, VG.Vector v a, VG.Vector v b, VGM.MVector vm a) + => Int -- ^ The interpolation factor+ -> Int -- ^ The decimation factor+ -> [b] -- ^ The filter coefficients+ -> IO (Resampler m v vm a) -- ^ The `Resampler` data structure+haskellResampler interpolationR decimationR coeffs = do+ let vCoeffs = VG.fromList coeffs+ evaluate vCoeffs+ let resampleOne = resampleHighLevel interpolationR decimationR vCoeffs+ resampleCross = resampleCrossHighLevel interpolationR decimationR vCoeffs+ numCoeffsR = length coeffs+ return $ Resampler {..}++{-+{-# INLINE fastResampler #-}+-- | Returns a fast Resampler data structure implemented in C using AVX instructions. For filtering real data with real coefficients.+fastResampler :: Int -- ^ The interpolation factor+ -> Int -- ^ The decimation factor+ -> [Float] -- ^ The filter coefficients+ -> IO (Resampler IO VS.Vector VS.MVector Float) -- ^ The `Resampler` data structure+fastResampler interpolationR decimationR coeffs = do+ let vCoeffs = VG.fromList coeffs+ evaluate vCoeffs+ resamp <- resampleCAVXRR interpolationR decimationR coeffs+ let resampleOne = resamp+ resampleCross = resampleCrossHighLevel interpolationR decimationR vCoeffs+ numCoeffsR = length coeffs+ return $ Resampler {..}+-}++data Buffer v a = Buffer {+ buffer :: v a,+ offset :: Int+}++space Buffer{..} = VGM.length buffer - offset++newBuffer :: (PrimMonad m, VGM.MVector vm a) => Int -> m (Buffer (vm (PrimState m)) a)+newBuffer size = do+ buf <- VGM.new size+ return $ Buffer buf 0++advanceOutBuf :: (PrimMonad m, VG.Vector v a) => Int -> Buffer (VG.Mutable v (PrimState m)) a -> Int -> Pipe b (v a) m (Buffer (VG.Mutable v (PrimState m)) a)+advanceOutBuf blockSizeOut buf@(Buffer bufOut offsetOut) count = + if count == space buf then do+ bufOutF <- lift $ VG.unsafeFreeze bufOut+ yield bufOutF+ lift $ newBuffer blockSizeOut+ else + return $ Buffer bufOut (offsetOut + count) ++-- | My own assert implementation since the GHC one doesnt seem to work even with optimisations disabled and using -fno-ignore-asserts+assert loc False = error loc+assert loc True = return ()++--Filtering+{-# INLINE firFilter #-}+{-| Create a pipe that performs filtering -}+firFilter :: (PrimMonad m, Functor m, VG.Vector v a, Num a) + => Filter m v (VG.Mutable v) a -- ^ The `Filter` data structure+ -> Int -- ^ The output block size+ -> Pipe (v a) (v a) m () -- ^ The `Pipe` that does the filtering+firFilter Filter{..} blockSizeOut = do+ inBuf <- await+ outBuf <- lift $ newBuffer blockSizeOut+ simple inBuf outBuf ++ where++ simple bufIn bufferOut@(Buffer bufOut offsetOut) = do+ assert "filter 1" (VG.length bufIn >= numCoeffsF)++ let count = min (VG.length bufIn - numCoeffsF + 1) (space bufferOut)+ lift $ filterOne count bufIn (VGM.unsafeDrop offsetOut bufOut)++ bufferOut' <- advanceOutBuf blockSizeOut bufferOut count+ let bufIn' = VG.drop count bufIn++ case VG.length bufIn' < numCoeffsF of+ False -> simple bufIn' bufferOut'+ True -> do+ next <- await+ crossover bufIn' next bufferOut'++ crossover bufLast bufNext bufferOut@(Buffer bufOut offsetOut) = do+ assert "filter 2" (VG.length bufLast < numCoeffsF) + assert "filter 3" (VG.length bufLast > 0) ++ let count = min (VG.length bufLast) (space bufferOut)+ lift $ filterCross count bufLast bufNext (VGM.unsafeDrop offsetOut bufOut)++ bufferOut' <- advanceOutBuf blockSizeOut bufferOut count++ case VG.length bufLast == count of + True -> simple bufNext bufferOut'+ False -> crossover (VG.drop count bufLast) bufNext bufferOut'++--Decimation+{-# INLINE firDecimator #-}+{-| Create a pipe that performs decimation -}+firDecimator :: (PrimMonad m, Functor m, VG.Vector v a, Num a) + => Decimator m v (VG.Mutable v) a -- ^ The `Decimator` data structure+ -> Int -- ^ The output block size+ -> Pipe (v a) (v a) m () -- ^ The `Pipe` that does the decimation+firDecimator Decimator{..} blockSizeOut = do+ inBuf <- await+ outBuf <- lift $ newBuffer blockSizeOut+ simple inBuf outBuf++ where++ simple bufIn bufferOut@(Buffer bufOut offsetOut) = do+ assert "decimate 1" (VG.length bufIn >= numCoeffsD)++ let count = min (((VG.length bufIn - numCoeffsD) `quot` decimationD) + 1) (space bufferOut)+ lift $ decimateOne count bufIn (VGM.unsafeDrop offsetOut bufOut)++ bufferOut' <- advanceOutBuf blockSizeOut bufferOut count+ let bufIn' = VG.drop (count * decimationD) bufIn++ case VG.length bufIn' < numCoeffsD of+ False -> simple bufIn' bufferOut'+ True -> do+ next <- await+ crossover bufIn' next bufferOut'++ crossover bufLast bufNext bufferOut@(Buffer bufOut offsetOut) = do+ assert "decimate 2" (VG.length bufLast < numCoeffsD) + assert "decimate 3" (VG.length bufLast > 0) ++ let count = min (VG.length bufLast `quotUp` decimationD) (space bufferOut)+ lift $ decimateCross count bufLast bufNext (VGM.unsafeDrop offsetOut bufOut)++ bufferOut' <- advanceOutBuf blockSizeOut bufferOut count++ case VG.length bufLast <= count * decimationD of + True -> simple (VG.drop (count * decimationD - VG.length bufLast) bufNext) bufferOut'+ False -> crossover (VG.drop (count * decimationD) bufLast) bufNext bufferOut'++{-+Rational Downsampling:++Input upsampled by 3: |**|**|**|**|**|**|**|**|**|**|**|+Output downsampled by 7: |******|******|******|******|*****++ Consider here ^+ Next output is here ^++Filter offset is 2++k is number of used inputs++filterOffset + k*interpolation = decimation + filterOffset'+where+ k > 0+ 0 <= filterOffset, filterOffset' < interpolation++k*interpolation - filterOffset' = decimation - filterOffset+k*interpolation - filterOffset' - 1 = decimation - filterOffset - 1++(k-1) * interpolation + (interpolation - filterOffset' - 1) = decimation - filterOffset - 1++k = (decimation - filterOffset - 1) / interpolation + 1+filterOffset' = interpolation - 1 - (decimation - filterOffset - 1) % interpolation++Only works if decimation > interpolation++-}++{-+Rational Upsampling:++Input upsampled by 7: |******|******|******|******|*****+Output downsampled by 3: |**|**|**|**|**|**|**|**|**|**|**|++ Consider Here ^ + Next sample is ^++Filter offset is 4++filterOffset + k * interpolation = decimation + filterOffset'+where+ k = {0, 1}+ 0 <= filterOffset, filterOffset' < interpolation++k * interpolation + (interpolation - filterOffset' - 1) = decimation - filterOffset + interpolation - 1++k = (decimation - filterOffset + interpolation - 1) / interpolation++============================++Or, equivalently, ++k = 0 | filterOffset >= decimation+ 1 | otherwise++o = o - decimation + k * interpolation++-}++--Rational resampling+quotUp q d = (q + (d - 1)) `quot` d++{-# INLINE firResampler #-}+{-| Create a pipe that performs resampling -}+firResampler :: (PrimMonad m, VG.Vector v a, Num a) + => Resampler m v (VG.Mutable v) a -- ^ The `Resampler` data structure+ -> Int -- ^ The output block size+ -> Pipe (v a) (v a) m () -- ^ The `Pipe` that does the resampling+firResampler Resampler{..} blockSizeOut = do+ inBuf <- await+ outBuf <- lift $ newBuffer blockSizeOut+ simple inBuf outBuf 0++ where++ simple bufIn bufferOut@(Buffer bufOut offsetOut) filterOffset = do+ assert "resample 1" (VG.length bufIn * interpolationR >= numCoeffsR - filterOffset)+ --available number of samples == interpolation * num_input+ --required number of samples == decimation * (num_output - 1) + filter_length - filter_offset+ let count = min (((VG.length bufIn * interpolationR - numCoeffsR + filterOffset) `quot` decimationR) + 1) (space bufferOut)+ endOffset <- lift $ resampleOne filterOffset count bufIn (VGM.unsafeDrop offsetOut bufOut)+ assert "resample 2" ((count * decimationR + endOffset - filterOffset) `rem` interpolationR == 0)+ bufferOut' <- advanceOutBuf blockSizeOut bufferOut count+ --samples no longer needed starting from filterOffset == count * decimation - filterOffset+ --inputs lying in this region == (count * decimation - filterOffset) / interpolation (rounding up)+ let usedInput = (count * decimationR - filterOffset) `quotUp` interpolationR + bufIn' = VG.drop usedInput bufIn++ case VG.length bufIn' * interpolationR < numCoeffsR - endOffset of+ False -> simple bufIn' bufferOut' endOffset+ True -> do+ next <- await+ --TODO: why is this not needed in filter and decimator+ case VG.length bufIn' == 0 of+ True -> simple next bufferOut' endOffset+ False -> crossover bufIn' next bufferOut' endOffset++ crossover bufLast bufNext bufferOut@(Buffer bufOut offsetOut) filterOffset = do+ assert "resample 3" (VG.length bufLast * interpolationR < numCoeffsR - filterOffset)+ --outputsComputable is the number of outputs that need to be computed for the last buffer to no longer be needed+ --outputsComputable * decimation == numInput * interpolation + filterOffset + k+ let outputsComputable = (VG.length bufLast * interpolationR + filterOffset) `quotUp` decimationR+ count = min outputsComputable (space bufferOut)+ assert "resample 4" (count /= 0)+ endOffset <- lift $ resampleCross filterOffset count bufLast bufNext (VGM.unsafeDrop offsetOut bufOut)+ assert "resample 5" ((count * decimationR + endOffset - filterOffset) `rem` interpolationR == 0)+ bufferOut' <- advanceOutBuf blockSizeOut bufferOut count++ let inputUsed = (count * decimationR - filterOffset) `quotUp` interpolationR++ case inputUsed >= VG.length bufLast of + True -> simple (VG.drop (inputUsed - VG.length bufLast) bufNext) bufferOut' endOffset+ False -> crossover (VG.drop inputUsed bufLast) bufNext bufferOut' endOffset++-- | A DC blocking filter+dcBlockingFilter :: Pipe (VS.Vector Float) (VS.Vector Float) IO ()+dcBlockingFilter = func 0 0 + where+ func lastSample lastOutput = do+ dat <- await+ out <- lift $ VGM.new (VG.length dat)+ (lastSample, lastOutput) <- lift $ dcBlocker (VG.length dat) lastSample lastOutput dat out+ outF <- lift $ VG.unsafeFreeze out+ yield outF+ func lastSample lastOutput+
+ hs_sources/SDR/FilterDesign.hs view
@@ -0,0 +1,76 @@+{-| Filter design and plotting of frequency responses. -}++module SDR.FilterDesign (+ -- * Sinc Function+ sinc,++ -- * Windows+ hanning,+ hamming,+ blackman,++ -- * Convenience Functions+ windowedSinc,+ + -- * Frequency Response Plot+ plotFrequency+ ) where++import Graphics.Rendering.Chart.Easy+import Graphics.Rendering.Chart.Backend.Cairo+import Data.Complex++import qualified Data.Vector.Generic as VG++-- | Compute a sinc function+sinc :: (Floating n, VG.Vector v n)+ => Int -- ^ The length. Must be odd.+ -> n -- ^ The cutoff frequency (from 0 to 1)+ -> v n+sinc size cutoff = VG.generate size (func . (-) ((size - 1) `quot` 2))+ where+ func 0 = cutoff+ func idx = sin (pi * cutoff * fromIntegral idx) / (fromIntegral idx * pi)++-- | Compute a Hanning window.+hanning :: (Floating n, VG.Vector v n) + => Int -- ^ The length of the window+ -> v n+hanning size = VG.generate size func+ where+ func idx = 0.5 * (1 - cos((2 * pi * fromIntegral idx) / (fromIntegral size - 1)))+ +-- | Compute a Hamming window. +hamming :: (Floating n, VG.Vector v n) + => Int -- ^ The length of the window+ -> v n+hamming size = VG.generate size func+ where+ func idx = 0.54 - 0.46 * cos((2 * pi * fromIntegral idx) / (fromIntegral size - 1))+ +-- | Compute a Blackman window.+blackman :: (Floating n, VG.Vector v n) + => Int -- ^ The length of the window+ -> v n+blackman size = VG.generate size func+ where+ func idx = 0.42 - 0.5 * cos((2 * pi * fromIntegral idx) / (fromIntegral size - 1)) + 0.08 * cos((4 * pi * fromIntegral idx) / (fromIntegral size - 1))++windowedSinc :: (Floating n, VG.Vector v n)+ => Int -- ^ The length+ -> n -- ^ The cutoff frequency (from 0 to 1)+ -> (Int -> v n) -- ^ The window function+ -> v n+windowedSinc size cutoff window = VG.zipWith (*) (sinc size cutoff) (window size)++signal :: [Double] -> [Double] -> [(Double, Double)]+signal coeffs xs = [ (x / pi, func x) | x <- xs ]+ where+ func phase = magnitude $ sum $ zipWith (\index mag -> mkPolar mag (phase * (- index))) (iterate (+ 1) (- ((fromIntegral (length coeffs) - 1) / 2))) coeffs++-- | Given filter coefficients, plot their frequency response and save the graph as "frequency_response.png".+plotFrequency :: [Double] -- ^ The filter coefficients+ -> IO ()+plotFrequency coeffs = toFile def "frequency_response.png" $ do+ layout_title .= "Frequency Response"+ plot (line "Frequency Response" [signal coeffs $ takeWhile (< pi) $ iterate (+ 0.01) 0])
+ hs_sources/SDR/FilterInternal.hs view
@@ -0,0 +1,401 @@+{-# LANGUAGE ScopedTypeVariables, BangPatterns, RecordWildCards #-}++{-| Functions used internally by the SDR.Filter module. Most of these are+ not actually used but exist for benchmarking purposes to determine the + fastest filter implementation.+-}+module SDR.FilterInternal where++import Control.Monad.Primitive +import Control.Monad+import Foreign.C.Types+import Foreign.Ptr+import Unsafe.Coerce+import Data.Complex+import Foreign.Marshal.Array+import Foreign.Marshal.Alloc+import Foreign.Storable++import qualified Data.Vector.Generic as VG+import qualified Data.Vector.Generic.Mutable as VGM+import qualified Data.Vector.Storable as VS+import qualified Data.Vector.Storable.Mutable as VSM+import qualified Data.Vector.Fusion.Stream as VFS+import qualified Data.Vector.Fusion.Stream.Monadic as VFSM++import SDR.VectorUtils+import SDR.Util++{-# INLINE filterHighLevel #-}+filterHighLevel :: (PrimMonad m, Functor m, Num a, Mult a b, VG.Vector v a, VG.Vector v b, VGM.MVector vm a) => v b -> Int -> v a -> vm (PrimState m) a -> m ()+filterHighLevel coeffs num inBuf outBuf = fill (VFSM.generate num dotProd) outBuf+ where+ dotProd offset = VG.sum $ VG.zipWith mult (VG.unsafeDrop offset inBuf) coeffs++{-# INLINE filterImperative1 #-}+filterImperative1 :: (PrimMonad m, Functor m, Num a, Mult a b, VG.Vector v a, VG.Vector v b, VGM.MVector vm a) => v b -> Int -> v a -> vm (PrimState m) a -> m ()+filterImperative1 coeffs num inBuf outBuf = go 0+ where+ go offset + | offset < num = do+ let res = dotProd offset+ VGM.unsafeWrite outBuf offset res+ go $ offset + 1+ | otherwise = return ()+ dotProd offset = VG.sum $ VG.zipWith mult (VG.unsafeDrop offset inBuf) coeffs++{-# INLINE filterImperative2 #-}+filterImperative2 :: (PrimMonad m, Functor m, Num a, Mult a b, VG.Vector v a, VG.Vector v b, VGM.MVector vm a) => v b -> Int -> v a -> vm (PrimState m) a -> m ()+filterImperative2 coeffs num inBuf outBuf = go 0+ where+ go offset + | offset < num = do+ let res = dotProd (VG.unsafeDrop offset inBuf)+ VGM.unsafeWrite outBuf offset res+ go $ offset + 1+ | otherwise = return ()+ dotProd buf = go 0 0+ where+ go !accum j + | j < VG.length coeffs = go (VG.unsafeIndex buf j `mult` VG.unsafeIndex coeffs j + accum) (j + 1)+ | otherwise = accum++type FilterCRR = CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()+type FilterRR = VS.Vector Float -> Int -> VS.Vector Float -> VS.MVector RealWorld Float -> IO ()+type FilterRC = VS.Vector Float -> Int -> VS.Vector (Complex Float) -> VS.MVector RealWorld (Complex Float) -> IO ()++filterFFIR :: FilterCRR -> FilterRR +filterFFIR func coeffs num inBuf outBuf = + VS.unsafeWith (unsafeCoerce coeffs) $ \cPtr -> + VS.unsafeWith (unsafeCoerce inBuf) $ \iPtr -> + VSM.unsafeWith (unsafeCoerce outBuf) $ \oPtr -> + func (fromIntegral num) (fromIntegral $ VG.length coeffs) cPtr iPtr oPtr++filterFFIC :: FilterCRR -> FilterRC +filterFFIC func coeffs num inBuf outBuf = + VS.unsafeWith (unsafeCoerce coeffs) $ \cPtr -> + VS.unsafeWith (unsafeCoerce inBuf) $ \iPtr -> + VSM.unsafeWith (unsafeCoerce outBuf) $ \oPtr -> + func (fromIntegral num) (fromIntegral $ VG.length coeffs) cPtr iPtr oPtr++foreign import ccall unsafe "filterRR"+ filterRR_c :: CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()++filterCRR :: FilterRR+filterCRR = filterFFIR filterRR_c ++foreign import ccall unsafe "filterRC"+ filterRC_c :: CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()++filterCRC :: FilterRC+filterCRC = filterFFIC filterRC_c++foreign import ccall unsafe "filterSSERR"+ filterSSERR_c :: CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()++filterCSSERR :: FilterRR+filterCSSERR = filterFFIR filterSSERR_c++foreign import ccall unsafe "filterSSESymmetricRR"+ filterSSESymmetricRR_c :: CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()++filterCSSESymmetricRR :: FilterRR+filterCSSESymmetricRR = filterFFIR filterSSESymmetricRR_c++foreign import ccall unsafe "filterSSERC"+ filterSSERC_c :: CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()++filterCSSERC :: FilterRC+filterCSSERC = filterFFIC filterSSERC_c++foreign import ccall unsafe "filterSSERC2"+ filterSSERC2_c :: CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()++filterCSSERC2 :: FilterRC+filterCSSERC2 = filterFFIC filterSSERC2_c++foreign import ccall unsafe "filterAVXRR"+ filterAVXRR_c :: CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()++filterCAVXRR :: FilterRR+filterCAVXRR = filterFFIR filterAVXRR_c++foreign import ccall unsafe "filterAVXSymmetricRR"+ filterAVXSymmetricRR_c :: CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()++filterCAVXSymmetricRR :: FilterRR+filterCAVXSymmetricRR = filterFFIR filterAVXSymmetricRR_c++foreign import ccall unsafe "filterAVXRC"+ filterAVXRC_c :: CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()++filterCAVXRC :: FilterRC+filterCAVXRC = filterFFIC filterAVXRC_c++foreign import ccall unsafe "filterAVXRC2"+ filterAVXRC2_c :: CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()++filterCAVXRC2 :: FilterRC+filterCAVXRC2 = filterFFIC filterAVXRC2_c++foreign import ccall unsafe "filterSSESymmetricRC"+ filterSSESymmetricRC_c :: CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()++filterCSSESymmetricRC :: FilterRC+filterCSSESymmetricRC = filterFFIC filterSSESymmetricRC_c++foreign import ccall unsafe "filterAVXSymmetricRC"+ filterAVXSymmetricRC_c :: CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()++filterCAVXSymmetricRC :: FilterRC+filterCAVXSymmetricRC = filterFFIC filterAVXSymmetricRC_c++-- Decimation++{-# INLINE decimateHighLevel #-}+decimateHighLevel :: (PrimMonad m, Functor m, Num a, Mult a b, VG.Vector v a, VG.Vector v b, VGM.MVector vm a) => Int -> v b -> Int -> v a -> vm (PrimState m) a -> m ()+decimateHighLevel factor coeffs num inBuf outBuf = fill x outBuf+ where + x = VFSM.map dotProd (VFSM.iterateN num (+ factor) 0)+ dotProd offset = VG.sum $ VG.zipWith mult (VG.unsafeDrop offset inBuf) coeffs++type DecimateCRR = CInt -> CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()+type DecimateRR = Int -> VS.Vector Float -> Int -> VS.Vector Float -> VS.MVector RealWorld Float -> IO ()+type DecimateRC = Int -> VS.Vector Float -> Int -> VS.Vector (Complex Float) -> VS.MVector RealWorld (Complex Float) -> IO ()++decimateFFIR :: DecimateCRR -> DecimateRR +decimateFFIR func factor coeffs num inBuf outBuf = + VS.unsafeWith (unsafeCoerce coeffs) $ \cPtr -> + VS.unsafeWith (unsafeCoerce inBuf) $ \iPtr -> + VSM.unsafeWith (unsafeCoerce outBuf) $ \oPtr -> + func (fromIntegral num) (fromIntegral factor) (fromIntegral $ VG.length coeffs) cPtr iPtr oPtr++decimateFFIC :: DecimateCRR -> DecimateRC +decimateFFIC func factor coeffs num inBuf outBuf = + VS.unsafeWith (unsafeCoerce coeffs) $ \cPtr -> + VS.unsafeWith (unsafeCoerce inBuf) $ \iPtr -> + VSM.unsafeWith (unsafeCoerce outBuf) $ \oPtr -> + func (fromIntegral num) (fromIntegral factor) (fromIntegral $ VG.length coeffs) cPtr iPtr oPtr++foreign import ccall unsafe "decimateRR"+ decimateCRR_c :: CInt -> CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()++decimateCRR :: DecimateRR+decimateCRR = decimateFFIR decimateCRR_c++foreign import ccall unsafe "decimateRC"+ decimateCRC_c :: CInt -> CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()++decimateCRC :: DecimateRC+decimateCRC = decimateFFIC decimateCRC_c++foreign import ccall unsafe "decimateSSERR"+ decimateSSERR_c :: CInt -> CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()++decimateCSSERR :: DecimateRR+decimateCSSERR = decimateFFIR decimateSSERR_c++foreign import ccall unsafe "decimateSSERC"+ decimateSSERC_c :: CInt -> CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()++decimateCSSERC :: DecimateRC+decimateCSSERC = decimateFFIC decimateSSERC_c++foreign import ccall unsafe "decimateSSERC2"+ decimateSSERC2_c :: CInt -> CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()++decimateCSSERC2 :: DecimateRC+decimateCSSERC2 = decimateFFIC decimateSSERC2_c++foreign import ccall unsafe "decimateAVXRR"+ decimateAVXRR_c :: CInt -> CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()++decimateCAVXRR :: DecimateRR+decimateCAVXRR = decimateFFIR decimateAVXRR_c++foreign import ccall unsafe "decimateAVXRC"+ decimateAVXRC_c :: CInt -> CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()++decimateCAVXRC :: DecimateRC+decimateCAVXRC = decimateFFIC decimateAVXRC_c++foreign import ccall unsafe "decimateAVXRC2"+ decimateAVXRC2_c :: CInt -> CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()++decimateCAVXRC2 :: DecimateRC+decimateCAVXRC2 = decimateFFIC decimateAVXRC2_c++foreign import ccall unsafe "decimateSSESymmetricRR"+ decimateSSESymmetricRR_c :: CInt -> CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()++decimateCSSESymmetricRR :: DecimateRR+decimateCSSESymmetricRR = decimateFFIR decimateSSESymmetricRR_c++foreign import ccall unsafe "decimateAVXSymmetricRR"+ decimateAVXSymmetricRR_c :: CInt -> CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()++decimateCAVXSymmetricRR :: DecimateRR+decimateCAVXSymmetricRR = decimateFFIR decimateAVXSymmetricRR_c++foreign import ccall unsafe "decimateSSESymmetricRC"+ decimateSSESymmetricRC_c :: CInt -> CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()++decimateCSSESymmetricRC :: DecimateRC+decimateCSSESymmetricRC = decimateFFIC decimateSSESymmetricRC_c++foreign import ccall unsafe "decimateAVXSymmetricRC"+ decimateAVXSymmetricRC_c :: CInt -> CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()++decimateCAVXSymmetricRC :: DecimateRC+decimateCAVXSymmetricRC = decimateFFIC decimateAVXSymmetricRC_c++-- Resampling+{-# INLINE resampleHighLevel #-}+resampleHighLevel :: (PrimMonad m, Num a, Mult a b, VG.Vector v a, VG.Vector v b, VGM.MVector vm a) => Int -> Int -> v b -> Int -> Int -> v a -> vm (PrimState m) a -> m Int+resampleHighLevel interpolation decimation coeffs filterOffset count inBuf outBuf = fill 0 filterOffset 0+ where+ fill i filterOffset inputOffset+ | i < count = do+ let dp = dotProd filterOffset inputOffset+ VGM.unsafeWrite outBuf i dp+ let (q, r) = divMod (decimation - filterOffset - 1) interpolation+ inputOffset' = inputOffset + q + 1+ filterOffset' = interpolation - 1 - r+ filterOffset' `seq` inputOffset' `seq` fill (i + 1) filterOffset' inputOffset'+ | otherwise = return filterOffset+ dotProd filterOffset offset = VG.sum $ VG.zipWith mult (VG.unsafeDrop offset inBuf) (stride interpolation (VG.unsafeDrop filterOffset coeffs))++foreign import ccall unsafe "resample"+ resample_c :: CInt -> CInt -> CInt -> CInt -> CInt -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()++resampleCRR :: Int -> Int -> Int -> Int -> VS.Vector Float -> VS.Vector Float -> VS.MVector RealWorld Float -> IO ()+resampleCRR num interpolation decimation offset coeffs inBuf outBuf = + VS.unsafeWith (unsafeCoerce coeffs) $ \cPtr -> + VS.unsafeWith (unsafeCoerce inBuf) $ \iPtr -> + VS.unsafeWith (unsafeCoerce outBuf) $ \oPtr -> + resample_c (fromIntegral num) (fromIntegral $ VG.length coeffs) (fromIntegral interpolation) (fromIntegral decimation) (fromIntegral offset) cPtr iPtr oPtr++pad :: a -> Int -> [a] -> [a]+pad with num list = list ++ replicate (num - length list) with ++strideList :: Int -> [a] -> [a]+strideList s xs = go 0 xs+ where+ go _ [] = []+ go 0 (x:xs) = x : go (s-1) xs+ go n (x:xs) = go (n - 1) xs++roundUp :: Int -> Int -> Int+roundUp num div = ((num + div - 1) `quot` div) * div++data Coeffs = Coeffs {+ numCoeffs :: Int,+ numGroups :: Int,+ increments :: [Int],+ groups :: [[Float]]+}++prepareCoeffs :: Int -> Int -> Int -> [Float] -> Coeffs+prepareCoeffs n interpolation decimation coeffs = Coeffs {..}+ where+ numCoeffs = maximum $ map (length . snd) dats+ numGroups = length groups+ increments = map fst dats++ groups :: [[Float]]+ groups = map (pad 0 (roundUp numCoeffs n)) $ map snd dats++ dats :: [(Int, [Float])]+ dats = func 0+ where++ func' 0 = []+ func' x = func x++ func :: Int -> [(Int, [Float])]+ func offset = (increment, strideList interpolation $ drop offset coeffs) : func' offset'+ where+ (q, r) = divMod (decimation - offset - 1) interpolation+ increment = q + 1+ offset' = interpolation - 1 - r++resampleFFIR :: (Ptr CFloat -> Ptr CFloat -> IO CInt) -> VS.Vector Float -> VSM.MVector RealWorld Float -> IO Int+resampleFFIR func inBuf outBuf = liftM fromIntegral $+ VS.unsafeWith (unsafeCoerce inBuf) $ \iPtr -> + VS.unsafeWith (unsafeCoerce outBuf) $ \oPtr -> + func iPtr oPtr++type ResampleR = CInt -> CInt -> CInt -> CInt -> Ptr CInt -> Ptr (Ptr CFloat) -> Ptr CFloat -> Ptr CFloat -> IO CInt++mkResampler :: ResampleR -> Int -> Int -> Int -> [Float] -> IO (Int -> Int -> VS.Vector Float -> VS.MVector RealWorld Float -> IO Int)+mkResampler func n interpolation decimation coeffs = do+ groupsP <- mapM newArray $ map (map realToFrac) groups+ groupsPP <- newArray groupsP+ incrementsP <- newArray $ map fromIntegral increments+ return $ \num offset -> resampleFFIR $ func (fromIntegral num) (fromIntegral numCoeffs) (fromIntegral offset) (fromIntegral numGroups) incrementsP groupsPP+ where+ Coeffs {..} = prepareCoeffs n interpolation decimation coeffs++foreign import ccall unsafe "resample2"+ resample2_c :: CInt -> CInt -> CInt -> CInt -> Ptr CInt -> Ptr (Ptr CFloat) -> Ptr CFloat -> Ptr CFloat -> IO CInt++resampleCRR2 :: Int -> Int -> [Float] -> IO (Int -> Int -> VS.Vector Float -> VS.MVector RealWorld Float -> IO Int)+resampleCRR2 = mkResampler resample2_c 1++foreign import ccall unsafe "resampleSSERR"+ resampleCSSERR_c :: CInt -> CInt -> CInt -> CInt -> Ptr CInt -> Ptr (Ptr CFloat) -> Ptr CFloat -> Ptr CFloat -> IO CInt++resampleCSSERR :: Int -> Int -> [Float] -> IO (Int -> Int -> VS.Vector Float -> VS.MVector RealWorld Float -> IO Int)+resampleCSSERR = mkResampler resampleCSSERR_c 4++foreign import ccall unsafe "resampleAVXRR"+ resampleAVXRR_c :: CInt -> CInt -> CInt -> CInt -> Ptr CInt -> Ptr (Ptr CFloat) -> Ptr CFloat -> Ptr CFloat -> IO CInt++resampleCAVXRR :: Int -> Int -> [Float] -> IO (Int -> Int -> VS.Vector Float -> VS.MVector RealWorld Float -> IO Int)+resampleCAVXRR = mkResampler resampleAVXRR_c 8++{-+ - Cross buffer+-}++{-# INLINE decimateCrossHighLevel #-}+decimateCrossHighLevel :: (PrimMonad m, Functor m, Num a, Mult a b, VG.Vector v a, VG.Vector v b, VGM.MVector vm a) => Int -> v b -> Int -> v a -> v a -> vm (PrimState m) a -> m ()+decimateCrossHighLevel factor coeffs num lastBuf nextBuf outBuf = fill x outBuf+ where+ x = VFSM.map dotProd (VFSM.iterateN num (+ factor) 0)+ dotProd i = VG.sum $ VG.zipWith mult (VG.unsafeDrop i lastBuf VG.++ nextBuf) coeffs++{-# INLINE filterCrossHighLevel #-}+filterCrossHighLevel :: (PrimMonad m, Functor m, Num a, Mult a b, VG.Vector v a, VG.Vector v b, VGM.MVector vm a) => v b -> Int -> v a -> v a -> vm (PrimState m) a -> m ()+filterCrossHighLevel coeffs num lastBuf nextBuf outBuf = fill (VFSM.generate num dotProd) outBuf+ where+ dotProd i = VG.sum $ VG.zipWith mult (VG.unsafeDrop i lastBuf VG.++ nextBuf) coeffs++{-# INLINE resampleCrossHighLevel #-}+resampleCrossHighLevel :: (PrimMonad m, Num a, Mult a b, VG.Vector v a, VG.Vector v b, VGM.MVector vm a) => Int -> Int -> v b -> Int -> Int -> v a -> v a -> vm (PrimState m) a -> m Int+resampleCrossHighLevel interpolation decimation coeffs filterOffset count lastBuf nextBuf outBuf = fill 0 filterOffset 0+ where+ fill i filterOffset inputOffset+ | i < count = do+ let dp = dotProd filterOffset inputOffset+ VGM.unsafeWrite outBuf i dp+ let (q, r) = divMod (decimation - filterOffset - 1) interpolation+ inputOffset' = inputOffset + q + 1+ filterOffset' = interpolation - 1 - r+ filterOffset' `seq` inputOffset' `seq` fill (i + 1) filterOffset' inputOffset'+ | otherwise = return filterOffset+ dotProd filterOffset i = VG.sum $ VG.zipWith mult (VG.unsafeDrop i lastBuf VG.++ nextBuf) (stride interpolation (VG.unsafeDrop filterOffset coeffs))++foreign import ccall unsafe "dcBlocker"+ c_dcBlocker :: CInt -> CFloat -> CFloat -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()++dcBlocker :: Int -> Float -> Float -> VS.Vector Float -> VS.MVector RealWorld Float -> IO (Float, Float)+dcBlocker num lastSample lastOutput inBuf outBuf = + alloca $ \fsp -> + alloca $ \fop -> + VS.unsafeWith (unsafeCoerce inBuf) $ \iPtr -> + VSM.unsafeWith (unsafeCoerce outBuf) $ \oPtr -> do+ c_dcBlocker (fromIntegral num) (realToFrac lastSample) (realToFrac lastOutput) fsp fop iPtr oPtr+ r1 <- peek fsp+ r2 <- peek fop+ return (realToFrac r1, realToFrac r2)
+ hs_sources/SDR/PipeUtils.hs view
@@ -0,0 +1,55 @@+{-| Pipes utility functions -}+module SDR.PipeUtils (+ fork, + combine,+ printStream,+ devnull,+ rate,+ ) where++import Data.Time.Clock+import Pipes+import Control.Monad++-- | Fork a pipe +fork :: Monad m => Producer a m r -> Producer a (Producer a m) r+fork prod = runEffect $ hoist (lift . lift) prod >-> fork' + where + fork' = forever $ do+ res <- await+ lift $ yield res+ lift $ lift $ yield res++-- | Combime two consumers into a single consumer+combine :: Monad m => Consumer a m r -> Consumer a m r -> Consumer a m r+combine x y = runEffect $ runEffect (fork func >-> hoist (lift . lift) x) >-> hoist lift y+ where+ func :: Monad m => Producer a (Consumer a m) r+ func = forever $ lift await >>= yield++-- | A consumer that prints everything to stdout+printStream :: (Show a) => Int -> Consumer a IO ()+printStream samples = for cat $ lift . print ++-- | A consumer that discards everything+devnull :: Monad m => Consumer a m ()+devnull = forever await++-- | Passthrough pipe that prints the sample rate+rate :: Int -> Pipe a a IO b+rate samples = do+ start <- lift getCurrentTime + let rate' buffers = do+ res <- await++ time <- lift getCurrentTime + let diff = diffUTCTime time start + diffSecs :: Double+ diffSecs = fromRational $ toRational diff++ lift $ print $ buffers * fromIntegral samples / diffSecs++ yield res+ rate' (buffers + 1)+ rate' 1+
+ hs_sources/SDR/Plot.hs view
@@ -0,0 +1,172 @@+{-| Create graphical plots of signals and their spectrums. Uses OpenGL. -}+module SDR.Plot (++ -- * Line Graphs+ plotLine,+ plotLineAxes,++ -- * Waterfalls+ plotWaterfall,+ --plotWaterfallAxes,++ -- * Filled In Line Graphs+ plotFill,+ plotFillAxes,++ -- * Axes+ zeroAxes,+ centeredAxes+ ) where++import Control.Monad.Trans.Either+import qualified Data.Vector.Storable as VS+import Graphics.Rendering.OpenGL+import Graphics.Rendering.Cairo++import Pipes +import Data.Colour.Names+import Graphics.Rendering.Pango++import Graphics.DynamicGraph.Line+import Graphics.DynamicGraph.Waterfall +import Graphics.DynamicGraph.FillLine +import Graphics.DynamicGraph.Axis+import Graphics.DynamicGraph.RenderCairo+import Graphics.DynamicGraph.Window++-- | Create a window and plot a dynamic line graph of the incoming data.+plotLine :: Int -- ^ Window width+ -> Int -- ^ Window height+ -> Int -- ^ Number of samples in each buffer+ -> Int -- ^ Number of vertices in graph+ -> EitherT String IO (Consumer (VS.Vector GLfloat) IO ())+plotLine width height samples resolution = window width height $ fmap pipeify $ renderLine samples resolution++-- | Create a window and plot a dynamic line graph of the incoming data. With Axes.+plotLineAxes :: Int -- ^ Window width+ -> Int -- ^ Window height+ -> Int -- ^ Number of samples in each buffer+ -> Int -- ^ Number of vertices in graph+ -> Render () -- ^ Cairo Render object that draws the axes+ -> EitherT String IO (Consumer (VS.Vector GLfloat) IO ())+plotLineAxes width height samples xResolution rm = window width height $ do+ --render the graph+ renderFunc <- renderLine samples xResolution++ --render the axes+ renderAxisFunc <- renderCairo rm width height++ return $ for cat $ \dat -> lift $ do+ blend $= Disabled++ viewport $= (Position 50 50, Size (fromIntegral width - 100) (fromIntegral height - 100))+ renderFunc dat++ blend $= Enabled+ blendFunc $= (SrcAlpha, OneMinusSrcAlpha)++ viewport $= (Position 0 0, Size (fromIntegral width) (fromIntegral height))+ renderAxisFunc++-- | Create a window and plot a waterfall of the incoming data.+plotWaterfall :: Int -- ^ Window width+ -> Int -- ^ Window height+ -> Int -- ^ Number of columns+ -> Int -- ^ Number of rows + -> [GLfloat] -- ^ The color map+ -> EitherT String IO (Consumer (VS.Vector GLfloat) IO ())+plotWaterfall windowWidth windowHeight width height colorMap = window windowWidth windowHeight $ renderWaterfall width height colorMap ++{-+-- | Create a window and plot a waterfall of the incoming data. With Axes. TODO: doesnt work.+plotWaterfallAxes :: Int -- ^ Window width + -> Int -- ^ Window height+ -> Int -- ^ Number of columns+ -> Int -- ^ Number of rows+ -> [GLfloat] -- ^ The color map+ -> Render () -- ^ Cairo Render object that draws the axes+ -> EitherT String IO (Consumer (VS.Vector GLfloat) IO ())+plotWaterfallAxes windowWidth windowHeight width height colorMap rm = window windowWidth windowHeight $ do+ renderPipe <- renderWaterfall width height colorMap+ + renderAxisFunc <- renderCairo rm width height++ return $ (>-> renderPipe) $ for cat $ \dat -> do+ lift $ viewport $= (Position 0 0, Size (fromIntegral width) (fromIntegral height))+ lift renderAxisFunc++ lift $ viewport $= (Position 50 50, Size (fromIntegral width - 100) (fromIntegral height - 100))++ yield dat+-}++-- | Create a window and plot a dynamic filled in line graph of the incoming data.+plotFill :: Int -- ^ Window width+ -> Int -- ^ Window height+ -> Int -- ^ Number of samples in each buffer+ -> [GLfloat] -- ^ The color map+ -> EitherT String IO (Consumer (VS.Vector GLfloat) IO ())+plotFill width height samples colorMap = window width height $ fmap pipeify $ renderFilledLine samples colorMap++-- | Create a window and plot a dynamic filled in line graph of the incoming data. With Axes.+plotFillAxes :: Int -- ^ Window width+ -> Int -- ^ Window height+ -> Int -- ^ Number of samples in each buffer+ -> [GLfloat] -- ^ The color map+ -> Render () -- ^ Cairo Render object that draws the axes+ -> EitherT String IO (Consumer (VS.Vector GLfloat) IO ())+plotFillAxes width height samples colorMap rm = window width height $ do+ renderFunc <- renderFilledLine samples colorMap+ + renderAxisFunc <- renderCairo rm width height++ return $ for cat $ \dat -> lift $ do+ viewport $= (Position 50 50, Size (fromIntegral width - 100) (fromIntegral height - 100))+ renderFunc dat++ blend $= Enabled+ blendFunc $= (SrcAlpha, OneMinusSrcAlpha)++ viewport $= (Position 0 0, Size (fromIntegral width) (fromIntegral height))+ renderAxisFunc++-- | Create a Cairo `Render` monad that draws a set of axes with 0 at the bottom left.+zeroAxes :: Int -- ^ Image width+ -> Int -- ^ Image height+ -> Double -- ^ X axis span+ -> Double -- ^ X axis grid interval+ -> Render ()+zeroAxes width height bandwidth interval = do+ blankCanvasAlpha black 0 (fromIntegral width) (fromIntegral height) + let xSeparation = (interval / bandwidth) * (fromIntegral width - 100)+ ySeparation = 0.2 * (fromIntegral height - 100)+ xCoords = takeWhile (< (fromIntegral width - 50)) $ iterate (+ xSeparation) 50+ yCoords = takeWhile (> 50) $ iterate (\x -> x - ySeparation) (fromIntegral height - 50)+ ctx <- liftIO $ cairoCreateContext Nothing+ xAxisLabels ctx white (map (\n -> show n ++ " KHz" ) (takeWhile (< bandwidth) $ iterate (+ interval) 0)) xCoords (fromIntegral height - 50)+ drawAxes (fromIntegral width) (fromIntegral height) 50 50 50 50 white 2+ xAxisGrid gray 1 [] 50 (fromIntegral height - 50) xCoords+ yAxisGrid gray 1 [4, 2] 50 (fromIntegral width - 50) yCoords++-- | Create a Cairo `Render` monad that draws a set of axes witb the X axis centered on a specified value.+centeredAxes :: Int -- ^ Image width+ -> Int -- ^ Image height+ -> Double -- ^ Center X value+ -> Double -- ^ X axis span+ -> Double -- ^ X axis grid interval+ -> Render ()+centeredAxes width height cFreq bandwidth interval = do+ blankCanvasAlpha black 0 (fromIntegral width) (fromIntegral height) + let xSeparation = (interval / bandwidth) * (fromIntegral width - 100)+ firstXLabel = fromIntegral (ceiling ((cFreq - (bandwidth / 2)) / interval)) * interval+ fract x = x - fromIntegral (floor x)+ xOffset = fract ((cFreq - (bandwidth / 2)) / interval) * xSeparation+ ySeparation = 0.2 * (fromIntegral height - 100)+ xCoords = takeWhile (< (fromIntegral width - 50)) $ iterate (+ xSeparation) (50 + xOffset)+ yCoords = takeWhile (> 50) $ iterate (\x -> x - ySeparation) (fromIntegral height - 50)+ ctx <- liftIO $ cairoCreateContext Nothing+ xAxisLabels ctx white (map (\n -> show n ++ " MHZ") (takeWhile (< (cFreq + bandwidth / 2)) $ iterate (+ interval) firstXLabel)) xCoords (fromIntegral height - 50)+ drawAxes (fromIntegral width) (fromIntegral height) 50 50 50 50 white 2+ xAxisGrid gray 1 [] 50 (fromIntegral height - 50) xCoords+ yAxisGrid gray 1 [4, 2] 50 (fromIntegral width - 50) yCoords+
+ hs_sources/SDR/Pulse.hs view
@@ -0,0 +1,33 @@+{-| Pulse Audio Pipes sink -}+module SDR.Pulse (+ pulseAudioSink,+ doPulse+ ) where++import Foreign.ForeignPtr+import Foreign.C.Types+import Control.Concurrent+import Data.ByteString.Internal+import Data.Vector.Storable as VS++import Sound.Pulse.Simple+import Pipes+import Pipes.Concurrent++-- | Returns a consumer that sends all incoming data to pulseaudio. Runs Pulse Audio output writing in a different thread. This is probably what you want as it does not block the entire pipline while the data is being played.+pulseAudioSink :: IO (Consumer (VS.Vector Float) IO ())+pulseAudioSink = do+ (output, input) <- spawn $ bounded 1+ doIt <- doPulse + forkOS $ runEffect $ fromInput input >-> doIt+ return $ toOutput output++-- | Returns a consumer that sends all incoming data to pulseaudio.+doPulse :: IO (Consumer (VS.Vector Float) IO ())+doPulse = do+ s <- simpleNew Nothing "Haskell SDR" Play Nothing "Software Defined Radio library" (SampleSpec (F32 LittleEndian) 48000 1) Nothing Nothing+ return $ for cat $ \buf -> + lift $ do+ let (fp, offset, length) = VS.unsafeToForeignPtr buf+ simpleWriteRaw s (PS (castForeignPtr fp) (offset * 4) (length * 4))+
+ hs_sources/SDR/RTLSDRStream.hs view
@@ -0,0 +1,53 @@+{-| Stream samples from a Realtek RTL2832U based device -}+module SDR.RTLSDRStream (+ sdrStream+ ) where++import Control.Monad+import Control.Monad.Trans.Either+import Data.Word+import Foreign.ForeignPtr+import Foreign.C.Types+import Control.Concurrent hiding (yield)+import Foreign.Marshal.Utils+import qualified Data.Vector.Storable as VS++import Pipes+import Pipes.Concurrent +import RTLSDR++-- | Returns a producer that streams data from a Realtek RTL2832U based device. You probably want to use `makeComplexBufferVect` to turn it into a list of complex Floats.+sdrStream :: Word32 -- ^ Frequency+ -> Word32 -- ^ Sample rate+ -> Word32 -- ^ Number of buffers+ -> Word32 -- ^ Buffer length+ -> EitherT String IO (Producer (VS.Vector CUChar) IO ()) -- ^ Either a string describing the error that occurred or the Producer+sdrStream frequency sampleRate bufNum bufLen = do+ lift $ putStrLn "Initializing RTLSDR device"++ dev' <- lift $ open 0+ dev <- maybe (left "Failed to open device") return dev'++ lift $ do+ t <- getTunerType dev+ putStrLn $ "Found a: " ++ show t++ setFreqCorrection dev 0+ setSampleRate dev sampleRate+ setCenterFreq dev frequency+ setTunerGainMode dev False++ resetBuffer dev++ (output, input) <- spawn unbounded++ forkOS $ void $ readAsync dev bufNum bufLen $ \dat num -> void $ do+ let numBytes = fromIntegral $ bufNum * bufLen+ fp <- mallocForeignPtrArray numBytes+ withForeignPtr fp $ \fpp -> moveBytes fpp dat numBytes+ let v = VS.unsafeFromForeignPtr0 fp numBytes+ atomically (send output v)++ return $ fromInput input++
+ hs_sources/SDR/Serialize.hs view
@@ -0,0 +1,84 @@+{-# LANGUAGE ScopedTypeVariables, FlexibleContexts #-}++{-| Utility functions for serializing and deserializing samples. -}+module SDR.Serialize (++ -- * Slow Serialization\/Deserialization+ -- | Slow functions for serializing\/deserializing vectors to\/from bytestrings using the Cereal library. There must be a better way to do this that doesn't involve copying.+ + -- ** Floats+ floatVecToByteString,+ floatVecFromByteString,++ -- ** Doubles+ doubleVecToByteString,+ doubleVecFromByteString,++ -- * Fast Serialization\/Deserialization+ -- | Fast functions for serializing\/deserializing storable vectors to\/from bytestrings.+ toByteString,+ fromByteString,++ -- * Pipes+ -- | Pipes that perform fast serialization/deserialization to a Handle.+ toHandle,+ fromHandle+ ) where++import Foreign.ForeignPtr+import Foreign.Storable+import Data.ByteString.Internal +import Data.ByteString as BS+import System.IO++import Data.Vector.Generic as VG hiding ((++))+import Data.Vector.Storable as VS hiding ((++))++import Pipes+import qualified Pipes.Prelude as P+import qualified Pipes.ByteString as PB+import Data.Serialize hiding (Done)+import qualified Data.Serialize as S++-- | Convert a Vector of Floats to a ByteString.+floatVecToByteString :: VG.Vector v Float => v Float -> ByteString+floatVecToByteString vect = runPut $ VG.mapM_ putFloat32le vect++-- | Convert a Vector of Doubles to a ByteString.+doubleVecToByteString :: VG.Vector v Double => v Double -> ByteString+doubleVecToByteString vect = runPut $ VG.mapM_ putFloat64le vect++-- | Convert a ByteString to a Vector of Floats.+floatVecFromByteString :: VG.Vector v Float => ByteString -> v Float+floatVecFromByteString bs = VG.unfoldrN (BS.length bs `div` 4) go bs+ where+ go bs = case runGetPartial getFloat32le bs of+ Fail _ _ -> Nothing+ Partial _ -> error "floatVecFromByteString: Partial"+ S.Done r b -> Just (r, b)++-- | Convert a ByteString to a Vector of Doubles.+doubleVecFromByteString :: VG.Vector v Double => ByteString -> v Double+doubleVecFromByteString bs = VG.unfoldrN (BS.length bs `div` 8) go bs+ where+ go bs = case runGetPartial getFloat64le bs of+ Fail _ _ -> Nothing+ Partial _ -> error "doubleVecFromByteString"+ S.Done r b -> Just (r, b)++-- | Convert a Vector of Storable values to a ByteString. This is fast as it is just a cast.+toByteString :: forall a. Storable a => VS.Vector a -> ByteString +toByteString dat = let (fp, o, sz) = VS.unsafeToForeignPtr dat in PS (castForeignPtr fp) o (sz * sizeOf (undefined :: a))++-- | Convert a ByteString to a Vector of Storable values. This is fast as it is just a cast.+fromByteString :: forall a. Storable a => ByteString -> VS.Vector a+fromByteString (PS fp o l) = VS.unsafeFromForeignPtr (castForeignPtr fp) o (l `quot` sizeOf (undefined :: a))++-- | Given a Handle, create a Consumer that dumps the Vectors written to it to a Handle.+toHandle :: (Storable a) => Handle -> Consumer (VS.Vector a) IO ()+toHandle handle = P.map toByteString >-> PB.toHandle handle ++-- | Given a Handle, create a Producer that creates Vectors from data read from the Handle.+fromHandle :: forall a. (Storable a) => Int -> Handle -> Producer (VS.Vector a) IO ()+fromHandle samples handle = PB.hGet (samples * sizeOf (undefined :: a)) handle >-> P.map fromByteString +
+ hs_sources/SDR/Util.hs view
@@ -0,0 +1,154 @@+{-# LANGUAGE FlexibleContexts, ScopedTypeVariables, MultiParamTypeClasses, FlexibleInstances #-}++{-| Various utiliy signal processing functions -}+module SDR.Util (+ -- * Classes+ Mult,+ mult,++ -- * Conversion to Floating Point+ makeComplexBufferVect,+ convertC, + convertCSSE,+ convertCAVX,++ -- * Scaling+ scaleC,+ scaleCSSE,+ scaleCAVX,++ -- * Misc Utils+ cplxMap,+ quarterBandUp+ ) where++import Foreign.C.Types+import Data.Complex+import Data.Vector.Generic as VG hiding ((++))+import qualified Data.Vector.Generic.Mutable as VGM+import Data.Vector.Storable as VS hiding ((++))+import Data.Vector.Storable.Mutable as VSM +import Control.Monad.Primitive+import Unsafe.Coerce+import Foreign.Ptr+import System.IO.Unsafe+import Foreign.Storable.Complex++-- | A class for things that can be multiplied by a scalar.+class Mult a b where+ mult :: a -> b -> a++instance (Num a) => Mult a a where+ mult = (*)++instance (Num a) => Mult (Complex a) a where+ mult (x :+ y) z = (x * z) :+ (y * z)++--TODO: none of these functions need the num argument++-- | Create a vector of complex float samples from a vector of interleaved I Q component bytes.+{-# INLINE makeComplexBufferVect #-}+makeComplexBufferVect :: (Num a, Integral a, Num b, Fractional b, VG.Vector v1 a, VG.Vector v2 (Complex b)) => v1 a -> v2 (Complex b)+makeComplexBufferVect input = VG.generate (VG.length input `quot` 2) convert+ where+ {-# INLINE convert #-}+ convert idx = convert' (input `VG.unsafeIndex` (2 * idx)) :+ convert' (input `VG.unsafeIndex` (2 * idx + 1))+ {-# INLINE convert' #-}+ convert' val = (fromIntegral val - 128) / 128++foreign import ccall unsafe "convertC"+ convertC_c :: CInt -> Ptr CUChar -> Ptr CFloat -> IO ()++-- | Same as `makeComplexBufferVect` but written in C and specialized for Floats+convertC :: VS.Vector CUChar -> VS.Vector (Complex Float)+convertC inBuf = unsafePerformIO $ do+ outBuf <- VGM.new $ VG.length inBuf `quot` 2+ VS.unsafeWith inBuf $ \iPtr -> + VSM.unsafeWith (unsafeCoerce outBuf) $ \oPtr -> + convertC_c (fromIntegral $ VG.length inBuf) iPtr oPtr+ VG.freeze outBuf++foreign import ccall unsafe "convertCSSE"+ convertCSSE_c :: CInt -> Ptr CUChar -> Ptr CFloat -> IO ()++-- | Same as `makeComplexBufferVect` but written in C using SSE intrinsics and specialized for Floats+convertCSSE :: VS.Vector CUChar -> VS.Vector (Complex Float)+convertCSSE inBuf = unsafePerformIO $ do+ outBuf <- VGM.new $ VG.length inBuf `quot` 2+ VS.unsafeWith inBuf $ \iPtr -> + VSM.unsafeWith (unsafeCoerce outBuf) $ \oPtr -> + convertCSSE_c (fromIntegral $ VG.length inBuf) iPtr oPtr+ VG.freeze outBuf++foreign import ccall unsafe "convertCAVX"+ convertCAVX_c :: CInt -> Ptr CUChar -> Ptr CFloat -> IO ()++-- | Same as `makeComplexBufferVect` but written in C using AVX intrinsics and specialized for Floats+convertCAVX :: VS.Vector CUChar -> VS.Vector (Complex Float)+convertCAVX inBuf = unsafePerformIO $ do+ outBuf <- VGM.new $ VG.length inBuf `quot` 2+ VS.unsafeWith inBuf $ \iPtr -> + VSM.unsafeWith (unsafeCoerce outBuf) $ \oPtr -> + convertCAVX_c (fromIntegral $ VG.length inBuf) iPtr oPtr+ VG.freeze outBuf++-- | Scaling+foreign import ccall unsafe "scale"+ scale_c :: CInt -> CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()++-- | Scale a vector, written in C+scaleC :: Float -- ^ Scale factor+ -> VS.Vector Float -- ^ Input vector+ -> VS.MVector RealWorld Float -- ^ Output vector+ -> IO ()+scaleC factor inBuf outBuf = + VS.unsafeWith (unsafeCoerce inBuf) $ \iPtr -> + VS.unsafeWith (unsafeCoerce outBuf) $ \oPtr -> + scale_c (fromIntegral (VG.length inBuf)) (unsafeCoerce factor) iPtr oPtr++foreign import ccall unsafe "scaleSSE"+ scaleSSE_c :: CInt -> CFloat -> Ptr CFloat -> Ptr CFloat-> IO ()++-- | Scale a vector, written in C using SSE intrinsics+scaleCSSE :: Float -- ^ Scale factor+ -> VS.Vector Float -- ^ Input vector+ -> VS.MVector RealWorld Float -- ^ Output vector+ -> IO ()+scaleCSSE factor inBuf outBuf = + VS.unsafeWith (unsafeCoerce inBuf) $ \iPtr -> + VS.unsafeWith (unsafeCoerce outBuf) $ \oPtr -> + scaleSSE_c (fromIntegral (VG.length inBuf)) (unsafeCoerce factor) iPtr oPtr++foreign import ccall unsafe "scaleAVX"+ scaleAVX_c :: CInt -> CFloat -> Ptr CFloat -> Ptr CFloat -> IO ()++-- | Scale a vector, written in C using AVX intrinsics+scaleCAVX :: Float -- ^ Scale factor+ -> VS.Vector Float -- ^ Input vector+ -> VS.MVector RealWorld Float -- ^ Output vector+ -> IO ()+scaleCAVX factor inBuf outBuf = + VS.unsafeWith (unsafeCoerce inBuf) $ \iPtr -> + VS.unsafeWith (unsafeCoerce outBuf) $ \oPtr -> + scaleAVX_c (fromIntegral (VG.length inBuf)) (unsafeCoerce factor) iPtr oPtr++-- | Apply a function to both parts of a complex number+cplxMap :: (a -> b) -- ^ The function+ -> Complex a -- ^ Input complex number+ -> Complex b -- ^ Output complex number+cplxMap f (x :+ y) = f x :+ f y++-- | Multiplication by this vector shifts all frequencies up by 1/4 of the sampling frequency+quarterBandUp :: (VG.Vector v (Complex n), Num n) + => Int -- ^ The length of the Vector+ -> v (Complex n)+quarterBandUp size = VG.generate size func+ where+ func idx + | m == 0 = 1 :+ 0+ | m == 1 = 0 :+ 1+ | m == 2 = (-1) :+ 0+ | m == 3 = 0 :+ (-1)+ where+ m = idx `mod` 4+
+ hs_sources/SDR/VectorUtils.hs view
@@ -0,0 +1,64 @@+{-# LANGUAGE BangPatterns #-}++{-| Various Vector based utility functions -}+module SDR.VectorUtils (+ mapAccumMV,+ stride,+ fill,+ ) where++import Control.Monad+import Control.Monad.Primitive++import Data.Vector.Generic as VG hiding ((++))+import qualified Data.Vector.Generic.Mutable as VGM+import Data.Vector.Fusion.Stream.Monadic hiding ((++))+import qualified Data.Vector.Fusion.Stream as VFS hiding ((++))+import qualified Data.Vector.Fusion.Stream.Monadic as VFSM hiding ((++))++{-| Like mapAccumL but monadic and over vectors. Doesn't return the+ accumulator at the end because it doesn't seem to be possible to do+ this with the Stream datatype, making this function pretty useless.+-}+mapAccumMV :: (Monad m) + => (acc -> x -> m (acc, y)) -- ^ The function+ -> acc -- ^ The initial accumulator+ -> Stream m x -- ^ The input stream+ -> Stream m y -- ^ The output stream+mapAccumMV func z (Stream step s sz) = Stream step' (s, z) sz+ where+ step' (s, acc) = do+ r <- step s+ case r of+ Yield y s' -> do+ (!acc', !res) <- func acc y + return $ Yield res (s', acc')+ Skip s' -> return $ Skip (s', acc)+ Done -> return Done++{-| Create a vector from another vector containing only the elements that+ occur every stride elements in the source vector.+-}+{-# INLINE stride #-}+stride :: VG.Vector v a + => Int -- ^ The stride+ -> v a -- ^ The input Vector+ -> v a -- ^ The output Vector+stride str inv = VG.unstream $ VFS.unfoldr func 0+ where+ len = VG.length inv+ func i | i >= len = Nothing+ | otherwise = Just (VG.unsafeIndex inv i, i + str)++-- | Fill a mutable vector from a monadic stream+{-# INLINE fill #-}+fill :: (PrimMonad m, Functor m, VGM.MVector vm a) + => VFS.MStream m a -- ^ The input Stream+ -> vm (PrimState m) a -- ^ The mutable Vector to stream into+ -> m ()+fill str outBuf = void $ VFSM.foldM' put 0 str+ where + put i x = do+ VGM.unsafeWrite outBuf i x+ return $ i + 1+
+ sdr.cabal view
@@ -0,0 +1,136 @@+name: sdr+version: 0.1.0.0+synopsis: A software defined radio library+description: + Write software defined radio applications in Haskell.+ .+ Features:+ .+ * Signal processing blocks can be chained together using the <https://hackage.haskell.org/package/pipes Pipes> library+ .+ * Zero copy design+ .+ * Signal processing functions are implemented in both Haskell and C (with SIMD acceleration)+ .+ * Can FIR filter, decimate and resample+ .+ * Helper functions for FIR filter design using window functions and plotting of the frequency response+ .+ * FFTs using <http://www.fftw.org/ FFTW>+ .+ * Line and waterfall plots using OpenGL+ .+ * FM demodulation+ .+ * PulseAudio sound sink+ .+ * <http://sdr.osmocom.org/trac/wiki/rtl-sdr rtl-sdr> based radio source supported and other sources are easily added+ .+ See <https://github.com/adamwalker/sdr> for more features and screenshots.+ .+ A collection of simple apps that use this library can be found <https://github.com/adamwalker/sdr-apps here>. These include an FM radio receiver, an OpenGL waterfall plotter and an AM radio receiver.++license: BSD3+license-file: LICENSE+author: Adam Walker+maintainer: adamwalker10@gmail.com+copyright: 2015 Adam Walker+category: Software Defined Radio+homepage: https://github.com/adamwalker/sdr+bug-reports: https://github.com/adamwalker/sdr/issues+build-type: Simple+-- extra-source-files: +cabal-version: >=1.10++source-repository head+ type: git+ location: https://github.com/adamwalker/sdr++library+ exposed-modules: + SDR.Pulse, + SDR.RTLSDRStream, + SDR.Util, + SDR.Plot, + SDR.Filter, + SDR.Demod, + SDR.FFT, + SDR.FilterInternal, + SDR.Serialize, + SDR.PipeUtils, + SDR.VectorUtils, + SDR.ArgUtils, + SDR.FilterDesign+ -- other-modules: + other-extensions: ScopedTypeVariables, GADTs+ build-depends: + base >=4.6 && <4.9, + fftwRaw >=0.1 && <0.2, + bytestring >=0.10 && <0.11, + pulse-simple >=0.1 && <0.2, + pipes >=4.1 && <4.2, + pipes-concurrency >=2.0 && <2.1, + either >=4.1 && <4.4, + time >=1.4 && <1.6, + rtlsdr >=0.1 && <0.2, + storable-complex >=0.2 && <0.3, + pipes-bytestring >=2.0 && <2.2, + dynamic-graph ==0.1.0.8,+ array >=0.4 && <0.6, + vector >=0.10 && <0.11,+ tuple >=0.2 && <0.4, + OpenGL >=2.11 && <2.13, + GLFW-b >=1.4.7 && <1.4.8,+ primitive >=0.5 && <0.7, + colour >=2.3 && <2.4, + pango >=0.13 && <0.14, + containers >=0.5 && <0.6, + cairo >=0.13 && <0.14, + cereal >=0.4 && <0.5, + optparse-applicative >=0.11 && <0.12, + Decimal >=0.4 && <0.5, + Chart >=1.3 && <1.5, + Chart-cairo >=1.3 && <1.5+ -- hs-source-dirs: + default-language: Haskell2010+ ghc-options: -O2+ includes: c_sources/common.h+ c-sources: + c_sources/filter.c, + c_sources/decimate.c, + c_sources/convert.c, + c_sources/resample.c, + c_sources/scale.c+ hs-source-dirs: hs_sources+ cc-options: -mavx2 -msse4 -g++Test-Suite test+ type: exitcode-stdio-1.0+ main-is: TestSuite.hs+ build-depends: + base >=4.6 && <4.9, + QuickCheck >=2.8 && <2.9, + vector >=0.10 && <0.11, + sdr ==0.1.0.0, + primitive >=0.5 && <0.7, + storable-complex >=0.2 && <0.3,+ test-framework >=0.8 && <0.9,+ test-framework-quickcheck2 >=0.3 && <0.4+ hs-source-dirs: tests+ ghc-options: -O2+ default-language: Haskell2010++benchmark benchmark+ type: exitcode-stdio-1.0+ main-is: Benchmarks.hs+ build-depends: + base >=4.6 && <4.9, + criterion >=1.0 && <1.2,+ vector >=0.10 && <0.11, + sdr ==0.1.0.0, + primitive >=0.5 && <0.7, + storable-complex >=0.2 && <0.3+ hs-source-dirs: benchmarks+ ghc-options: -O2+ default-language: Haskell2010+
+ tests/TestSuite.hs view
@@ -0,0 +1,213 @@+{-# LANGUAGE ScopedTypeVariables #-}++import Control.Monad.Primitive +import Data.Complex++import qualified Data.Vector.Generic as VG+import qualified Data.Vector.Generic.Mutable as VGM+import qualified Data.Vector.Storable as VS+import qualified Data.Vector.Storable.Mutable as VSM++import Test.QuickCheck+import Test.QuickCheck.Monadic+import Test.Framework (defaultMain, testGroup)+import Test.Framework.Providers.QuickCheck2 (testProperty)++import SDR.FilterInternal+import SDR.Util++tests = [+ testGroup "filters" [+ testProperty "real" propFiltersReal,+ testProperty "complex" propFiltersComplex+ ],+ testGroup "decimators" [+ testProperty "real" propDecimationReal,+ testProperty "complex" propDecimationComplex+ ],+ testGroup "resamplers" [+ testProperty "real" propResamplingReal+ ],+ testProperty "conversion" propConversion,+ testProperty "scaling" propScaleReal+ ]+ where+ sizes = elements [1024, 2048, 4096, 8192, 16384, 32768, 65536]+ numCoeffs = elements [32, 64, 128, 256, 512]+ factors = elements [1, 2, 3, 4, 7, 9, 12, 15, 21]+ factors' = [1, 2, 3, 4, 7, 9, 12, 15, 21]++ propFiltersReal = forAll sizes $ \size -> + forAll (vectorOf size (choose (-10, 10))) $ \inBuf -> + forAll numCoeffs $ \numCoeffs -> + forAll (vectorOf numCoeffs (choose (-10, 10))) $ \coeffs -> + testFiltersReal size numCoeffs coeffs inBuf++ testFiltersReal :: Int -> Int -> [Float] -> [Float] -> Property+ testFiltersReal size numCoeffs coeffs inBuf = monadicIO $ do+ let vCoeffsHalf = VS.fromList coeffs+ vCoeffs = VS.fromList $ coeffs ++ reverse coeffs+ vInput = VS.fromList inBuf+ num = size - numCoeffs*2 + 1++ r1 <- run $ getResult num $ filterHighLevel vCoeffs num vInput+ r2 <- run $ getResult num $ filterImperative1 vCoeffs num vInput+ r3 <- run $ getResult num $ filterImperative2 vCoeffs num vInput+ r4 <- run $ getResult num $ filterCRR vCoeffs num vInput+ r5 <- run $ getResult num $ filterCSSERR vCoeffs num vInput+ r6 <- run $ getResult num $ filterCAVXRR vCoeffs num vInput+ r7 <- run $ getResult num $ filterCSSESymmetricRR vCoeffsHalf num vInput+ r8 <- run $ getResult num $ filterCAVXSymmetricRR vCoeffsHalf num vInput++ assert $ and $ map (r1 `eqDelta`) [r2, r3, r4, r5, r6, r7, r8]++ propFiltersComplex = forAll sizes $ \size -> + forAll (vectorOf size (choose (-10, 10))) $ \inBufR -> + forAll (vectorOf size (choose (-10, 10))) $ \inBufI -> + forAll numCoeffs $ \numCoeffs -> + forAll (vectorOf numCoeffs (choose (-10, 10))) $ \coeffs -> + testFiltersComplex size numCoeffs coeffs $ zipWith (:+) inBufR inBufI++ testFiltersComplex :: Int -> Int -> [Float] -> [Complex Float] -> Property+ testFiltersComplex size numCoeffs coeffs inBuf = monadicIO $ do+ let vCoeffsHalf = VS.fromList coeffs+ vCoeffs = VS.fromList $ coeffs ++ reverse coeffs+ vInput = VS.fromList inBuf+ num = size - numCoeffs*2 + 1+ vCoeffs2 = VG.fromList $ duplicate $ coeffs ++ reverse coeffs++ r1 <- run $ getResult num $ filterHighLevel vCoeffs num vInput+ r2 <- run $ getResult num $ filterCRC vCoeffs num vInput+ r3 <- run $ getResult num $ filterCSSERC vCoeffs2 num vInput+ r4 <- run $ getResult num $ filterCSSERC2 vCoeffs num vInput+ r5 <- run $ getResult num $ filterCAVXRC vCoeffs2 num vInput+ r6 <- run $ getResult num $ filterCSSESymmetricRC vCoeffsHalf num vInput+ r7 <- run $ getResult num $ filterCAVXSymmetricRC vCoeffsHalf num vInput++ assert $ and $ map (r1 `eqDeltaC`) [r2, r3, r4, r5, r6, r7]++ propDecimationReal = forAll sizes $ \size -> + forAll (vectorOf size (choose (-10, 10))) $ \inBuf -> + forAll numCoeffs $ \numCoeffs -> + forAll (vectorOf numCoeffs (choose (-10, 10))) $ \coeffs -> + forAll factors $ \factor -> + testDecimationReal size numCoeffs factor coeffs inBuf++ testDecimationReal :: Int -> Int -> Int -> [Float] -> [Float] -> Property+ testDecimationReal size numCoeffs factor coeffs inBuf = monadicIO $ do+ let vCoeffsHalf = VS.fromList coeffs+ vCoeffs = VS.fromList $ coeffs ++ reverse coeffs+ vInput = VS.fromList inBuf+ num = (size - numCoeffs*2 + 1) `quot` factor++ r1 <- run $ getResult num $ decimateHighLevel factor vCoeffs num vInput+ r2 <- run $ getResult num $ decimateCRR factor vCoeffs num vInput+ r3 <- run $ getResult num $ decimateCSSERR factor vCoeffs num vInput+ r4 <- run $ getResult num $ decimateCAVXRR factor vCoeffs num vInput+ r5 <- run $ getResult num $ decimateCSSESymmetricRR factor vCoeffsHalf num vInput+ r6 <- run $ getResult num $ decimateCAVXSymmetricRR factor vCoeffsHalf num vInput++ assert $ and $ map (r1 `eqDelta`) [r2, r3, r4, r5, r6]++ propDecimationComplex = forAll sizes $ \size -> + forAll (vectorOf size (choose (-10, 10))) $ \inBufR -> + forAll (vectorOf size (choose (-10, 10))) $ \inBufI -> + forAll numCoeffs $ \numCoeffs -> + forAll (vectorOf numCoeffs (choose (-10, 10))) $ \coeffs -> + forAll factors $ \factor -> + testDecimationComplex size numCoeffs factor coeffs $ zipWith (:+) inBufR inBufI++ testDecimationComplex :: Int -> Int -> Int -> [Float] -> [Complex Float] -> Property+ testDecimationComplex size numCoeffs factor coeffs inBuf = monadicIO $ do+ let vCoeffsHalf = VS.fromList coeffs+ vCoeffs = VS.fromList $ coeffs ++ reverse coeffs+ vInput = VS.fromList inBuf+ num = (size - numCoeffs*2 + 1) `quot` factor+ vCoeffs2 = VG.fromList $ duplicate $ coeffs ++ reverse coeffs++ r1 <- run $ getResult num $ decimateHighLevel factor vCoeffs num vInput+ r2 <- run $ getResult num $ decimateCRC factor vCoeffs num vInput+ r3 <- run $ getResult num $ decimateCSSERC factor vCoeffs2 num vInput+ r4 <- run $ getResult num $ decimateCSSERC2 factor vCoeffs num vInput+ r5 <- run $ getResult num $ decimateCAVXRC factor vCoeffs2 num vInput+ r6 <- run $ getResult num $ decimateCSSESymmetricRC factor vCoeffsHalf num vInput+ r7 <- run $ getResult num $ decimateCAVXRC2 factor vCoeffs num vInput+ r8 <- run $ getResult num $ decimateCAVXSymmetricRC factor vCoeffsHalf num vInput++ assert $ and $ map (r1 `eqDeltaC`) [r2, r3, r4, r5, r6, r7, r8]++ propResamplingReal = forAll sizes $ \size -> + forAll (vectorOf size (choose (-10, 10))) $ \inBuf -> + forAll numCoeffs $ \numCoeffs -> + forAll (vectorOf numCoeffs (choose (-10, 10))) $ \coeffs -> + forAll (elements $ tail factors') $ \decimation -> + forAll (elements $ filter (< decimation) factors') $ \interpolation -> + testResamplingReal size numCoeffs interpolation decimation coeffs inBuf++ testResamplingReal :: Int -> Int -> Int -> Int -> [Float] -> [Float] -> Property+ testResamplingReal size numCoeffs interpolation decimation coeffs inBuf = monadicIO $ do+ let vCoeffsHalf = VS.fromList coeffs+ vCoeffs = VS.fromList $ coeffs ++ reverse coeffs+ vInput = VS.fromList inBuf+ num = (size - numCoeffs*2 + 1) `quot` decimation++ resampler3 <- run $ resampleCRR2 interpolation decimation (coeffs ++ reverse coeffs)+ resampler4 <- run $ resampleCSSERR interpolation decimation (coeffs ++ reverse coeffs)+ resampler5 <- run $ resampleCAVXRR interpolation decimation (coeffs ++ reverse coeffs)++ r1 <- run $ getResult num $ resampleHighLevel interpolation decimation vCoeffs 0 num vInput+ r2 <- run $ getResult num $ resampleCRR num interpolation decimation 0 vCoeffs vInput+ r3 <- run $ getResult num $ resampler3 num 0 vInput+ r4 <- run $ getResult num $ resampler4 num 0 vInput+ r5 <- run $ getResult num $ resampler5 num 0 vInput++ assert $ and $ map (r1 `eqDelta`) [r2, r3, r4, r5]++ getResult :: (VSM.Storable a) => Int -> (VS.MVector RealWorld a -> IO b) -> IO [a]+ getResult size func = do+ outBuf <- VGM.new size+ func outBuf+ out :: VS.Vector a <- VG.freeze outBuf+ return $ VG.toList out++ propConversion = forAll sizes $ \size -> + forAll (vectorOf (2 * size) (choose (-10, 10))) $ \inBuf -> + testConversion size inBuf+ testConversion :: Int -> [Int] -> Property+ testConversion size inBuf = monadicIO $ do+ let vInput = VS.fromList $ map fromIntegral inBuf++ let r1 = VG.toList $ (makeComplexBufferVect vInput :: VS.Vector (Complex Float))+ r2 = VG.toList $ convertC vInput+ r3 = VG.toList $ convertCSSE vInput+ r4 = VG.toList $ convertCAVX vInput++ assert $ and $ map (r1 `eqDeltaC`) [r2, r3, r4]++ scales = elements [0.1, 0.5, 1, 2, 10]+ propScaleReal = forAll sizes $ \size -> + forAll (vectorOf size (choose (-10, 10))) $ \inBuf -> + forAll scales $ \factor -> + testScaleReal size inBuf factor+ testScaleReal :: Int -> [Float] -> Float -> Property+ testScaleReal size inBuf factor = monadicIO $ do+ let vInput = VS.fromList inBuf++ r1 <- run $ getResult size $ scaleC factor vInput+ r2 <- run $ getResult size $ scaleCSSE factor vInput+ r3 <- run $ getResult size $ scaleCAVX factor vInput++ assert $ and $ map (r1 `eqDelta`) [r2, r3]++ eqDelta x y = and $ map (uncurry eqDelta') $ zip x y+ where+ eqDelta' x y = abs (x - y) < 0.01+ eqDeltaC x y = and $ map (uncurry eqDelta') $ zip x y+ where+ eqDelta' x y = magnitude (x - y) < 0.01+ duplicate :: [a] -> [a]+ duplicate = concat . map func + where func x = [x, x]++main = defaultMain tests+