Implementation FFT (#44)

* fft

* implemented rfft backward forward (#39)

* rfft_forward

* passed real_forward

* irfft

* modified to_complex

* implemented irfft but memory leak occurred... (#39)

* passed rfft and irfft test (#39)

* fixed bug (#42)

* modified vDSP rFFT

* add fft

---------

Co-authored-by: jjjkkkjjj-mizuno <[email protected]>

.swiftpm/xcode/xcshareddata/xcschemes/Matft.xcscheme

@@ -34,6 +34,20 @@
                ReferencedContainer = "container:">
             </BuildableReference>
          </BuildActionEntry>
+         <BuildActionEntry
+            buildForTesting = "YES"
+            buildForRunning = "YES"
+            buildForProfiling = "YES"
+            buildForArchiving = "YES"
+            buildForAnalyzing = "YES">
+            <BuildableReference
+               BuildableIdentifier = "primary"
+               BlueprintIdentifier = "pocketFFT"
+               BuildableName = "pocketFFT"
+               BlueprintName = "pocketFFT"
+               ReferencedContainer = "container:">
+            </BuildableReference>
+         </BuildActionEntry>
       </BuildActionEntries>
    </BuildAction>
    <TestAction

MatftDemo/MatftDemo/accelerate.playground/Contents.swift

@@ -0,0 +1,68 @@
+import Accelerate
+
+
+var sreal: [Float] = [Float(0), 1, 0, 0]
+var simag: [Float] = Array(repeating: Float.zero, count: sreal.count)
+
+var real = Array(repeating: Float.zero, count: 4)
+var imag = Array(repeating: Float.zero, count: 4)
+let log2N = 2
+
+let setup = vDSP_create_fftsetup(vDSP_Length(log2N), FFTRadix(kFFTRadix2))!
+let direction = kFFTDirection_Forward
+
+sreal.withUnsafeMutableBufferPointer{
+    srealptr in
+    simag.withUnsafeMutableBufferPointer{
+        simagptr in
+        real.withUnsafeMutableBufferPointer{
+            realptr in
+            imag.withUnsafeMutableBufferPointer{
+                imagptr in
+                var src = DSPSplitComplex(realp: srealptr.baseAddress!, imagp: simagptr.baseAddress!)
+                var dst = DSPSplitComplex(realp: realptr.baseAddress!, imagp: imagptr.baseAddress!)
+                //vDSP_fft_zrop(setup, &src, vDSP_Stride(1), &dst, vDSP_Stride(1), vDSP_Length(log2N), FFTDirection(direction))
+                vDSP_fft_zrip(setup, &src, vDSP_Stride(1), vDSP_Length(log2N), FFTDirection(direction))
+
+            }
+        }
+    }
+
+}
+
+vDSP_destroy_fftsetup(setup)
+
+print(sreal)
+print(simag)
+
+sreal = [Float(0), 1, 0, 0]
+simag = Array(repeating: Float.zero, count: sreal.count)
+
+real = Array(repeating: Float.zero, count: 4)
+imag = Array(repeating: Float.zero, count: 4)
+
+let setup_dft = vDSP_DFT_zop_CreateSetup(nil, vDSP_Length(4), vDSP_DFT_Direction.FORWARD)
+
+sreal.withUnsafeMutableBufferPointer{
+    srealptr in
+    simag.withUnsafeMutableBufferPointer{
+        simagptr in
+        real.withUnsafeMutableBufferPointer{
+            realptr in
+            imag.withUnsafeMutableBufferPointer{
+                imagptr in
+                var src = DSPSplitComplex(realp: srealptr.baseAddress!, imagp: simagptr.baseAddress!)
+                var dst = DSPSplitComplex(realp: realptr.baseAddress!, imagp: imagptr.baseAddress!)
+                //vDSP_fft_zrop(setup, &src, vDSP_Stride(1), &dst, vDSP_Stride(1), vDSP_Length(log2N), FFTDirection(direction))
+                vDSP_DFT_Execute(setup, <#T##__Ir: UnsafePointer<Float>##UnsafePointer<Float>#>, <#T##__Ii: UnsafePointer<Float>##UnsafePointer<Float>#>, <#T##__Or: UnsafeMutablePointer<Float>##UnsafeMutablePointer<Float>#>, <#T##__Oi: UnsafeMutablePointer<Float>##UnsafeMutablePointer<Float>#>)(setup, &src, vDSP_Stride(1), vDSP_Length(log2N), FFTDirection(direction))
+
+            }
+        }
+    }
+
+}
+vDSP.FFT(log2n: <#T##vDSP_Length#>, radix: <#T##vDSP.Radix#>, ofType: <#T##_.Type#>).fo
+vDSP_destroy_fftsetup(setup)
+
+print(sreal)
+print(simag)

MatftDemo/MatftDemo/accelerate.playground/Sources/New File.swift

@@ -0,0 +1,117 @@
+import UIKit
+import Matft
+import Accelerate
+/*
+let rgb = UIImage(named: "rena.jpeg")!
+
+func convertToGrayScale(image: UIImage) -> UIImage{
+    let width = Int(image.size.width)
+    let height = Int(image.size.height)
+    let channel = Int(image.cgImage!.bitsPerPixel/8)
+
+    let imageRect: CGRect = CGRect(x:0, y:0, width: width, height: height)
+
+    let colorSpace = CGColorSpaceCreateDeviceGray()
+    let bitmapInfo = CGBitmapInfo(rawValue: CGImageAlphaInfo.none.rawValue | CGImageByteOrderInfo.orderDefault.rawValue)
+    let context = CGContext(data: nil, width: Int(width), height: Int(height), bitsPerComponent: 8, bytesPerRow: 0, space: colorSpace, bitmapInfo: bitmapInfo.rawValue)
+
+    context?.draw(image.cgImage!, in: imageRect)
+    let imageRef = context!.makeImage()
+    let newImage = UIImage(cgImage: imageRef!)
+    return newImage
+}
+
+let gray = convertToGrayScale(image: rgb)
+
+
+var mf_rgb = Matft.image.cgimage2mfarray(rgb.cgImage!)
+var mf_gray = Matft.image.cgimage2mfarray(gray.cgImage!)
+
+
+mf_rgb = mf_rgb[Matft.reverse].to_contiguous(mforder: .Row)
+mf_gray = mf_gray[Matft.reverse]
+
+Matft.image.mfarray2cgimage(mf_rgb)
+Matft.image.mfarray2cgimage(mf_gray)
+
+/*********
+    Row Contiguous
+*/
+/*
+// resize
+// ref: https://developer.apple.com/documentation/accelerate/1509208-vimagescale_argbffff
+let newdata = MfData(size: 150*150*4, mftype: .Float)
+let newstructure = MfStructure(shape: [150, 150, 4], mforder: .Row)
+newdata.withUnsafeMutableStartRawPointer{
+    dstptr in
+    mf_rgb.withUnsafeMutableStartRawPointer{
+        srcptr in
+        var src_buffer = vImage_Buffer(data: srcptr, height: vImagePixelCount(225), width: vImagePixelCount(225), rowBytes: 225*4*4)
+        var dst_buffer = vImage_Buffer(data: dstptr, height: vImagePixelCount(150), width: vImagePixelCount(150), rowBytes: 150*4*4)
+
+        vImageScale_ARGBFFFF(&src_buffer, &dst_buffer, nil, vImage_Flags(kvImageHighQualityResampling))
+    }
+
+}
+
+let resize = MfArray(mfdata: newdata, mfstructure: newstructure)
+print(resize)
+Matft.image.mfarray2cgimage(resize)
+*/
+
+/*********
+    Column Contiguous
+*/
+mf_rgb = mf_rgb.to_contiguous(mforder: .Column)
+
+// resize
+// ref: https://developer.apple.com/documentation/accelerate/1509208-vimagescale_argbffff
+let newdata = MfData(size: 150*150*4, mftype: .Float)
+let newstructure = MfStructure(shape: [150, 150, 4], mforder: .Column)
+newdata.withUnsafeMutableStartRawPointer{
+    dstptr in
+    mf_rgb.withUnsafeMutableStartRawPointer{
+        srcptr in
+        for i in 0..<4{
+            var src_buffer = vImage_Buffer(data: srcptr + i*225*225*4, height: vImagePixelCount(225), width: vImagePixelCount(225), rowBytes: 225*4)
+            var dst_buffer = vImage_Buffer(data: dstptr + i*150*150*4, height: vImagePixelCount(150), width: vImagePixelCount(150), rowBytes: 150*4)
+
+            vImageScale_PlanarF(&src_buffer, &dst_buffer, nil, vImage_Flags(kvImageHighQualityResampling))
+        }
+
+    }
+
+}
+
+var resize = MfArray(mfdata: newdata, mfstructure: newstructure)
+Matft.image.mfarray2cgimage(resize)
+
+resize.swapaxes(axis1: -1, axis2: 0)[3] = MfArray([0.5] as [Float])
+print(resize)
+let alpha = resize[Matft.all, Matft.all, 3~<4]
+resize[0~<, 0~<, 0~<3] = (resize[Matft.all, Matft.all, 0~<3]*alpha + (1 - alpha) * MfArray([1, 1, 1], mftype: resize.mftype))
+print(resize)
+
+
+
+let grayed = Matft.image.color(resize, exclude_alpha: false)
+Matft.image.mfarray2cgimage(grayed)
+
+
+/*
+let colorSpace = CGColorSpaceCreateDeviceRGB()
+let bitmapInfo = CGBitmapInfo(rawValue: CGImageAlphaInfo.premultipliedFirst.rawValue | CGImageByteOrderInfo.order32Little.rawValue | CGBitmapInfo.floatComponents.rawValue)
+print(resize.strides)
+resize.withUnsafeMutableStartRawPointer{
+    srcptr -> CGImage in
+    // NOTE: Force cast to UInt8
+    let provider = CGDataProvider(data: CFDataCreate(kCFAllocatorDefault, srcptr.assumingMemoryBound(to: UInt8.self), 150*150*4*4))
+    let cgimage =  CGImage(width: 150, height: 150, bitsPerComponent: 8*4, bitsPerPixel: 8*4*150*150, bytesPerRow: 4, space: colorSpace, bitmapInfo: bitmapInfo, provider: provider!, decode: nil, shouldInterpolate: false, intent: CGColorRenderingIntent.defaultIntent)!
+
+    return cgimage
+}
+
+*/
+
+
+*/

MatftDemo/MatftDemo/accelerate.playground/contents.xcplayground

@@ -0,0 +1,4 @@
+<?xml version="1.0" encoding="UTF-8" standalone="yes"?>
+<playground version='5.0' target-platform='ios' buildActiveScheme='true' importAppTypes='true'>
+    <timeline fileName='timeline.xctimeline'/>
+</playground>

Package.swift

@@ -21,14 +21,18 @@ let package = Package(
     targets: [
         // Targets are the basic building blocks of a package. A target can define a module or a test suite.
         // Targets can depend on other targets in this package, and on products in packages which this package depends on.
+        .target(
+            name: "pocketFFT"
+            ),
         .target(
             name: "Matft",
-            dependencies: ["Collections"]),
+            dependencies: ["Collections", "pocketFFT"]),
         .testTarget(
             name: "MatftTests",
             dependencies: ["Matft"]),
         .testTarget(
             name: "PerformanceTests",
             dependencies: ["Matft"]),
     ]
+    //cxxLanguageStandard: .gnucxx1z
 )

README.md

@@ -687,6 +687,13 @@ Below is Matft's function list. As I mentioned above, almost functions are simil
 | Matft.complex.conjugate         | numpy.conj / numpy.conjugate   |
 | Matft.complex.abs         | numpy.abs / numpy.absolute   |
 
+- FFT
+
+| Matft                            | Numpy              |
+| -------------------------------- | ----------------- |
+| Matft.fft.rfft       | numpy.fft.rfft |
+| Matft.fft.irfft         | numpy.fft.irfft   |
+
 - Interpolation
 
 Matft supports only natural cubic spline. I'll implement other boundary condition later.

Sources/Matft/Matft.swift

@@ -45,6 +45,11 @@ public class Matft{
      */
     public class image{}
 
+    /**
+       FFT
+     */
+    public class fft{}
+
     /**
        Complex
      */

Sources/Matft/core/object/mfarray.swift

@@ -212,6 +212,15 @@ extension MfArray{
     public var isComplex: Bool{
         return !self.mfdata._isReal
     }
+
+    internal var mfdata_base: MfData{
+        if self.isView{
+            return self.base!.mfdata_base
+        }
+        else{
+            return self.mfdata
+        }
+    }
 }
 
 

Sources/Matft/core/object/mfdata.swift

@@ -16,7 +16,7 @@ internal enum MfDataSource{
 
 public class MfData: MfDataProtocol{
     internal var _base: MfDataBasable? // must be referenced because refdata could be freed automatically?
-    private var _fromOtherDataSource: Bool = false
+    internal var _fromOtherDataSource: Bool = false
     internal var data_real: UnsafeMutableRawPointer
     internal var data_imag: UnsafeMutableRawPointer?
 

Sources/Matft/core/protocol/mftypeProtocol.swift

@@ -194,8 +194,8 @@ extension Bool: MfBinary, StoredFloat {
 }
 
 public protocol MfStorable: MfTypable, FloatingPoint{
-    associatedtype vDSPType: vDSP_ComplexTypable
-    associatedtype blasType: blas_ComplexTypable
+    associatedtype vDSPComplexType: vDSP_ComplexTypable
+    associatedtype blasComplexType: blas_ComplexTypable
 
     static func num(_ number: Int) -> Self
     static func from<T: MfTypable>(_ value: T) -> Self
@@ -205,8 +205,8 @@ public protocol MfStorable: MfTypable, FloatingPoint{
 }
 
 extension Float: MfStorable{
-    public typealias vDSPType = DSPSplitComplex
-    public typealias blasType = DSPComplex
+    public typealias vDSPComplexType = DSPSplitComplex
+    public typealias blasComplexType = DSPComplex
 
     public static func num(_ number: Int) -> Float {
         return Float(number)
@@ -252,8 +252,8 @@ extension Float: MfStorable{
     }
 }
 extension Double: MfStorable{
-    public typealias vDSPType = DSPDoubleSplitComplex
-    public typealias blasType = DSPDoubleComplex
+    public typealias vDSPComplexType = DSPDoubleSplitComplex
+    public typealias blasComplexType = DSPDoubleComplex
 
     public static func num(_ number: Int) -> Double {
         return Double(number)

Sources/Matft/core/util/pointer/withptr.swift

@@ -36,7 +36,7 @@ extension MfArray{
         return ret
     }
 
-    public func withUnsafeMutablevDSPPointer<T: vDSP_ComplexTypable, R>(datatype: T.Type, _ body: (UnsafeMutablePointer<T>) throws -> R) rethrows -> R{
+    public func withUnsafeMutablevDSPComplexPointer<T: vDSP_ComplexTypable, R>(datatype: T.Type, _ body: (UnsafeMutablePointer<T>) throws -> R) rethrows -> R{
 
         let ret = try self.withUnsafeMutableStartPointer(datatype: T.T.self){ ptrrT in
             return try self.withUnsafeMutableStartImagPointer(datatype: T.T.self){
@@ -48,9 +48,9 @@ extension MfArray{
 
         return ret
     }
-    internal func withUnsafeMutableblasPointer<T: blas_ComplexTypable, R>(datatype: T.Type, vDSP_func: vDSP_convert_func<T.vDSPType, T>, _ body: (UnsafeMutablePointer<T>) throws -> R) rethrows -> R{
+    internal func withUnsafeMutableBlasComplexPointer<T: blas_ComplexTypable, R>(datatype: T.Type, vDSP_func: vDSP_convert_func<T.vDSPType, T>, _ body: (UnsafeMutablePointer<T>) throws -> R) rethrows -> R{
 
-        let ret = try self.withUnsafeMutablevDSPPointer(datatype: T.vDSPType.self){ [unowned self](ptr) -> R in
+        let ret = try self.withUnsafeMutablevDSPComplexPointer(datatype: T.vDSPType.self){ [unowned self](ptr) -> R in
             var arr = Array(repeating: T(real: T.T.zero, imag: T.T.zero), count: self.storedSize)
             wrap_vDSP_convert(arr.count, ptr, 1, &arr, 1, vDSP_func)
             return try body(&arr)
@@ -119,7 +119,7 @@ extension MfData{
         return ret
     }
 
-    public func withUnsafeMutablevDSPPointer<T: vDSP_ComplexTypable, R>(datatype: T.Type, _ body: (UnsafeMutablePointer<T>) throws -> R) rethrows -> R{
+    public func withUnsafeMutablevDSPComplexPointer<T: vDSP_ComplexTypable, R>(datatype: T.Type, _ body: (UnsafeMutablePointer<T>) throws -> R) rethrows -> R{
 
         let ret = try self.withUnsafeMutableStartPointer(datatype: T.T.self){ ptrrT in
             return try self.withUnsafeMutableStartImagPointer(datatype: T.T.self){
@@ -131,9 +131,9 @@ extension MfData{
 
         return ret
     }
-    internal func withUnsafeMutableblasPointer<T: blas_ComplexTypable, R>(datatype: T.Type, vDSP_func: vDSP_convert_func<T.vDSPType, T>, _ body: (UnsafeMutablePointer<T>) throws -> R) rethrows -> R{
+    internal func withUnsafeBlasComplexPointer<T: blas_ComplexTypable, R>(datatype: T.Type, vDSP_func: vDSP_convert_func<T.vDSPType, T>, _ body: (UnsafeMutablePointer<T>) throws -> R) rethrows -> R{
 
-        let ret = try self.withUnsafeMutablevDSPPointer(datatype: T.vDSPType.self){ [unowned self](ptr) -> R in
+        let ret = try self.withUnsafeMutablevDSPComplexPointer(datatype: T.vDSPType.self){ [unowned self](ptr) -> R in
             var arr = Array(repeating: T(real: T.T.zero, imag: T.T.zero), count: self.storedSize)
             wrap_vDSP_convert(arr.count, ptr, 1, &arr, 1, vDSP_func)
             return try body(&arr)

Sources/Matft/function/method/conversion+method.swift

@@ -25,6 +25,8 @@ extension MfArray{
             - isplace: Whether to operate in-place
     */
     internal func to_complex(_ inplace: Bool = true) -> MfArray{
+        precondition(!self.mfdata._fromOtherDataSource, "Other data source couldn't be converted into Complex.")
+
         if self.isComplex{
             return self
         }
@@ -40,9 +42,11 @@ extension MfArray{
         switch mfarray.storedType{
         case .Float:
             let ptri = allocate_unsafeMRPtr(type: Float.self, count: mfarray.storedSize)
+            mfarray.mfdata_base.data_imag = ptri
             mfarray.mfdata.data_imag = ptri
         case .Double:
             let ptri = allocate_unsafeMRPtr(type: Double.self, count: mfarray.storedSize)
+            mfarray.mfdata_base.data_imag = ptri
             mfarray.mfdata.data_imag = ptri
         }
         return mfarray