将 Matlab 中的快速 MCLT 算法转换为 java

Question

我有 java 代码，它为我提供了真实输入的 FFT 输出。我需要执行 MCLT。目前我有以下格式的 FFT 输出。我见过一些用 Matlab 编码的快速 MCLT 算法 (https://www.microsoft.com/en-us/research/wp-content/uploads/2016/02/tr-2005-02.pdf)，但不能完全理解。谁能帮我写对应的java代码。

Java 代码起点：

int dtLength =  data.length/2;
double[] realPart = new double[dtLength];
double[] imagPart = new double[dtLength];

Matlab 代码：

function X = fmclt(x)
% FMCLT - Compute MCLT of a vector via double-length FFT
%
% H. Malvar, September 2001 -- (c) 1998-2001 Microsoft Corp.
%
% Syntax: X = fmclt(x)
%
% Input: x : real-valued input vector of length 2*M
%
% Output: X : complex-valued MCLT coefficients, M subbands
% in Matlab, by default j = sqrt(-1)
% determine # of subbands, M
L = length(x);
M = L/2;
% normalized FFT of input
U = sqrt(1/(2*M)) * fft(x);
% compute modulation function
k = [0:M]';
c = W(8,2*k+1) .* W(4*M,k);
% modulate U into V
V = c .* U(1:M+1);
% compute MCLT coefficients
X = j * V(1:M) + V(2:M+1);
return;
% Local function: complex exponential
function w = W(M,r)
w = exp(-j*2*pi*r/M);
return; 

Answer 1

尽管这个问题对于 SO 来说有点边界，但这篇论文很有趣，所以我决定花一些时间阅读它并尝试将 Matlab 代码转换为 Java。结果如下：

import org.apache.commons.math3.complex.Complex;

public class MCLT
{
    public static void main(String args[])
    {
        Complex[] x = new Complex[16];

        for (int i = 1; i <= 16; ++i)
            x[(i - 1)] = new Complex((double)i, 0.0d);

        Complex[] result = fmclt(x);

        for (int i = 0; i < result.length; ++i)
            System.out.println(result[i]);
    }

    public static Complex[] fmclt(Complex[] x)
    {
        int L = x.length;
        int M = L / 2;

        double z = Math.sqrt(1.0d / (2.0d * M));

        Complex[] F = fft(x);
        Complex[] U = new Complex[F.length];

        for (int i = 0; i < F.length; ++i)
            U[i] = F[i].multiply(z);

        double[] k = new double[(M + 1)];

        for (int i = 0; i <= M; ++i)
            k[i] = (double)i;

        Complex[] c = new Complex[(M + 1)];

        for (int i = 0; i <= M; ++i)
            c[i] = W(8.0d, ((2.0d * k[i]) + 1.0d)).multiply(W((4.0d * M), k[i]));

        Complex[][] V = new Complex[(M + 1)][];

        for (int i = 0; i <= M; ++i)
        {
            V[i] = new Complex[(M + 1)];

            for (int j = 0; j <= M; ++j)
                V[i][j] = c[i].multiply(U[j]);
        }

        Complex[] V1 = new Complex[M];

        for (int i = 0; i < M; ++i)
            V1[i] = V[i][0];

        Complex[] V2 = new Complex[M];

        for (int i = 1; i <= M; ++i)
            V2[(i - 1)] = V[i][0];

        Complex b = new Complex(0.0d, 1.0d);
        Complex[] result = new Complex[M];

        for (int i = 0; i < M; ++i)
            result[i] = b.multiply(V1[i]).add(V2[i]); 

        return result;
    }

    public static Complex[] fft(Complex[] x)
    {
        int n = x.length;

        if (n == 1)
            return new Complex[] { x[0] };

        if ((n % 2) != 0)
            throw new IllegalArgumentException("Invalid length.");

        int nh = n / 2;

        Complex[] even = new Complex[nh];

        for (int i = 0; i < nh; ++i)
            even[i] = x[(2 * i)];

        Complex[] q = fft(even);

        Complex[] odd  = even;

        for (int i = 0; i < nh; ++i)
            odd[i] = x[((2 * i) + 1)];

        Complex[] r = fft(odd);

        Complex[] y = new Complex[n];

        for (int i = 0; i < nh; ++i)
        {
            double kth = -2.0d * i * (Math.PI / n);
            Complex wk = new Complex(Math.cos(kth), Math.sin(kth));

            y[i] = q[i].add(wk.multiply(r[i]));
            y[(i + nh)] = q[i].subtract(wk.multiply(r[i]));
        }

        return y;
    }

    public static Complex W(double M, double r)
    {
        Complex j = (new Complex(0.0d, 1.0d)).multiply(-1.0d);
        double z = 2.0d * Math.PI * (r / M);

        return j.multiply(z).exp();
    }
}

在我看来，对实部和虚部使用单独的双数组并不是一个好的设计选择，所以我决定将我的代码基于 Apache Commons 库的 Complex 类。

为了计算快速傅里叶变换，我决定使用一些现成的代码。我的fft函数是基于this implementation，看起来很可靠，利用了前面提到的Complex类。

使用相同的值向量，Matlab 和 Java 代码返回相同的输出。您可以通过将代码复制粘贴到this website 来在线测试代码，但您还需要安装Apache Commons 库才能成功运行它。点击底部的Add External Library (from Maven Repo)按钮，然后在输入表单中插入以下参数：

<!-- https://mvnrepository.com/artifact/org.apache.commons/commons-math3 -->
<dependency>
    <groupId>org.apache.commons</groupId>
    <artifactId>commons-math3</artifactId>
    <version>3.6.1</version>
</dependency>