Swift 中的 DFT 结果与 MATLAB 的结果不同


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")
}

输入文件"正弦.txt"在以下链接中https://dpaste.de/M1VD

输入文件数据由频率为 50 和 120 的两个正弦波组成。Matlab 代码生成以下链接中给出的正确输出:

https://dpaste.de/2mdK

当对 Matlab 的结果进行缩放并获取幅度时,它正确地显示频率为 50 时的振幅为 0.7,频率为 120 时的振幅为 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)|')

与 Matlab 输出相比,Swift 代码输出完全不同且无法识别,无论应用了什么比例因子,也无论是否应用了实数到复杂或复杂到复杂的转换:

https://dpaste.de/MUwB

知道为什么会这样吗?

您的 2 FFT 的长度不同,因此,结果当然不会匹配。 您还将不同数量的数据传递给您的 2 个 FFT。

打印出FFT长度和输入数据向量以调试代码。 在比较结果之前,请确保输入匹配。

此外,Apple 的 Accelerate/vDSP FFT 可以使用 2 次幂以外的长度(也允许系数为 3 或 5 的长度)。

另外,请注意,Matlab 索引从 1 开始的数组,而不是 0,这在 C 和 Swift 函数中更为典型。

事实上,FFT结果不匹配的问题是由于输入大小不匹配造成的。将输入限制为 2 的特定倍数会极大地限制 FFT 在加速框架中的使用。一个建议是用 0 填充输入,直到它达到适当的长度。无论是用 0 填充还是截断输入,使大小是 2 的幂的特定倍数,加速框架的结果都将与 MATLAB 等程序的结果不同。对此的解决方案是执行 Chirp-z 转换,如 Martin R 指定的链接中所述。线性调频-z变换本身产生与FFT相同的结果,并且可以在任意大小的输入上执行。

最新更新