diff --git a/Readme.md b/Readme.md new file mode 100644 --- /dev/null +++ b/Readme.md @@ -0,0 +1,124 @@ +# SDR + +A Software Defined Radio library written in Haskell + +# Features +* Write software defined radio applications in Haskell +* Signal processing blocks can be chained together using the [Pipes](https://hackage.haskell.org/package/pipes) library +* Zero copy design +* Signal processing functions are implemented in both Haskell and C: + * Optimised C implementations of signal processing functions that utilise SIMD instructions + * Performance of Haskell signal processing functions within a factor of 2 of C (without SIMD) thanks to the vector library, stream fusion and ghc's LLVM backend +* Can filter, decimate and resample +* Helper functions for FIR filter design using window functions and plotting of the frequency response +* FFTs using [FFTW](http://www.fftw.org) +* Line and waterfall plots using OpenGL +* FM demodulation +* PulseAudio sound sink +* [rtl-sdr](http://sdr.osmocom.org/trac/wiki/rtl-sdr) based radio source supported and other sources are easily added +* Extensive benchmark and test suites of signal processing functions + +See https://github.com/adamwalker/sdr-apps for a collection of simple apps built on the library and https://github.com/adamwalker/sdr-demo for a demo application. + +# Screenshot +A chunk of the FM broadcast spectrum. Captured with an RTLSDR device and drawn as a waterfall using the [Plot](https://github.com/adamwalker/sdr/blob/master/hs_sources/SDR/Plot.hs) module. + +![Screenshot](../screenshots/screenshots/screenshot.png?raw=true) + + +# Getting Started + +## Installation +Building with cabal sandboxes is recommended: + +``` +cabal sandbox init +git clone https://github.com/adamwalker/dynamic-graph +git clone https://github.com/adamwalker/haskell-fftw-simple +git clone https://github.com/adamwalker/sdr +cabal sandbox add-source dynamic-graph haskell-fftw-simple sdr +cabal install sdr +``` + +## Example Applications + +A collection of simple apps can be found [here](https://github.com/adamwalker/sdr-apps). These include an FM radio receiver, an OpenGL waterfall plotter and an AM radio receiver that can be used to listen to [Airband](https://en.wikipedia.org/wiki/Airband). + +Clone and build: + +``` +git clone https://github.com/adamwalker/sdr-apps +cabal sandbox add-source sdr-apps +cabal install sdr-apps +``` + +To run the FM receiver: +``` +.cabal-sandbox/bin/fm -f +``` + +To run the waterfall plot: +``` +.cabal-sandbox/bin/waterfall -f
-r +``` + +To run the AM receiver: +``` +.cabal-sandbox/bin/am -f
+``` + +# Usage + +An FM receiver: + +```haskell +import Control.Monad.Trans.Either +import Data.Vector.Generic as VG +import Pipes +import qualified Pipes.Prelude as P +import Foreign.Storable.Complex + +import SDR.Filter +import SDR.RTLSDRStream +import SDR.Util +import SDR.Demod +import SDR.Pulse +import SDR.CPUID + +--The filter coefficients are stored in another module +import Coeffs + +samples = 8192 +frequency = 105700000 + +main = eitherT putStrLn return $ do + + info <- lift getCPUInfo + + str <- sdrStream frequency 1280000 1 (fromIntegral samples * 2) + + lift $ do + + sink <- pulseAudioSink + + deci <- fastDecimatorC info 8 coeffsRFDecim + resp <- fastResamplerR info 3 10 coeffsAudioResampler + filt <- fastFilterSymR info coeffsAudioFilter + + runEffect $ str + >-> P.map convertCAVX + >-> firDecimator deci samples + >-> fmDemod + >-> firResampler resp samples + >-> firFilter filt samples + >-> P.map (VG.map (* 0.2)) + >-> sink +``` + +# Disclaimer +I started this project to learn about signal processing. I still have no idea what I'm doing. + +Only tested on Arch Linux. + +If you actually use this library for anything, let me know: adamwalker10@gmail.com + diff --git a/benchmarks/Benchmarks.hs b/benchmarks/Benchmarks.hs --- a/benchmarks/Benchmarks.hs +++ b/benchmarks/Benchmarks.hs @@ -1,4 +1,4 @@ -{-# LANGUAGE ScopedTypeVariables, BangPatterns #-} +{-# LANGUAGE ScopedTypeVariables #-} import Control.Monad.Primitive import Control.Monad @@ -19,6 +19,7 @@ import SDR.FilterInternal import SDR.Util +import SDR.CPUID theBench :: IO () theBench = do @@ -57,6 +58,11 @@ outBuf :: VS.MVector RealWorld Float <- VGM.new size outBufComplex :: VS.MVector RealWorld (Complex Float) <- VGM.new size + info <- getCPUInfo + + let hasFeatures :: [(CPUInfo -> Bool, a)] -> [a] + hasFeatures = map snd . filter (($ info) . fst) + resampler3 <- resampleCRR2 interpolation decimation coeffsList resampler4 <- resampleCSSERR interpolation decimation coeffsList resampler5 <- resampleCAVXRR interpolation decimation coeffsList @@ -64,63 +70,63 @@ --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 "real" $ hasFeatures [ + (const True, bench "highLevel" $ nfIO $ filterHighLevel coeffs num inBuf outBuf), + (const True, bench "imperative1" $ nfIO $ filterImperative1 coeffs num inBuf outBuf), + (const True, bench "imperative2" $ nfIO $ filterImperative2 coeffs num inBuf outBuf), + (const True, bench "c" $ nfIO $ filterCRR coeffs num inBuf outBuf), + (hasSSE42, bench "cSSE" $ nfIO $ filterCSSERR coeffs num inBuf outBuf), + (hasSSE42, bench "cSSESym" $ nfIO $ filterCSSESymmetricRR coeffsSym num inBuf outBuf), + (hasAVX, bench "cAVX" $ nfIO $ filterCAVXRR coeffs num inBuf outBuf), + (hasAVX, 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 "complex" $ hasFeatures [ + (const True, bench "highLevel" $ nfIO $ filterHighLevel coeffs num inBufComplex outBufComplex), + (const True, bench "c" $ nfIO $ filterCRC coeffs num inBufComplex outBufComplex), + (hasSSE42, bench "cSSE" $ nfIO $ filterCSSERC coeffs2 num inBufComplex outBufComplex), + (hasSSE42, bench "cSSE2" $ nfIO $ filterCSSERC2 coeffs num inBufComplex outBufComplex), + (hasSSE42, bench "cSSESym" $ nfIO $ filterCSSESymmetricRC coeffsSym num inBufComplex outBufComplex), + (hasAVX, bench "cAVX" $ nfIO $ filterCAVXRC coeffs2 num inBufComplex outBufComplex), + (hasAVX, bench "cAVX2" $ nfIO $ filterCAVXRC2 coeffs num inBufComplex outBufComplex), + (hasAVX, 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 "real" $ hasFeatures [ + (const True, bench "highLevel" $ nfIO $ decimateHighLevel decimation coeffs (num `quot` decimation) inBuf outBuf), + (const True, bench "c" $ nfIO $ decimateCRR decimation coeffs (num `quot` decimation) inBuf outBuf), + (hasSSE42, bench "cSSE" $ nfIO $ decimateCSSERR decimation coeffs (num `quot` decimation) inBuf outBuf), + (hasSSE42, bench "cSSESym" $ nfIO $ decimateCSSESymmetricRR decimation coeffsSym (num `quot` decimation) inBuf outBuf), + (hasAVX, bench "cAVX" $ nfIO $ decimateCAVXRR decimation coeffs (num `quot` decimation) inBuf outBuf), + (hasAVX, 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 "complex" $ hasFeatures [ + (const True, bench "highLevel" $ nfIO $ decimateHighLevel decimation coeffs (num `quot` decimation) inBufComplex outBufComplex), + (const True, bench "c" $ nfIO $ decimateCRC decimation coeffs (num `quot` decimation) inBufComplex outBufComplex), + (hasSSE42, bench "cSSE" $ nfIO $ decimateCSSERC decimation coeffs2 (num `quot` decimation) inBufComplex outBufComplex), + (hasSSE42, bench "cSSE2" $ nfIO $ decimateCSSERC2 decimation coeffs (num `quot` decimation) inBufComplex outBufComplex), + (hasSSE42, bench "cSSESym" $ nfIO $ decimateCSSESymmetricRC decimation coeffsSym (num `quot` decimation) inBufComplex outBufComplex), + (hasAVX, bench "cAVX" $ nfIO $ decimateCAVXRC decimation coeffs2 (num `quot` decimation) inBufComplex outBufComplex), + (hasAVX, bench "cAVX2" $ nfIO $ decimateCAVXRC2 decimation coeffs (num `quot` decimation) inBufComplex outBufComplex), + (hasAVX, 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 "real" $ hasFeatures [ + (const True, bench "highLevel" $ nfIO $ resampleHighLevel interpolation decimation coeffs 0 (num `quot` decimation) inBuf outBuf), + (const True, bench "c" $ nfIO $ resampleCRR (num `quot` decimation) interpolation decimation 0 coeffs inBuf outBuf), + (const True, bench "c2" $ nfIO $ resampler3 (num `quot` decimation) 0 inBuf outBuf), + (hasSSE42, bench "cSSE" $ nfIO $ resampler4 (num `quot` decimation) 0 inBuf outBuf), + (hasAVX, 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 "complex" $ hasFeatures [ + (const True, 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 + bgroup "scaling" $ hasFeatures [ + (const True, bench "c" $ nfIO $ scaleC 0.3 inBuf outBuf), + (hasSSE42, bench "cSSE" $ nfIO $ scaleCSSE 0.3 inBuf outBuf), + (hasAVX, bench "cAVX" $ nfIO $ scaleCAVX 0.3 inBuf outBuf) ] ] diff --git a/c_sources/cpuid.c b/c_sources/cpuid.c new file mode 100644 --- /dev/null +++ b/c_sources/cpuid.c @@ -0,0 +1,17 @@ +#include + +void cpuid(uint32_t op, uint32_t *a, uint32_t *b, uint32_t *c, uint32_t *d){ + asm volatile( + "cpuid;" + : "=a"(*a), "=b"(*b), "=c"(*c), "=d"(*d) + : "a"(op) + ); +} + +void cpuid_extended(uint32_t op, uint32_t sub_op, uint32_t *a, uint32_t *b, uint32_t *c, uint32_t *d){ + asm volatile( + "cpuid;" + : "=a"(*a), "=b"(*b), "=c"(*c), "=d"(*d) + : "a"(op), "c"(sub_op) + ); +} diff --git a/hs_sources/SDR/CPUID.hs b/hs_sources/SDR/CPUID.hs new file mode 100644 --- /dev/null +++ b/hs_sources/SDR/CPUID.hs @@ -0,0 +1,104 @@ +{-# LANGUAGE RecordWildCards #-} + +{-| This module is for detecting which SIMD instruction sets your CPU supports. In particular, it can detect SSE4.2, AVX and AVX2. -} +module SDR.CPUID ( + -- * Raw CPUID + cpuid, + cpuidExtended, + + -- * High level CPU capabilities + CPUInfo(..), + getCPUInfo, + + -- * Features + hasSSE42, + hasAVX, + hasAVX2, + + -- * Convenience functions + featureSelect + ) where + +import Data.Word +import Data.Bits +import Foreign.Ptr +import Foreign.Marshal.Alloc +import Foreign.Storable +import Data.List +import Data.Maybe + +foreign import ccall unsafe "cpuid" + cpuid_c :: Word32 -> Ptr Word32 -> Ptr Word32 -> Ptr Word32 -> Ptr Word32 -> IO () + +foreign import ccall unsafe "cpuid_extended" + cpuidExtended_c :: Word32 -> Word32 -> Ptr Word32 -> Ptr Word32 -> Ptr Word32 -> Ptr Word32 -> IO () + +-- | Execute the CPUID instruction +cpuid :: Word32 -- ^ Operation (EAX) + -> IO (Word32, Word32, Word32, Word32) -- ^ Result (EAX, EBX, ECX, EDX) +cpuid x = + alloca $ \p1 -> + alloca $ \p2 -> + alloca $ \p3 -> + alloca $ \p4 -> do + cpuid_c x p1 p2 p3 p4 + (,,,) <$> peek p1 <*> peek p2 <*> peek p3 <*> peek p4 + +-- | Execute the CPUID instruction setting ECX as well +cpuidExtended :: Word32 -- ^ Operation (EAX) + -> Word32 -- ^ ECX + -> IO (Word32, Word32, Word32, Word32) -- ^ Result (EAX, EBX, ECX, EDX) +cpuidExtended x y = + alloca $ \p1 -> + alloca $ \p2 -> + alloca $ \p3 -> + alloca $ \p4 -> do + cpuidExtended_c x y p1 p2 p3 p4 + (,,,) <$> peek p1 <*> peek p2 <*> peek p3 <*> peek p4 + +-- | Information about the features supported by your CPU +data CPUInfo = CPUInfo { + features :: Word32, + extendedFeatures :: Maybe Word32 +} + +-- | Get a `CPUInfo` +getCPUInfo :: IO CPUInfo +getCPUInfo = do + (x, _, _, _) <- cpuid 0 + (_, _, f, _) <- cpuid 1 + if x < 7 then + return $ CPUInfo f Nothing + else do + (_, e, _, _) <- cpuidExtended 7 0 + return $ CPUInfo f (Just e) + +-- | Feature bit for SSE4.2 +sse42 = 20 + +-- | Feature bit for AVX +avx = 28 + +-- | Extended feature bit for AVX2 +avx2 = 5 + +-- | Check if the CPU supports SSE4.2 +hasSSE42 :: CPUInfo -> Bool +hasSSE42 CPUInfo{..} = testBit features sse42 + +-- | Check if the CPU supports AVX +hasAVX :: CPUInfo -> Bool +hasAVX CPUInfo{..} = testBit features avx + +-- | Check if the CPU supports AVX2 +hasAVX2 :: CPUInfo -> Bool +hasAVX2 (CPUInfo _ Nothing) = False +hasAVX2 (CPUInfo _ (Just f)) = testBit f avx2 + +-- | Convenience function for selecting a function based on the features that the CPU supports +featureSelect :: CPUInfo -- ^ The CPU features + -> a -- ^ Default implementation + -> [(CPUInfo -> Bool, a)] -- ^ List of (feature, implementation) pairs + -> a -- ^ The selected implementation +featureSelect info def list = maybe def snd $ find (($ info) . fst) list + diff --git a/hs_sources/SDR/Demod.hs b/hs_sources/SDR/Demod.hs --- a/hs_sources/SDR/Demod.hs +++ b/hs_sources/SDR/Demod.hs @@ -33,6 +33,7 @@ -> v a -- ^ The output Vector fmDemodVec init = unstream . fmDemodStr init . stream +-- | Pipe that performs FM demodulation {-# INLINE fmDemod #-} fmDemod :: (RealFloat a, Vector v (Complex a), Vector v a) => Pipe (v (Complex a)) (v a) IO () fmDemod = func 0 diff --git a/hs_sources/SDR/Filter.hs b/hs_sources/SDR/Filter.hs --- a/hs_sources/SDR/Filter.hs +++ b/hs_sources/SDR/Filter.hs @@ -1,4 +1,4 @@ -{-# LANGUAGE RecordWildCards, FlexibleContexts, GADTs #-} +{-# LANGUAGE RecordWildCards, FlexibleContexts, GADTs, ExistentialQuantification #-} {-| FIR filtering, decimation and resampling. @@ -8,11 +8,11 @@ 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 + > decimatorStruct <- fastDecimatorC cpuInfo 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`. + 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 `fastFilterSymR`. 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. @@ -29,20 +29,53 @@ -- * Helper Functions -- ** Filters haskellFilter, + + -- *** Real Data + fastFilterCR, + fastFilterSSER, + fastFilterAVXR, fastFilterR, + + -- *** Complex Data + fastFilterCC, + fastFilterSSEC, + fastFilterAVXC, fastFilterC, - fastSymmetricFilterR, + -- *** Linear Phase Real Data + fastFilterSymSSER, + fastFilterSymAVXR, + fastFilterSymR, + -- ** Decimators haskellDecimator, + + -- *** Real Data + fastDecimatorCR, + fastDecimatorSSER, + fastDecimatorAVXR, fastDecimatorR, + + -- *** Complex Data + fastDecimatorCC, + fastDecimatorSSEC, + fastDecimatorAVXC, fastDecimatorC, - fastSymmetricDecimatorR, + -- *** Linear Phase Real Data + fastDecimatorSymSSER, + fastDecimatorSymAVXR, + fastDecimatorSymR, + -- ** Resamplers haskellResampler, - --fastResampler, + -- *** Real Data + fastResamplerCR, + fastResamplerSSER, + fastResamplerAVXR, + fastResamplerR, + -- * Filter firFilter, @@ -66,7 +99,8 @@ import Pipes import SDR.Util -import SDR.FilterInternal +import SDR.FilterInternal hiding (mkResampler) +import SDR.CPUID {- | A `Filter` contains all of the information needed by the `filterr` function to perform filtering. i.e. it contains the filter coefficients @@ -93,12 +127,13 @@ 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 { +data Resampler m v vm a = forall dat. 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 + startDat :: dat, + resampleOne :: dat -> Int -> v a -> vm (PrimState m) a -> m (dat, Int), + resampleCross :: dat -> Int -> v a -> v a -> vm (PrimState m) a -> m (dat, Int) } duplicate :: [a] -> [a] @@ -118,51 +153,106 @@ 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 +mkFilter :: Int + -> FilterRR + -> [Float] + -> IO (Filter IO VS.Vector VS.MVector Float) +mkFilter sizeMultiple filterFunc coeffs = do let l = length coeffs - ru = (l + 8 - 1) `quot` 8 - numCoeffsF = ru * 8 + numCoeffsF = roundUp l sizeMultiple diff = numCoeffsF - l vCoeffs = VG.fromList $ coeffs ++ replicate diff 0 evaluate vCoeffs - let filterOne = filterCAVXRR vCoeffs + let filterOne = filterFunc 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 +-- | Returns a fast Filter data structure implemented in C. For filtering real data with real coefficients. +fastFilterCR :: [Float] -- ^ The filter coefficients + -> IO (Filter IO VS.Vector VS.MVector Float) -- ^ The `Filter` data structure +fastFilterCR = mkFilter 1 filterCRR + +-- | Returns a fast Filter data structure implemented in C using SSE instructions. For filtering real data with real coefficients. +fastFilterSSER :: [Float] -- ^ The filter coefficients + -> IO (Filter IO VS.Vector VS.MVector Float) -- ^ The `Filter` data structure +fastFilterSSER = mkFilter 4 filterCSSERR + +-- | Returns a fast Filter data structure implemented in C using AVX instructions. For filtering real data with real coefficients. +fastFilterAVXR :: [Float] -- ^ The filter coefficients + -> IO (Filter IO VS.Vector VS.MVector Float) -- ^ The `Filter` data structure +fastFilterAVXR = mkFilter 8 filterCAVXRR + +-- | Returns a fast Filter data structure implemented in C using the fastest SIMD instruction set your processor supports. For filtering real data with real coefficients. +fastFilterR :: CPUInfo -- ^ The CPU's capabilities + -> [Float] -- ^ The filter coefficients + -> IO (Filter IO VS.Vector VS.MVector Float) -- ^ The `Filter` data structure +fastFilterR info = featureSelect info fastFilterCR [(hasAVX, fastFilterAVXR), (hasSSE42, fastFilterSSER)] + +mkFilterC :: Int + -> FilterRC + -> [Float] + -> IO (Filter IO VS.Vector VS.MVector (Complex Float)) +mkFilterC sizeMultiple filterFunc coeffs = do let l = length coeffs - ru = (l + 4 - 1) `quot` 4 - numCoeffsF = ru * 4 + numCoeffsF = roundUp sizeMultiple l diff = numCoeffsF - l vCoeffs = VG.fromList $ duplicate $ coeffs ++ replicate diff 0 vCoeffs2 = VG.fromList $ coeffs ++ replicate diff 0 evaluate vCoeffs - let filterOne = filterCAVXRC vCoeffs + let filterOne = filterFunc 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 +-- | Returns a fast Filter data structure implemented in C For filtering complex data with real coefficients. +fastFilterCC :: [Float] -- ^ The filter coefficients + -> IO (Filter IO VS.Vector VS.MVector (Complex Float)) -- ^ The `Filter` data structure +fastFilterCC = mkFilterC 1 filterCRC + +-- | Returns a fast Filter data structure implemented in C using SSE instructions. For filtering complex data with real coefficients. +fastFilterSSEC :: [Float] -- ^ The filter coefficients + -> IO (Filter IO VS.Vector VS.MVector (Complex Float)) -- ^ The `Filter` data structure +fastFilterSSEC = mkFilterC 2 filterCSSERC + +-- | Returns a fast Filter data structure implemented in C using AVX instructions. For filtering complex data with real coefficients. +fastFilterAVXC :: [Float] -- ^ The filter coefficients + -> IO (Filter IO VS.Vector VS.MVector (Complex Float)) -- ^ The `Filter` data structure +fastFilterAVXC = mkFilterC 4 filterCAVXRC + +-- | Returns a fast Filter data structure implemented in C using the fastest SIMD instruction set your processor supports. For filtering complex data with real coefficients. +fastFilterC :: CPUInfo -- ^ The CPU's capabilities + -> [Float] -- ^ The filter coefficients + -> IO (Filter IO VS.Vector VS.MVector (Complex Float)) -- ^ The `Filter` data structure +fastFilterC info = featureSelect info fastFilterCC [(hasAVX, fastFilterAVXC), (hasSSE42, fastFilterSSEC)] + +mkFilterSymR :: FilterRR + -> [Float] + -> IO (Filter IO VS.Vector VS.MVector Float) +mkFilterSymR filterFunc 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 + let filterOne = filterFunc vCoeffs + filterCross = filterCrossHighLevel vCoeffs2 numCoeffsF = length coeffs * 2 return $ Filter {..} +-- | Returns a fast Filter data structure implemented in C using SSE 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. +fastFilterSymSSER :: [Float] -- ^ The first half of the filter coefficients + -> IO (Filter IO VS.Vector VS.MVector Float) -- ^ The `Filter` data structure +fastFilterSymSSER = mkFilterSymR filterCSSESymmetricRR + +-- | 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. +fastFilterSymAVXR :: [Float] -- ^ The first half of the filter coefficients + -> IO (Filter IO VS.Vector VS.MVector Float) -- ^ The `Filter` data structure +fastFilterSymAVXR = mkFilterSymR filterCAVXSymmetricRR + +-- | Returns a fast Filter data structure implemented in C using the fastest SIMD instruction set your processor supports. For filtering complex data with real coefficients. For filters with symmetric coefficients, i.e. 'linear phase'. Coefficient length must be a multiple of 4. +fastFilterSymR :: CPUInfo -- ^ The CPU's capabilities + -> [Float] -- ^ The filter coefficients + -> IO (Filter IO VS.Vector VS.MVector Float) -- ^ The `Filter` data structure +fastFilterSymR info = featureSelect info (error "At least SSE4.2 required") [(hasAVX, fastFilterSymAVXR), (hasSSE42, fastFilterSymSSER)] + {-# 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) @@ -177,54 +267,120 @@ 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 +mkDecimator :: Int + -> DecimateRR + -> Int + -> [Float] + -> IO (Decimator IO VS.Vector VS.MVector Float) +mkDecimator sizeMultiple filterFunc decimationD coeffs = do let l = length coeffs - ru = (l + 8 - 1) `quot` 8 - numCoeffsD = ru * 8 + numCoeffsD = roundUp l sizeMultiple diff = numCoeffsD - l vCoeffs = VG.fromList $ coeffs ++ replicate diff 0 evaluate vCoeffs - let decimateOne = decimateCAVXRR decimationD vCoeffs + let decimateOne = filterFunc 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 +-- | Returns a fast Decimator data structure implemented in C. For decimating real data with real coefficients. +fastDecimatorCR :: Int -- ^ The decimation factor + -> [Float] -- ^ The filter coefficients + -> IO (Decimator IO VS.Vector VS.MVector Float) -- ^ The `Decimator` data structure +fastDecimatorCR = mkDecimator 1 decimateCRR + +-- | Returns a fast Decimator data structure implemented in C using SSE instructions. For decimating real data with real coefficients. +fastDecimatorSSER :: Int -- ^ The decimation factor + -> [Float] -- ^ The filter coefficients + -> IO (Decimator IO VS.Vector VS.MVector Float) -- ^ The `Decimator` data structure +fastDecimatorSSER = mkDecimator 4 decimateCSSERR + +-- | Returns a fast Decimator data structure implemented in C using AVX instructions. For decimating real data with real coefficients. +fastDecimatorAVXR :: Int -- ^ The decimation factor + -> [Float] -- ^ The filter coefficients + -> IO (Decimator IO VS.Vector VS.MVector Float) -- ^ The `Decimator` data structure +fastDecimatorAVXR = mkDecimator 8 decimateCAVXRR + +-- | Returns a fast Decimator data structure implemented in C using the fastest SIMD instruction set your processor supports. For decimating real data with real coefficients. +fastDecimatorR :: CPUInfo -- ^ The CPU's capabilities + -> Int -- ^ The decimation factor + -> [Float] -- ^ The filter coefficients + -> IO (Decimator IO VS.Vector VS.MVector Float) -- ^ The `Decimator` data structure +fastDecimatorR info = featureSelect info fastDecimatorCR [(hasAVX, fastDecimatorAVXR), (hasSSE42, fastDecimatorSSER)] + +mkDecimatorC :: Int + -> DecimateRC + -> Int + -> [Float] + -> IO (Decimator IO VS.Vector VS.MVector (Complex Float)) +mkDecimatorC sizeMultiple filterFunc decimationD coeffs = do let l = length coeffs - ru = (l + 4 - 1) `quot` 4 - numCoeffsD = ru * 4 + numCoeffsD = roundUp l sizeMultiple 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 + let decimateOne = filterFunc 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 +-- | Returns a fast Decimator data structure implemented in C. For decimating complex data with real coefficients. +fastDecimatorCC :: Int -- ^ The decimation factor + -> [Float] -- ^ The filter coefficients + -> IO (Decimator IO VS.Vector VS.MVector (Complex Float)) -- ^ The `Decimator` data structure +fastDecimatorCC = mkDecimatorC 1 decimateCRC + +-- | Returns a fast Decimator data structure implemented in C using SSE instructions. For decimating complex data with real coefficients. +fastDecimatorSSEC :: Int -- ^ The decimation factor + -> [Float] -- ^ The filter coefficients + -> IO (Decimator IO VS.Vector VS.MVector (Complex Float)) -- ^ The `Decimator` data structure +fastDecimatorSSEC = mkDecimatorC 2 decimateCSSERC + +-- | Returns a fast Decimator data structure implemented in C using AVX instructions. For decimating complex data with real coefficients. +fastDecimatorAVXC :: Int -- ^ The decimation factor + -> [Float] -- ^ The filter coefficients + -> IO (Decimator IO VS.Vector VS.MVector (Complex Float)) -- ^ The `Decimator` data structure +fastDecimatorAVXC = mkDecimatorC 4 decimateCAVXRC + +-- | Returns a fast Decimator data structure implemented in C using the fastest SIMD instruction set your processor supports. For decimating complex data with real coefficients. +fastDecimatorC :: CPUInfo -- ^ The CPU's capabilities + -> Int -- ^ The decimation factor + -> [Float] -- ^ The filter coefficients + -> IO (Decimator IO VS.Vector VS.MVector (Complex Float)) -- ^ The `Decimator` data structure +fastDecimatorC info = featureSelect info fastDecimatorCC [(hasAVX, fastDecimatorAVXC), (hasSSE42, fastDecimatorSSEC)] + +mkDecimatorSymR :: DecimateRR + -> Int + -> [Float] + -> IO (Decimator IO VS.Vector VS.MVector Float) +mkDecimatorSymR filterFunc 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 + let decimateOne = filterFunc decimationD vCoeffs + decimateCross = decimateCrossHighLevel decimationD vCoeffs2 numCoeffsD = length coeffs * 2 return $ Decimator {..} + +-- | Returns a fast Decimator data structure implemented in C using SSE 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. +fastDecimatorSymSSER :: Int -- ^ The decimation factor + -> [Float] -- ^ The first half of the filter coefficients + -> IO (Decimator IO VS.Vector VS.MVector Float) -- ^ The `Decimator` data structure +fastDecimatorSymSSER = mkDecimatorSymR decimateCSSESymmetricRR +-- | 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. +fastDecimatorSymAVXR :: Int -- ^ The decimation factor + -> [Float] -- ^ The first half of the filter coefficients + -> IO (Decimator IO VS.Vector VS.MVector Float) -- ^ The `Decimator` data structure +fastDecimatorSymAVXR = mkDecimatorSymR decimateCAVXSymmetricRR + +-- | Returns a fast Decimator data structure implemented in C using the fastest SIMD instruction set your processor supports. For decimating real data with real coefficients. For decimators with symmetric coefficients, i.e. 'linear phase'. Coefficient length must be a multiple of 4. +fastDecimatorSymR :: CPUInfo -- ^ The CPU's capabilities + -> Int -- ^ The decimation factor + -> [Float] -- ^ The filter coefficients + -> IO (Decimator IO VS.Vector VS.MVector Float) -- ^ The `Decimator` data structure +fastDecimatorSymR info = featureSelect info (error "at least AVX required") [(hasAVX, fastDecimatorSymAVXR), (hasSSE42, fastDecimatorSymSSER)] + {-# 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) @@ -235,28 +391,61 @@ 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 + let resampleOne v w x y = func <$> resampleHighLevel interpolationR decimationR vCoeffs v w x y + resampleCross v w x y z = func <$> resampleCrossHighLevel interpolationR decimationR vCoeffs v w x y z + numCoeffsR = length coeffs + func x = (x, x) + startDat = 0 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 +mkResampler :: Int + -> ResampleRR + -> Int + -> Int + -> [Float] + -> IO (Resampler IO VS.Vector VS.MVector Float) +mkResampler sizeMultiple filterFunc 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 + resamp <- filterFunc interpolationR decimationR coeffs + let resampleOne v w x y = func1 <$> resamp (fst v) w x y + resampleCross (group, offset) count x y z = do + offset' <- resampleCrossHighLevel interpolationR decimationR vCoeffs offset count x y z + return (((group + count) `mod` interpolationR, offset'), offset') + numCoeffsR = roundUp (length coeffs) (interpolationR * sizeMultiple) + func1 group = let offset = interpolationR - 1 - ((interpolationR + group * decimationR - 1) `mod` interpolationR) in ((group, offset), offset) + startDat = (0, 0) return $ Resampler {..} --} +-- | Returns a fast Resampler data structure implemented in C. For filtering real data with real coefficients. +fastResamplerCR :: Int -- ^ The interpolation factor + -> Int -- ^ The decimation factor + -> [Float] -- ^ The filter coefficients + -> IO (Resampler IO VS.Vector VS.MVector Float) -- ^ The `Resampler` data structure +fastResamplerCR = mkResampler 1 resampleCRR2 + +-- | Returns a fast Resampler data structure implemented in C using SSE instructions. For filtering real data with real coefficients. +fastResamplerSSER :: Int -- ^ The interpolation factor + -> Int -- ^ The decimation factor + -> [Float] -- ^ The filter coefficients + -> IO (Resampler IO VS.Vector VS.MVector Float) -- ^ The `Resampler` data structure +fastResamplerSSER = mkResampler 4 resampleCSSERR + +-- | Returns a fast Resampler data structure implemented in C using AVX instructions. For filtering real data with real coefficients. +fastResamplerAVXR :: Int -- ^ The interpolation factor + -> Int -- ^ The decimation factor + -> [Float] -- ^ The filter coefficients + -> IO (Resampler IO VS.Vector VS.MVector Float) -- ^ The `Resampler` data structure +fastResamplerAVXR = mkResampler 8 resampleCAVXRR + +-- | Returns a fast Resampler data structure implemented in C using the fastest SIMD instruction set your processor supports. For resampling real data with real coefficients. +fastResamplerR :: CPUInfo -- ^ The CPU's capabilities + -> Int -- ^ The interpolation factor + -> Int -- ^ The decimation factor + -> [Float] -- ^ The filter coefficients + -> IO (Resampler IO VS.Vector VS.MVector Float) -- ^ The `Resampler` data structure +fastResamplerR info = featureSelect info fastResamplerCR [(hasAVX, fastResamplerAVXR), (hasSSE42, fastResamplerSSER)] + data Buffer v a = Buffer { buffer :: v a, offset :: Int @@ -439,16 +628,16 @@ firResampler Resampler{..} blockSizeOut = do inBuf <- await outBuf <- lift $ newBuffer blockSizeOut - simple inBuf outBuf 0 + simple inBuf outBuf startDat 0 where - simple bufIn bufferOut@(Buffer bufOut offsetOut) filterOffset = do + simple bufIn bufferOut@(Buffer bufOut offsetOut) dat 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) + (dat, endOffset) <- lift $ resampleOne dat 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 @@ -457,30 +646,30 @@ bufIn' = VG.drop usedInput bufIn case VG.length bufIn' * interpolationR < numCoeffsR - endOffset of - False -> simple bufIn' bufferOut' endOffset + False -> simple bufIn' bufferOut' dat 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 + True -> simple next bufferOut' dat endOffset + False -> crossover bufIn' next bufferOut' dat endOffset - crossover bufLast bufNext bufferOut@(Buffer bufOut offsetOut) filterOffset = do + crossover bufLast bufNext bufferOut@(Buffer bufOut offsetOut) dat 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) + (dat, endOffset) <- lift $ resampleCross dat 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 + True -> simple (VG.drop (inputUsed - VG.length bufLast) bufNext) bufferOut' dat endOffset + False -> crossover (VG.drop inputUsed bufLast) bufNext bufferOut' dat endOffset -- | A DC blocking filter dcBlockingFilter :: Pipe (VS.Vector Float) (VS.Vector Float) IO () diff --git a/hs_sources/SDR/FilterInternal.hs b/hs_sources/SDR/FilterInternal.hs --- a/hs_sources/SDR/FilterInternal.hs +++ b/hs_sources/SDR/FilterInternal.hs @@ -332,26 +332,28 @@ 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 + return $ \offset num -> resampleFFIR $ func (fromIntegral num) (fromIntegral numCoeffs) (fromIntegral offset) (fromIntegral numGroups) incrementsP groupsPP where Coeffs {..} = prepareCoeffs n interpolation decimation coeffs +type ResampleRR = Int -> Int -> [Float] -> IO (Int -> Int -> VS.Vector Float -> VS.MVector RealWorld Float -> IO Int) + 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 :: ResampleRR 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 :: ResampleRR 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 :: ResampleRR resampleCAVXRR = mkResampler resampleAVXRR_c 8 {- diff --git a/hs_sources/SDR/Plot.hs b/hs_sources/SDR/Plot.hs --- a/hs_sources/SDR/Plot.hs +++ b/hs_sources/SDR/Plot.hs @@ -148,7 +148,7 @@ 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. +-- | Create a Cairo `Render` monad that draws a set of axes with the X axis centered on a specified value. centeredAxes :: Int -- ^ Image width -> Int -- ^ Image height -> Double -- ^ Center X value diff --git a/hs_sources/SDR/Util.hs b/hs_sources/SDR/Util.hs --- a/hs_sources/SDR/Util.hs +++ b/hs_sources/SDR/Util.hs @@ -11,11 +11,13 @@ convertC, convertCSSE, convertCAVX, + convertFast, -- * Scaling scaleC, scaleCSSE, scaleCAVX, + scaleFast, -- * Misc Utils cplxMap, @@ -34,6 +36,8 @@ import System.IO.Unsafe import Foreign.Storable.Complex +import SDR.CPUID + -- | A class for things that can be multiplied by a scalar. class Mult a b where mult :: a -> b -> a @@ -92,6 +96,10 @@ convertCAVX_c (fromIntegral $ VG.length inBuf) iPtr oPtr VG.freeze outBuf +-- | Create a vector of complex float samples from a vector of interleaved I Q component bytes. Uses the fastest SIMD instruction set your processor supports. +convertFast :: CPUInfo -> VS.Vector CUChar -> VS.Vector (Complex Float) +convertFast info = featureSelect info convertC [(hasAVX2, convertCAVX), (hasSSE42, convertCSSE)] + -- | Scaling foreign import ccall unsafe "scale" scale_c :: CInt -> CFloat -> Ptr CFloat -> Ptr CFloat -> IO () @@ -131,6 +139,10 @@ VS.unsafeWith (unsafeCoerce inBuf) $ \iPtr -> VS.unsafeWith (unsafeCoerce outBuf) $ \oPtr -> scaleAVX_c (fromIntegral (VG.length inBuf)) (unsafeCoerce factor) iPtr oPtr + +-- | Scale a vector. Uses the fastest SIMD instruction set your processor supports. +scaleFast :: CPUInfo -> Float -> VS.Vector Float -> VS.MVector RealWorld Float -> IO () +scaleFast info = featureSelect info scaleC [(hasAVX, scaleCAVX), (hasSSE42, scaleCSSE)] -- | Apply a function to both parts of a complex number cplxMap :: (a -> b) -- ^ The function diff --git a/sdr.cabal b/sdr.cabal --- a/sdr.cabal +++ b/sdr.cabal @@ -1,5 +1,5 @@ name: sdr -version: 0.1.0.1 +version: 0.1.0.2 synopsis: A software defined radio library description: Write software defined radio applications in Haskell. @@ -39,7 +39,7 @@ homepage: https://github.com/adamwalker/sdr bug-reports: https://github.com/adamwalker/sdr/issues build-type: Simple --- extra-source-files: +extra-source-files: Readme.md cabal-version: >=1.10 source-repository head @@ -60,7 +60,8 @@ SDR.PipeUtils, SDR.VectorUtils, SDR.ArgUtils, - SDR.FilterDesign + SDR.FilterDesign, + SDR.CPUID -- other-modules: other-extensions: ScopedTypeVariables, GADTs build-depends: @@ -101,7 +102,8 @@ c_sources/decimate.c, c_sources/convert.c, c_sources/resample.c, - c_sources/scale.c + c_sources/scale.c, + c_sources/cpuid.c hs-source-dirs: hs_sources cc-options: -mavx2 -msse4 -g @@ -112,7 +114,7 @@ base >=4.6 && <4.9, QuickCheck >=2.8 && <2.9, vector >=0.10 && <0.11, - sdr ==0.1.0.0, + sdr ==0.1.0.2, primitive >=0.5 && <0.7, storable-complex >=0.2 && <0.3, test-framework >=0.8 && <0.9, @@ -128,7 +130,7 @@ base >=4.6 && <4.9, criterion >=1.0 && <1.2, vector >=0.10 && <0.11, - sdr ==0.1.0.0, + sdr ==0.1.0.2, primitive >=0.5 && <0.7, storable-complex >=0.2 && <0.3 hs-source-dirs: benchmarks diff --git a/tests/TestSuite.hs b/tests/TestSuite.hs --- a/tests/TestSuite.hs +++ b/tests/TestSuite.hs @@ -2,6 +2,7 @@ import Control.Monad.Primitive import Data.Complex +import Control.Monad import qualified Data.Vector.Generic as VG import qualified Data.Vector.Generic.Mutable as VGM @@ -15,8 +16,20 @@ import SDR.FilterInternal import SDR.Util +import SDR.CPUID -tests = [ +sameResultM :: Monad m => (a -> a -> Bool) -> [m a] -> m Bool +sameResultM _ [] = return True +sameResultM eq (x:xs) = do + res <- x + ress <- sequence xs + return $ and $ map (eq res) ress + +sameResult :: (a -> a -> Bool) -> [a] -> Bool +sameResult _ [] = True +sameResult eq (x:xs) = and $ map (eq x) xs + +tests info = [ testGroup "filters" [ testProperty "real" propFiltersReal, testProperty "complex" propFiltersComplex @@ -32,16 +45,20 @@ testProperty "scaling" propScaleReal ] where + hasFeatures :: [(CPUInfo -> Bool, a)] -> [a] + hasFeatures = map snd . filter (($ info) . fst) + 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] + factors' = [1, 2, 3, 5, 7, 11, 13, 17, 23] + factors = elements factors' - 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 + 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 @@ -50,23 +67,25 @@ 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] + res <- run $ sameResultM eqDelta $ map (getResult num $) $ hasFeatures [ + (const True, filterHighLevel vCoeffs num vInput), + (const True, filterImperative1 vCoeffs num vInput), + (const True, filterImperative2 vCoeffs num vInput), + (const True, filterCRR vCoeffs num vInput), + (hasSSE42, filterCSSERR vCoeffs num vInput), + (hasAVX, filterCAVXRR vCoeffs num vInput), + (hasSSE42, filterCSSESymmetricRR vCoeffsHalf num vInput), + (hasAVX, filterCAVXSymmetricRR vCoeffsHalf num vInput) + ] + assert res - 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 + 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 @@ -76,22 +95,24 @@ 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] + res <- run $ sameResultM eqDeltaC $ map (getResult num $) $ hasFeatures [ + (const True, filterHighLevel vCoeffs num vInput), + (const True, filterCRC vCoeffs num vInput), + (hasSSE42, filterCSSERC vCoeffs2 num vInput), + (hasSSE42, filterCSSERC2 vCoeffs num vInput), + (hasAVX, filterCAVXRC vCoeffs2 num vInput), + (hasSSE42, filterCSSESymmetricRC vCoeffsHalf num vInput), + (hasAVX, filterCAVXSymmetricRC vCoeffsHalf num vInput) + ] + assert res - 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 + 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 @@ -100,22 +121,24 @@ 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] + res <- run $ sameResultM eqDelta $ map (getResult num $) $ hasFeatures [ + (const True, decimateHighLevel factor vCoeffs num vInput), + (const True, decimateCRR factor vCoeffs num vInput), + (hasSSE42, decimateCSSERR factor vCoeffs num vInput), + (hasAVX, decimateCAVXRR factor vCoeffs num vInput), + (hasSSE42, decimateCSSESymmetricRR factor vCoeffsHalf num vInput), + (hasAVX, decimateCAVXSymmetricRR factor vCoeffsHalf num vInput) + ] + assert res - 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 + 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 @@ -125,80 +148,88 @@ 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] + res <- run $ sameResultM eqDeltaC $ map (getResult num $) $ hasFeatures [ + (const True, decimateHighLevel factor vCoeffs num vInput), + (const True, decimateCRC factor vCoeffs num vInput), + (hasSSE42, decimateCSSERC factor vCoeffs2 num vInput), + (hasSSE42, decimateCSSERC2 factor vCoeffs num vInput), + (hasAVX, decimateCAVXRC factor vCoeffs2 num vInput), + (hasSSE42, decimateCSSESymmetricRC factor vCoeffsHalf num vInput), + (hasAVX, decimateCAVXRC2 factor vCoeffs num vInput), + (hasAVX, decimateCAVXSymmetricRC factor vCoeffsHalf num vInput) + ] + assert res - 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 + propResamplingReal = + forAll sizes $ \size -> + forAll (vectorOf size (choose (-10, 10))) $ \inBuf -> + forAll (elements [32 .. 512]) $ \numCoeffs -> + forAll (vectorOf numCoeffs (choose (-10, 10))) $ \coeffs -> + forAll (elements $ tail factors') $ \decimation -> + forAll (elements $ filter (< decimation) factors') $ \interpolation -> + forAll (elements [0..interpolation - 1]) $ \group -> + testResamplingReal size group numCoeffs interpolation decimation coeffs inBuf - testResamplingReal :: Int -> Int -> Int -> Int -> [Float] -> [Float] -> Property - testResamplingReal size numCoeffs interpolation decimation coeffs inBuf = monadicIO $ do + testResamplingReal :: Int -> Int -> Int -> Int -> Int -> [Float] -> [Float] -> Property + testResamplingReal size group 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 + offset = interpolation - 1 - ((interpolation + group * decimation - 1) `mod` interpolation) 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 + res <- run $ sameResultM eqDelta $ map (getResult num $) $ hasFeatures [ + (const True, void . resampleHighLevel interpolation decimation vCoeffs offset num vInput), + (const True, void . resampleCRR num interpolation decimation offset vCoeffs vInput), + (const True, void . resampler3 group num vInput), + (hasSSE42, void . resampler4 group num vInput), + (hasAVX, void . resampler5 group num vInput) + ] + assert res - propConversion = forAll sizes $ \size -> - forAll (vectorOf (2 * size) (choose (-10, 10))) $ \inBuf -> - testConversion size inBuf + 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] + let res = sameResult eqDeltaC $ map (VG.toList $) $ hasFeatures [ + (const True, (makeComplexBufferVect vInput :: VS.Vector (Complex Float))), + (const True, convertC vInput), + (hasSSE42, convertCSSE vInput), + (hasAVX2, convertCAVX vInput) + ] + assert res - 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 + 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] + res <- run $ sameResultM eqDelta $ map (getResult size $) $ hasFeatures [ + (const True, scaleC factor vInput), + (hasSSE42, scaleCSSE factor vInput), + (hasAVX, scaleCAVX factor vInput) + ] + assert res + 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 eqDelta x y = and $ map (uncurry eqDelta') $ zip x y where eqDelta' x y = abs (x - y) < 0.01 @@ -209,5 +240,7 @@ duplicate = concat . map func where func x = [x, x] -main = defaultMain tests +main = do + info <- getCPUInfo + defaultMain $ tests info