Python转C语言

Question

import numpy as np

N = 512  # 数据长度
fs = 200  # 采样Frequency
t = np.arange(N) / fs  # 时间刻度

# 构成信号序列
x = 10 * np.sin(2 * np.pi * 32 * t) + 10 * np.sin(2 * np.pi * 50 * t) + \
    20 * np.sin(2 * np.pi * 54 * t) + 20 * np.sin(2 * np.pi * 56 * t) + \
    30 * np.sin(2 * np.pi * 59 * t) + 20 * np.sin(2 * np.pi * 83 * t)

nfft = 128  # FFT长度
fe = 55  # 中心Frequency
D = 20  # 细化倍数

# 细化分析
def exzfft_ma(x, fe, fs, nfft, D):
    nt = len(x)  # 计算读入数据长度
    fi = fe - fs / (2 * D)  # 计算细化截止频率下限
    fa = fi + fs / D  # 计算细化截止频率上限
    na = round(0.5 * nt / D + 1)  # 确定低通滤波器截止频率对应的谱线条数

    # 频移
    n = np.arange(nt)  # 序列索引号
    b = n * np.pi * (fi + fa) / fs  # 设置单位旋转因子
    y = x * np.exp(-1j * b)  # 进行频移

    b = np.fft.fft(y, nt)  # FFT

    # 低通滤波和下采样
    a = np.zeros(nt, dtype=complex)
    a[:na] = b[:na]  # 取正频率部分的低频成分
    a[nt - na + 1:] = b[nt - na + 1:]  # 取负频率部分的低频成分

    b = np.fft.ifft(a, nt)  # IFFT
    c = b[:nt:D]  # 下采样

    # 求细化频谱
    y = np.fft.fft(c, nfft) * 2 / nfft  # 再一次FFT
    y = np.fft.fftshift(y)  # 重新排列
    freq = fi + np.arange(nfft) * fs / (D * nfft)  # 频率设置
    return y * D , freq

y, freq = exzfft_ma(x, fe, fs, nfft, D)
print(y, freq)

将上面Python代码转换为C语言.

Answer 1

这个试试看

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#if 1


#define N 512   //数据长度
#define fs 200  //采样频率
#define D 20    //细化倍数
#define fe 55   //中心Frequency
#define nfft 128    //FFT长度

#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

typedef struct {
    float real;
    float imag;
} Complex;


Complex complex_add(Complex a, Complex b) {
    Complex c = { a.real + b.real, a.imag + b.imag };
    return c;
}

Complex complex_subtract(Complex a, Complex b) {
    Complex c = { a.real - b.real, a.imag - b.imag };
    return c;
}

Complex complex_multiply(Complex a, Complex b) {
    Complex c = { a.real * b.real - a.imag * b.imag, a.real * b.imag + a.imag * b.real };
    return c;
}

void fft(Complex* x, int n) {
    if (n <= 1) {
        return;
    }

    Complex* even = (Complex*)malloc(n / 2 * sizeof(Complex));
    Complex* odd = (Complex*)malloc(n / 2 * sizeof(Complex));

    for (int i = 0; i < n / 2; i++) {
        even[i] = x[2 * i];
        odd[i] = x[2 * i + 1];
    }

    fft(even, n / 2);
    fft(odd, n / 2);

    for (int k = 0; k < n / 2; k++) {
        Complex t = { cos(2 * M_PI * k / n), -sin(2 * M_PI * k / n) };
        Complex u = complex_multiply(t, odd[k]);

        x[k] = complex_add(even[k], u);
        x[k + n / 2] = complex_subtract(even[k], u);
    }

    free(even);
    free(odd);
}

void ifft(Complex* x, int n) {
    for (int i = 0; i < n; i++) {
        x[i].imag = -x[i].imag;
    }

    fft(x, n);

    for (int i = 0; i < n; i++) {
        x[i].real /= n;
        x[i].imag = -x[i].imag / n;
    }
}
void exzfft_ma(float* x, float* y, float* freq, int nt, float fi, float fa, int na) {
    Complex* b = (Complex*)malloc(nt * sizeof(Complex));
    Complex* a = (Complex*)malloc(nt * sizeof(Complex));

    for (int i = 0; i < nt; i++) {
        float n = i;
        float b_val = n * M_PI * (fi + fa) / fs;
        Complex exp_val = { cos(-b_val), sin(-b_val) };
        Complex x_val = { x[i], 0 };
        a[i] = complex_multiply(x_val, exp_val);
    }

    fft(a, nt);

    for (int i = 0; i < nt; i++) {
        if (i < na || i >= nt - na + 1) {
            b[i] = a[i];
        }
        else {
            b[i].real = 0;
            b[i].imag = 0;
        }
    }

    ifft(b, nt);

    for (int i = 0; i < nt; i++) {
        y[i] = b[i].real;
    }

    int mid = nfft / 2;
    int shift = nt / 2;
    for (int i = 0; i < nfft; i++) {
        int idx = (i + shift) % nt;
        y[i] = y[idx];
    }

    float step = fs / (D * nfft);
    for (int i = 0; i < nfft; i++) {
        freq[i] = fi + i * step;
    }

    free(b);
    free(a);
}




int main() {
    float t[N];
    float x[N];
    Complex* y = (Complex*)malloc(nfft * sizeof(Complex));
    float freq[N];

    for (int i = 0; i < N; i++) {
        t[i] = i / (float)fs;
    }

    for (int i = 0; i < N; i++) {
        x[i] = 10 * sin(2 * M_PI * 32 * t[i]) + 10 * sin(2 * M_PI * 50 * t[i]) +
            20 * sin(2 * M_PI * 54 * t[i]) + 20 * sin(2 * M_PI * 56 * t[i]) +
            30 * sin(2 * M_PI * 59 * t[i]) + 20 * sin(2 * M_PI * 83 * t[i]);
    }

    int nt = N;
    float fi = fe - fs / (2 * D);
    float fa = fi + fs / D;
    int na = round(0.5 * nt / D + 1);
    exzfft_ma(x, y, freq, nt, fi, fa, na);

    // Print the results
    for (int i = 0; i < N; i++) {
        printf("y[%d] = %f\n", i, y[i]);
    }

    for (int i = 0; i < N; i++) {
        printf("freq[%d] = %f\n", i, freq[i]);
    }

    return 0;
}
#endif // 0