3
import Cocoa
import Accelerate

let filePath = Bundle.main.path(forResource: "sinusoid", ofType: "txt")
let contentData = FileManager.default.contents(atPath: filePath!)
var content = NSString(data: contentData!, encoding: String.Encoding.utf8.rawValue) as? String

var idx = content?.characters.index(of: "\n")
idx = content?.index(after: idx!)

repeat {
    //let fromIndex = index(from: )
    content = content?.substring(from: idx!)
    idx = content?.characters.index(of: "\n")
    idx = content?.index(after: idx!)
} while content!.characters.contains("%")

let regex = try? NSRegularExpression(pattern: "[ ]+", options:[])

let delimiter = ","
var modifiedString = regex?.stringByReplacingMatches(in: content!, options: [], range: NSRange(location: 0, length: (content! as NSString).length), withTemplate: delimiter)

let lines = modifiedString?.components(separatedBy: "\n")

var s = [Double]()

for var line in lines! {
    if !line.isEmpty {
        let data = line.components(separatedBy: ",")
        s.append(Double(data[1])!)
    }
}

let length = vDSP_Length(pow(2, floor(log2(Float(s.count)))))
let L = Int(length)

// zrop or zop? 
// zrop covers real to complex, and zop covers complex
// length must be a power of 2 or specific multiples of powers of 2 if size is at least 4
let setup = vDSP_DFT_zrop_CreateSetupD(nil, length, vDSP_DFT_Direction.FORWARD)

var inputReal = UnsafeMutablePointer<Double>.allocate(capacity: L)
var inputImaginary = UnsafeMutablePointer<Double>.allocate(capacity: L)
var outputReal = UnsafeMutablePointer<Double>.allocate(capacity: L)
var outputImaginary = UnsafeMutablePointer<Double>.allocate(capacity: L)

for i in 0..<L {
    inputReal[i] = s[i]
    inputImaginary[i] = 0.0
}

vDSP_DFT_ExecuteD(setup!, inputReal, inputImaginary, outputReal, outputImaginary)

for i in 0..<L {
    print("\(outputReal[i]) + \(outputImaginary[i])i")
}

The input file "sinusoid.txt" is in the following link https://dpaste.de/M1VD

The input file data consists of two sine waves at frequencies of 50 and 120. The Matlab code produces the correct output given in the following link:

https://dpaste.de/2mdK

When the result from Matlab is scaled and the magnitude is taken, it correctly shows that the amplitude at a frequency of 50 is 0.7 and the amplitude at a frequency of 120 is 1.

clear all; close all; clc;
data = load('sinusoid.txt');
S = data(:,2);
Fs = 1000;
Y = fft(S);
L = length(S);
P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = Fs*(0:(L/2))/L;
plot(f,P1)
title('Single-Sided Amplitude Spectrum of X(t)')
xlabel('f (Hz)')
ylabel('|P1(f)|')

The Swift code output is entirely different and unrecognizable when compared to the Matlab output, regardless of whatever scaling factors are applied and whether or not a real-to-complex or complex-to-complex transformation is applied:

https://dpaste.de/MUwB

Any ideas why this is?

SwiftMatt
  • 949
  • 1
  • 9
  • 19
  • 2
    Your input data has 1000 numbers, but only 512 of them are passed into the DSP function and the remaining 488 values are simply ignored. How can that work?? – Martin R Jan 07 '17 at 09:45
  • What do you mean how can that work? The FFT function in Swift requires a power of 2, so it needs to be 512 – SwiftMatt Jan 07 '17 at 09:55
  • Trying it with 1024 – SwiftMatt Jan 07 '17 at 09:56
  • Where is the MATLAB code for comparison ? – Paul R Jan 07 '17 at 09:57
  • Ok I added the MATLAB code as well – SwiftMatt Jan 07 '17 at 09:59
  • 1
    I am not familiar with MATLAB, but it seems to be that the 1000 element data set is *not* truncated to 512 in your MATLAB code. So you cannot expect the same result. – Martin R Jan 07 '17 at 10:01
  • Ok well if rather than truncating, I pad with 0s to 1024 data elements, there is still no same result Martin R – SwiftMatt Jan 07 '17 at 10:04
  • 1
    As I understand it, the FFT functions from the Accelerate framework *can only* handle sample sizes which are a power of two times 1, 3, 5, or 15, while MATLAB handles arbitrary sizes. Related: http://stackoverflow.com/q/12134208/1187415, http://stackoverflow.com/q/10708667/1187415 and http://math.stackexchange.com/q/77118/42969 from MSE. As you can see from the last reference, transforming the "arbitrary size" case to a "power of two" case is not trivial, you cannot simply pad the samples with zeros. – Martin R Jan 07 '17 at 16:12
  • Hmm, okay thanks, I will try this – SwiftMatt Jan 07 '17 at 18:35

2 Answers2

1

The lengths of your 2 FFT are different, and so, of course the results won't match. You are also passing different amounts of data to your 2 FFTs.

Print out the FFT lengths and the input data vector to debug your code. Make sure the inputs match before comparing results.

Also, Apple's Accelerate/vDSP FFTs can use lengths other than just powers of 2 (lengths with factors of 3 or 5 are also allowed).

Also, note that Matlab indexes arrays starting at 1, not 0, as is more typical in C and Swift functions.

hotpaw2
  • 70,107
  • 14
  • 90
  • 153
0

Indeed, the problem with the FFT result mismatch was due to the input size mismatch. Restricting input to be specific multiples of powers of 2 greatly constrains the FFT usage in the Accelerate framework. One suggestion was to pad the input with 0s until it is of an appropriate length. Whether you pad with 0s or truncate the input such that the size is a specific multiple of a power of 2, the results from the Accelerate framework will be different than the results from a program such as MATLAB. The solution to this is to perform a chirp-z transform as mentioned in a link specified by Martin R. The chirp-z transform itself yields identical results to the FFT and may be performed on inputs of arbitrary size.

SwiftMatt
  • 949
  • 1
  • 9
  • 19