DFT result in Swift is different than that of MATL

2019-02-24 11:17发布

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?

2条回答
Summer. ? 凉城
2楼-- · 2019-02-24 11:47

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.

查看更多
可以哭但决不认输i
3楼-- · 2019-02-24 11:51

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.

查看更多
登录 后发表回答