diff --git a/benchmarks/Benchmarks.hs b/benchmarks/Benchmarks.hs --- a/benchmarks/Benchmarks.hs +++ b/benchmarks/Benchmarks.hs @@ -71,6 +71,10 @@ resampler4 <- resampleCSSERR interpolation decimation coeffsList resampler5 <- resampleCAVXRR interpolation decimation coeffsList + resampler3C <- resampleCRC interpolation decimation coeffsList + resampler4C <- resampleCSSERC interpolation decimation coeffsList + resampler5C <- resampleCAVXRC interpolation decimation coeffsList + --Benchmarks defaultMain [ bgroup "filter" [ @@ -124,7 +128,10 @@ (hasAVX, bench "cAVX" $ nfIO $ resampler5 (num `quot` decimation) 0 inBuf outBuf) ], bgroup "complex" $ hasFeatures [ - (const True, bench "highLevel" $ nfIO $ resampleHighLevel interpolation decimation coeffs 0 (num `quot` decimation) inBufComplex outBufComplex) + (const True, bench "highLevel" $ nfIO $ resampleHighLevel interpolation decimation coeffs 0 (num `quot` decimation) inBufComplex outBufComplex), + (const True, bench "c" $ nfIO $ resampler3C (num `quot` decimation) 0 inBufComplex outBufComplex), + (const True, bench "SSE" $ nfIO $ resampler4C (num `quot` decimation) 0 inBufComplex outBufComplex), + (const True, bench "AVX" $ nfIO $ resampler5C (num `quot` decimation) 0 inBufComplex outBufComplex) ] ], bgroup "scaling" $ hasFeatures [ diff --git a/c_sources/convert.c b/c_sources/convert.c --- a/c_sources/convert.c +++ b/c_sources/convert.c @@ -84,3 +84,19 @@ } } +void convertBladeRFTransmit(int num, float *in, int16_t *out){ + int i; + for(i=0; i 2047) res = 2047; + if(res < -2048) res = -2048; + + out[i] = res; + } +} + diff --git a/c_sources/resample.c b/c_sources/resample.c --- a/c_sources/resample.c +++ b/c_sources/resample.c @@ -13,7 +13,7 @@ /* * 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){ +void resampleRR(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 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 + -> VFSM.Stream m (Complex a) -- ^ The input stream + -> VFSM.Stream m a -- ^ The output stream fmDemodStr = mapAccumMV func where {-# INLINE func #-} @@ -27,15 +29,15 @@ -- | FM demodulate a vector of complex samples {-# INLINE fmDemodVec #-} -fmDemodVec :: (RealFloat a, Vector v (Complex a), Vector v a) +fmDemodVec :: (RealFloat a, VG.Vector v (Complex a), VG.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 +fmDemodVec init inp = VG.unstream $ flip VFBM.fromStream (VFBS.Exact $ VG.length inp) $ fmDemodStr init $ VFBM.elements $ VG.stream inp -- | 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 :: (RealFloat a, VG.Vector v (Complex a), VG.Vector v a) => Pipe (v (Complex a)) (v a) IO () fmDemod = func 0 where func lastSample = do diff --git a/hs_sources/SDR/FFT.hs b/hs_sources/SDR/FFT.hs --- a/hs_sources/SDR/FFT.hs +++ b/hs_sources/SDR/FFT.hs @@ -3,7 +3,6 @@ {-| Fast FFTs using FFTW -} module SDR.FFT ( -- * Windows - fftFixup, hamming, hanning, blackman, @@ -23,32 +22,21 @@ 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 qualified Data.Vector.Generic as VG +import qualified Data.Vector.Storable as VS +import qualified Data.Vector.Storable.Mutable as VSM import Pipes import Numeric.FFTW import SDR.FilterDesign +import SDR.VectorUtils 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 @@ -66,7 +54,7 @@ ina <- lift $ mallocForeignBufferAligned samples let inv = VSM.unsafeFromForeignPtr0 ina samples - lift $ VGM.fill inv $ VFS.mapM return $ VG.stream inv' + lift $ copyInto inv inv' let (fp, offset, length) = VSM.unsafeToForeignPtr inv @@ -94,7 +82,7 @@ ina <- lift $ mallocForeignBufferAligned samples let inv = VSM.unsafeFromForeignPtr0 ina samples - lift $ VGM.fill inv $ VFS.mapM return $ VG.stream inv' + lift $ copyInto inv inv' let (fp, offset, length) = VSM.unsafeToForeignPtr inv lift $ withForeignPtr fp $ \fpp -> @@ -132,7 +120,7 @@ ina <- mallocForeignBufferAligned samples let inv = VSM.unsafeFromForeignPtr0 ina samples - VGM.fill inv $ VFS.mapM return $ VG.stream res + copyInto inv res let (fp, offset, length) = VSM.unsafeToForeignPtr inv 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 @@ -76,6 +76,12 @@ fastResamplerAVXR, fastResamplerR, + -- *** Complex Data + fastResamplerCC, + fastResamplerSSEC, + fastResamplerAVXC, + fastResamplerC, + -- * Filter firFilter, @@ -100,7 +106,7 @@ import Pipes import SDR.Util -import SDR.FilterInternal hiding (mkResampler) +import SDR.FilterInternal hiding (mkResampler, mkResamplerC) import SDR.CPUID {- | A `Filter` contains all of the information needed by the `filterr` @@ -418,6 +424,25 @@ startDat = (0, 0) return $ Resampler {..} +mkResamplerC :: Int + -> ResampleRC + -> Int + -> Int + -> [Float] + -> IO (Resampler IO VS.Vector VS.MVector (Complex Float)) +mkResamplerC sizeMultiple filterFunc interpolationR decimationR coeffs = do + let vCoeffs = VG.fromList coeffs + evaluate vCoeffs + 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 @@ -446,6 +471,35 @@ -> [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)] + +-- | Returns a fast Resampler data structure implemented in C. For filtering complex data with real coefficients. +fastResamplerCC :: Int -- ^ The interpolation factor + -> Int -- ^ The decimation factor + -> [Float] -- ^ The filter coefficients + -> IO (Resampler IO VS.Vector VS.MVector (Complex Float)) -- ^ The `Resampler` data structure +fastResamplerCC = mkResamplerC 1 resampleCRC + +-- | Returns a fast Resampler data structure implemented in C using SSE instructions. For filtering complex data with real coefficients. +fastResamplerSSEC :: Int -- ^ The interpolation factor + -> Int -- ^ The decimation factor + -> [Float] -- ^ The filter coefficients + -> IO (Resampler IO VS.Vector VS.MVector (Complex Float)) -- ^ The `Resampler` data structure +fastResamplerSSEC = mkResamplerC 4 resampleCSSERC + +-- | Returns a fast Resampler data structure implemented in C using AVX instructions. For filtering complex data with real coefficients. +fastResamplerAVXC :: Int -- ^ The interpolation factor + -> Int -- ^ The decimation factor + -> [Float] -- ^ The filter coefficients + -> IO (Resampler IO VS.Vector VS.MVector (Complex Float)) -- ^ The `Resampler` data structure +fastResamplerAVXC = mkResamplerC 8 resampleCAVXRC + +-- | Returns a fast Resampler data structure implemented in C using the fastest SIMD instruction set your processor supports. For resampling complex data with real coefficients. +fastResamplerC :: CPUInfo -- ^ The CPU's capabilities + -> Int -- ^ The interpolation factor + -> Int -- ^ The decimation factor + -> [Float] -- ^ The filter coefficients + -> IO (Resampler IO VS.Vector VS.MVector (Complex Float)) -- ^ The `Resampler` data structure +fastResamplerC info = featureSelect info fastResamplerCC [(hasAVX, fastResamplerAVXC), (hasSSE42, fastResamplerSSEC)] data Buffer v a = Buffer { buffer :: v a, diff --git a/hs_sources/SDR/FilterDesign.hs b/hs_sources/SDR/FilterDesign.hs --- a/hs_sources/SDR/FilterDesign.hs +++ b/hs_sources/SDR/FilterDesign.hs @@ -80,7 +80,11 @@ --ts is really the ratio ts / sampling period -- | Square root raised cosine -srrc :: (Ord a, Floating a) => Int -> Int -> a -> [a] +srrc :: (Ord a, Floating a) + => Int -- ^ size: from [-n .. n] + -> Int -- ^ sampling period + -> a -- ^ beta + -> [a] srrc n ts beta = map func [(-n) .. n] where func x 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 @@ -20,15 +20,14 @@ 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 qualified Data.Vector.Fusion.Bundle as VFB 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 +filterHighLevel coeffs num inBuf outBuf = fill (VFB.generate num dotProd) outBuf where dotProd offset = VG.sum $ VG.zipWith mult (VG.unsafeDrop offset inBuf) coeffs @@ -156,7 +155,7 @@ 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) + x = VFB.map dotProd (VFB.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 () @@ -265,7 +264,7 @@ | 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" +foreign import ccall unsafe "resampleRR" 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 () @@ -325,6 +324,12 @@ VS.unsafeWith (unsafeCoerce outBuf) $ \oPtr -> func iPtr oPtr +resampleFFIC :: (Ptr CFloat -> Ptr CFloat -> IO CInt) -> VS.Vector (Complex Float) -> VSM.MVector RealWorld (Complex Float) -> IO Int +resampleFFIC 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) @@ -338,7 +343,7 @@ type ResampleRR = Int -> Int -> [Float] -> IO (Int -> Int -> VS.Vector Float -> VS.MVector RealWorld Float -> IO Int) -foreign import ccall unsafe "resample2" +foreign import ccall unsafe "resample2RR" resample2_c :: CInt -> CInt -> CInt -> CInt -> Ptr CInt -> Ptr (Ptr CFloat) -> Ptr CFloat -> Ptr CFloat -> IO CInt resampleCRR2 :: ResampleRR @@ -356,6 +361,35 @@ resampleCAVXRR :: ResampleRR resampleCAVXRR = mkResampler resampleAVXRR_c 8 +type ResampleRC = Int -> Int -> [Float] -> IO (Int -> Int -> VS.Vector (Complex Float) -> VS.MVector RealWorld (Complex Float) -> IO Int) + +mkResamplerC :: ResampleR -> Int -> Int -> Int -> [Float] -> IO (Int -> Int -> VS.Vector (Complex Float) -> VS.MVector RealWorld (Complex Float) -> IO Int) +mkResamplerC func n interpolation decimation coeffs = do + groupsP <- mapM newArray $ map (map realToFrac) groups + groupsPP <- newArray groupsP + incrementsP <- newArray $ map fromIntegral increments + return $ \offset num -> resampleFFIC $ func (fromIntegral num) (fromIntegral numCoeffs) (fromIntegral offset) (fromIntegral numGroups) incrementsP groupsPP + where + Coeffs {..} = prepareCoeffs n interpolation decimation coeffs + +foreign import ccall unsafe "resample2RC" + resample2RC_c :: CInt -> CInt -> CInt -> CInt -> Ptr CInt -> Ptr (Ptr CFloat) -> Ptr CFloat -> Ptr CFloat -> IO CInt + +resampleCRC :: ResampleRC +resampleCRC = mkResamplerC resample2RC_c 1 + +foreign import ccall unsafe "resampleSSERC" + resampleCSSERC_c :: CInt -> CInt -> CInt -> CInt -> Ptr CInt -> Ptr (Ptr CFloat) -> Ptr CFloat -> Ptr CFloat -> IO CInt + +resampleCSSERC :: ResampleRC +resampleCSSERC = mkResamplerC resampleCSSERC_c 4 + +foreign import ccall unsafe "resampleAVXRC" + resampleAVXRC_c :: CInt -> CInt -> CInt -> CInt -> Ptr CInt -> Ptr (Ptr CFloat) -> Ptr CFloat -> Ptr CFloat -> IO CInt + +resampleCAVXRC :: ResampleRC +resampleCAVXRC = mkResamplerC resampleAVXRC_c 8 + {- - Cross buffer -} @@ -364,12 +398,12 @@ 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) + x = VFB.map dotProd (VFB.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 +filterCrossHighLevel coeffs num lastBuf nextBuf outBuf = fill (VFB.generate num dotProd) outBuf where dotProd i = VG.sum $ VG.zipWith mult (VG.unsafeDrop i lastBuf VG.++ nextBuf) coeffs diff --git a/hs_sources/SDR/PipeUtils.hs b/hs_sources/SDR/PipeUtils.hs --- a/hs_sources/SDR/PipeUtils.hs +++ b/hs_sources/SDR/PipeUtils.hs @@ -54,7 +54,11 @@ rate' (buffers + 1) rate' 1 -pMapAccum :: (Monad m) => (acc -> x -> (acc, y)) -> acc -> Pipe x y m () +-- | mapAccum for Pipes +pMapAccum :: (Monad m) + => (acc -> x -> (acc, y)) -- ^ Accumulating function + -> acc -- ^ Initial value of the accumulator + -> Pipe x y m () pMapAccum func acc = go acc where go acc = do diff --git a/hs_sources/SDR/Serialize.hs b/hs_sources/SDR/Serialize.hs --- a/hs_sources/SDR/Serialize.hs +++ b/hs_sources/SDR/Serialize.hs @@ -31,8 +31,8 @@ import Data.ByteString as BS import System.IO -import Data.Vector.Generic as VG hiding ((++)) -import Data.Vector.Storable as VS hiding ((++)) +import qualified Data.Vector.Generic as VG hiding ((++)) +import qualified Data.Vector.Storable as VS hiding ((++)) import Pipes import qualified Pipes.Prelude as P 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 @@ -6,42 +6,72 @@ Mult, mult, - -- * Conversion to Floating Point + -- * Conversion to floating point for reception + -- ** RTLSDR interleavedIQUnsigned256ToFloat, interleavedIQUnsignedByteToFloat, interleavedIQUnsignedByteToFloatSSE, interleavedIQUnsignedByteToFloatAVX, interleavedIQUnsignedByteToFloatFast, + -- ** BladeRF interleavedIQSigned2048ToFloat, interleavedIQSignedWordToFloat, interleavedIQSignedWordToFloatSSE, interleavedIQSignedWordToFloatAVX, + interleavedIQSignedWordToFloatFast, + -- * Conversion from floating point for transmission + -- ** BladeRF + complexFloatToInterleavedIQSigned2048, + complexFloatToInterleavedIQSignedWord, + -- * Scaling scaleC, scaleCSSE, scaleCAVX, scaleFast, - -- * Misc Utils + -- * Mapping over complex numbers cplxMap, - quarterBandUp + + -- * Frequency shifting + halfBandUp, + quarterBandUp, + + -- * Data streams + streamString, + streamRandom, + + -- * Automatic gain control + agc, + agcPipe ) where import Foreign.C.Types import Data.Complex -import Data.Vector.Generic as VG hiding ((++)) +import qualified Data.Vector.Generic as VG import qualified Data.Vector.Generic.Mutable as VGM -import Data.Vector.Storable as VS hiding ((++)) -import Data.Vector.Storable.Mutable as VSM +import qualified Data.Vector.Storable as VS +import qualified Data.Vector.Storable.Mutable as VSM import Control.Monad.Primitive import Unsafe.Coerce import Foreign.Ptr import System.IO.Unsafe import Foreign.Storable.Complex +import Control.Monad +import qualified System.Random.MWC as R +import Data.Bits +import Pipes +import qualified Pipes.Prelude as P +import Data.Word +import Foreign.Storable +import Control.Arrow as A +import Data.Tuple import SDR.CPUID +import SDR.VectorUtils +import SDR.PipeUtils -- | A class for things that can be multiplied by a scalar. class Mult a b where @@ -116,7 +146,7 @@ foreign import ccall unsafe "convertCBladeRF" convertCBladeRF_c :: CInt -> Ptr CShort -> Ptr CFloat -> IO () --- | Same as `interleavedIQUnsigned256ToFloat` but written in C and specialized for unsigned byte inputs and Float outputs. +-- | Same as `interleavedIQUnsigned256ToFloat` but written in C and specialized for signed short inputs and Float outputs. interleavedIQSignedWordToFloat :: VS.Vector CShort -> VS.Vector (Complex Float) interleavedIQSignedWordToFloat inBuf = unsafePerformIO $ do outBuf <- VGM.new $ VG.length inBuf `quot` 2 @@ -128,7 +158,7 @@ foreign import ccall unsafe "convertCSSEBladeRF" convertCSSEBladeRF_c :: CInt -> Ptr CShort -> Ptr CFloat -> IO () --- | Same as `interleavedIQUnsigned256ToFloat` but written in C using SSE intrinsics and specialized for unsigned byte inputs and Float outputs. +-- | Same as `interleavedIQUnsigned256ToFloat` but written in C using SSE intrinsics and specialized for signed short inputs and Float outputs. interleavedIQSignedWordToFloatSSE :: VS.Vector CShort -> VS.Vector (Complex Float) interleavedIQSignedWordToFloatSSE inBuf = unsafePerformIO $ do outBuf <- VGM.new $ VG.length inBuf `quot` 2 @@ -140,7 +170,7 @@ foreign import ccall unsafe "convertCAVXBladeRF" convertCAVXBladeRF_c :: CInt -> Ptr CShort -> Ptr CFloat -> IO () --- | Same as `interleavedIQUnsigned256ToFloat` but written in C using AVX intrinsics and specialized for unsigned byte inputs and Float outputs. +-- | Same as `interleavedIQUnsigned256ToFloat` but written in C using AVX intrinsics and specialized for signed short inputs and Float outputs. interleavedIQSignedWordToFloatAVX :: VS.Vector CShort -> VS.Vector (Complex Float) interleavedIQSignedWordToFloatAVX inBuf = unsafePerformIO $ do outBuf <- VGM.new $ VG.length inBuf `quot` 2 @@ -149,6 +179,33 @@ convertCAVXBladeRF_c (fromIntegral $ VG.length inBuf) iPtr oPtr VG.freeze outBuf +-- | Same as `interleavedIQSigned2048ToFloat` but uses the fastest SIMD instruction set your processor supports and specialized for signed short inputs and Float outputs. +interleavedIQSignedWordToFloatFast :: CPUInfo -> VS.Vector CShort -> VS.Vector (Complex Float) +interleavedIQSignedWordToFloatFast info = featureSelect info interleavedIQSignedWordToFloat [(hasAVX2, interleavedIQSignedWordToFloatAVX), (hasSSE42, interleavedIQSignedWordToFloatSSE)] + +-- | Create a vector of interleaved I Q component integral samples from a vector of complex Floats. Each input ranges from -2048 to 2047. This is the format the BladeRF uses. +complexFloatToInterleavedIQSigned2048 :: (Integral b, RealFrac a, VG.Vector v1 (Complex a), VG.Vector v2 b) => v1 (Complex a) -> v2 b +complexFloatToInterleavedIQSigned2048 input = VG.generate (VG.length input * 2) convert + where + {-# INLINE convert #-} + convert idx + | even idx = convert' $ realPart (input `VG.unsafeIndex` (idx `quot` 2)) + | odd idx = convert' $ imagPart (input `VG.unsafeIndex` (idx `quot` 2)) + {-# INLINE convert' #-} + convert' val = round $ val * 2048 + +foreign import ccall unsafe "convertBladeRFTransmit" + convertBladeRFTransmit_c :: CInt -> Ptr CFloat -> Ptr CShort -> IO () + +-- | Same as `complexFloatToInterleavedIQSigned2048` but written in C and specialized for Float inputs and signed short outputs. +complexFloatToInterleavedIQSignedWord :: VS.Vector (Complex Float) -> VS.Vector CShort +complexFloatToInterleavedIQSignedWord inBuf = unsafePerformIO $ do + outBuf <- VGM.new $ VG.length inBuf * 2 + VS.unsafeWith (unsafeCoerce inBuf) $ \iPtr -> + VSM.unsafeWith outBuf $ \oPtr -> + convertBladeRFTransmit_c (fromIntegral $ VG.length inBuf * 2) iPtr oPtr + VG.freeze outBuf + -- | Scaling foreign import ccall unsafe "scale" scale_c :: CInt -> CFloat -> Ptr CFloat -> Ptr CFloat -> IO () @@ -199,6 +256,16 @@ -> Complex b -- ^ Output complex number cplxMap f (x :+ y) = f x :+ f y +-- | Multiplication by this vector shifts all frequencies up by 1/2 of the sampling frequency +halfBandUp :: (VG.Vector v n, Num n) + => Int -- ^ The length of the Vector + -> v n +halfBandUp size = VG.generate size func + where + func idx + | even idx = 1 + | otherwise = -1 + -- | 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 @@ -212,4 +279,67 @@ | m == 3 = 0 :+ (-1) where m = idx `mod` 4 + +-- | A Producer that streams vectors of the bits that make up the string argument concatenated repeatedly. Each bit is encoded as a float with value (+1) for 1 and (-1) for 0. +streamString :: forall m b. (FiniteBits b, Monad m) + => [b] -- ^ The string whose bits are to be streamed + -> Int -- ^ The size of each streamed vector + -> Producer (VS.Vector Float) m () +streamString str size = P.unfoldr (return . Right . func) (str, 0) + where + bitsPerChar = finiteBitSize (undefined :: b) + toFloat :: Bool -> Float + toFloat x = if x then 1 else (-1) + func = vUnfoldr size funcy + where + funcy ([], offsetChar) = funcy (str, 0) + funcy (rem@(x:xs), offsetChar) + | offsetChar == bitsPerChar = funcy (xs, 0) + | otherwise = (toFloat $ testBit x offsetChar, (rem, offsetChar + 1)) + +-- | A Producer that streams vectors of random bits. Each bit is encoded as a float with value (+1) for 1 and (-1) for 0. +streamRandom :: forall m. PrimMonad m + => Int -- ^ The size of each streamed vector + -> Producer (VS.Vector Float) m () +streamRandom size = do + gen <- lift R.create + start <- lift $ R.uniform gen + P.unfoldr (liftM Right . func gen) (start, 0) + where + toFloat :: Bool -> Float + toFloat x = if x then 1 else (-1) + func :: R.Gen (PrimState m) -> (Word64, Int) -> m (VS.Vector Float, (Word64, Int)) + func gen = vUnfoldrM size funcy + where + funcy (current, offset) = do + let res = toFloat $ testBit current offset + if offset == 63 then do + current' <- R.uniform gen + return (res, (current', 0)) + else return (res, (current, offset+1)) + +(a :+ b) `cdiv` y = (a/y) :+ (b/y) +(a :+ b) `cmul` y = (a*y) :+ (b*y) + +-- | Simple automatic gain control +agc :: (Num a, Storable a, RealFloat a) + => a -- ^ a + -> a -- ^ reference + -> a -- ^ initial state + -> VS.Vector (Complex a) -- ^ input vector + -> (a, VS.Vector (Complex a)) -- ^ (final state, output vector) +agc mu reference state input = A.first snd $ swap $ vUnfoldr (VS.length input) go (0, state) + where + go (offset, state) = + let + corrected = (input VS.! offset) `cmul` state + state' = state + mu * (reference - magnitude corrected) + in (corrected, (offset + 1, state')) + +-- | Simple automatic gain control pipe +agcPipe :: (Num a, Storable a, RealFloat a, Monad m) + => a -- ^ a + -> a -- ^ reference + -> Pipe (VS.Vector (Complex a)) (VS.Vector (Complex a)) m () +agcPipe mu reference = pMapAccum (agc mu reference) 1 diff --git a/hs_sources/SDR/VectorUtils.hs b/hs_sources/SDR/VectorUtils.hs --- a/hs_sources/SDR/VectorUtils.hs +++ b/hs_sources/SDR/VectorUtils.hs @@ -5,18 +5,19 @@ mapAccumMV, stride, fill, - vUnfoldr + copyInto, + vUnfoldr, + vUnfoldrM ) where import Control.Monad import Control.Monad.Primitive import Control.Monad.ST -import Data.Vector.Generic as VG hiding ((++)) +import qualified Data.Vector.Generic as VG 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 ((++)) +import qualified Data.Vector.Fusion.Bundle as VFB +import Data.Vector.Fusion.Stream.Monadic as VFSM {-| 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 @@ -25,18 +26,18 @@ mapAccumMV :: (Monad m) => (acc -> x -> m (acc, y)) -- ^ The function -> acc -- ^ The initial accumulator - -> Stream m x -- ^ The input stream + -> VFSM.Stream m x -- ^ The input stream -> Stream m y -- ^ The output stream -mapAccumMV func z (Stream step s sz) = Stream step' (s, z) sz +mapAccumMV func z (VFSM.Stream step s) = VFSM.Stream step' (s, z) where step' (s, acc) = do r <- step s case r of - Yield y s' -> do + VFB.Yield y s' -> do (!acc', !res) <- func acc y - return $ Yield res (s', acc') - Skip s' -> return $ Skip (s', acc) - Done -> return Done + return $ VFB.Yield res (s', acc') + VFB.Skip s' -> return $ VFB.Skip (s', acc) + VFB.Done -> return VFB.Done {-| Create a vector from another vector containing only the elements that occur every stride elements in the source vector. @@ -46,25 +47,39 @@ => Int -- ^ The stride -> v a -- ^ The input Vector -> v a -- ^ The output Vector -stride str inv = VG.unstream $ VFS.unfoldr func 0 +stride str inv = VG.unstream $ VFB.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 +-- | Fill a mutable vector from a monadic stream. This appears to be missing from the Vector library. {-# INLINE fill #-} fill :: (PrimMonad m, Functor m, VGM.MVector vm a) - => VFS.MStream m a -- ^ The input Stream + => VFB.Bundle v a -- ^ The input Stream -> vm (PrimState m) a -- ^ The mutable Vector to stream into -> m () -fill str outBuf = void $ VFSM.foldM' put 0 str +fill str outBuf = void $ VFB.foldM' put 0 str where put i x = do VGM.unsafeWrite outBuf i x return $ i + 1 -vUnfoldr :: VG.Vector v x => Int -> (acc -> (x, acc)) -> acc -> (v x, acc) +-- | Copy a Vector into a mutable vector +{-# INLINE copyInto #-} +copyInto :: (PrimMonad m, VGM.MVector vm a, VG.Vector v a) + => vm (PrimState m) a -- ^ The destination + -> v a -- ^ The source + -> m () +copyInto dst src = fill (VG.stream src) dst + +-- | Similar to unfoldrN from the vector package but the generator function cannot terminate and it returns the final value of the seed in addition to the vector +{-# INLINE vUnfoldr #-} +vUnfoldr :: VG.Vector v x + => Int -- ^ Generates a vector with this size + -> (acc -> (x, acc)) -- ^ The generator function + -> acc -- ^ The initial value of the seed + -> (v x, acc) -- ^ The (vector, final value of seed) result vUnfoldr size func acc = runST $ do vect <- VGM.new size acc' <- go vect 0 acc @@ -79,3 +94,26 @@ let (res, acc') = func acc VGM.write vect offset res go' (offset + 1) acc' + +-- | The same as `vUnfoldr` but the generator function is monadic +{-# INLINE vUnfoldrM #-} +vUnfoldrM :: (PrimMonad m, VG.Vector v x) + => Int -- ^ Generates a vector with this size + -> (acc -> m (x, acc)) -- ^ The monadic generator function + -> acc -- ^ The initial value of the seed + -> m (v x, acc) -- ^ The (vector, final value of seed) result +vUnfoldrM size func acc = do + vect <- VGM.new size + acc' <- go vect 0 acc + vect' <- VG.unsafeFreeze vect + return (vect', acc') + where + go vect offset acc = go' offset acc + where + go' offset acc + | offset == size = return acc + | otherwise = do + (res, acc') <- func acc + VGM.write vect offset res + go' (offset + 1) acc' + diff --git a/sdr.cabal b/sdr.cabal --- a/sdr.cabal +++ b/sdr.cabal @@ -1,5 +1,5 @@ name: sdr -version: 0.1.0.5 +version: 0.1.0.6 synopsis: A software defined radio library description: Write software defined radio applications in Haskell. @@ -24,7 +24,7 @@ . * PulseAudio sound sink . - * based radio source supported and other sources are easily added + * and based radio sources/sinks supported and other sources are easily added . See for more features and screenshots. . @@ -80,7 +80,7 @@ pipes-bytestring >=2.0 && <2.2, dynamic-graph ==0.1.0.8, array >=0.4 && <0.6, - vector >=0.10 && <0.11, + vector >=0.11 && <0.12, tuple >=0.2 && <0.4, OpenGL >=2.11 && <2.13, GLFW-b >=1.4.7 && <1.4.8, @@ -89,11 +89,12 @@ 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, + cereal >=0.4 && <0.6, + optparse-applicative >=0.11 && <0.13, Decimal >=0.4 && <0.5, - Chart >=1.3 && <1.5, - Chart-cairo >=1.3 && <1.5 + Chart >=1.3 && <1.6, + Chart-cairo >=1.3 && <1.6, + mwc-random -- hs-source-dirs: default-language: Haskell2010 ghc-options: -O2 @@ -115,8 +116,8 @@ build-depends: base >=4.6 && <4.9, QuickCheck >=2.8 && <2.9, - vector >=0.10 && <0.11, - sdr ==0.1.0.5, + vector >=0.11 && <0.12, + sdr ==0.1.0.6, primitive >=0.5 && <0.7, storable-complex >=0.2 && <0.3, test-framework >=0.8 && <0.9, @@ -131,8 +132,8 @@ build-depends: base >=4.6 && <4.9, criterion >=1.0 && <1.2, - vector >=0.10 && <0.11, - sdr ==0.1.0.5, + vector >=0.11 && <0.12, + sdr ==0.1.0.6, 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 @@ -39,7 +39,8 @@ testProperty "complex" propDecimationComplex ], testGroup "resamplers" [ - testProperty "real" propResamplingReal + testProperty "real" propResamplingReal, + testProperty "complex" propResamplingComplex ], testGroup "conversion" [ testProperty "rtlsdr" propConversionRTLSDR, @@ -188,6 +189,37 @@ 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 + + propResamplingComplex = + forAll sizes $ \size -> + forAll (vectorOf size (choose (-10, 10))) $ \inBufR -> + forAll (vectorOf size (choose (-10, 10))) $ \inBufI -> + 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 -> + testResamplingComplex size group numCoeffs interpolation decimation coeffs $ zipWith (:+) inBufR inBufI + + testResamplingComplex :: Int -> Int -> Int -> Int -> Int -> [Float] -> [Complex Float] -> Property + testResamplingComplex 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 $ resampleCRC interpolation decimation (coeffs ++ reverse coeffs) + resampler4 <- run $ resampleCSSERC interpolation decimation (coeffs ++ reverse coeffs) + resampler5 <- run $ resampleCAVXRC interpolation decimation (coeffs ++ reverse coeffs) + + res <- run $ sameResultM eqDeltaC $ map (getResult num $) $ hasFeatures [ + (const True, void . resampleHighLevel interpolation decimation vCoeffs offset num vInput), (const True, void . resampler3 group num vInput), (hasSSE42, void . resampler4 group num vInput), (hasAVX, void . resampler5 group num vInput)