C++ 傅里叶频谱的计算以及应用

一维傅里叶频谱的计算

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

#define pi 3.1415926
#define rows 3
#define colums 5

typedef struct {
    float re;// really 
    float im;// imaginary 
}complex, * pcomplex;

complex complexadd(complex a, complex b) //复数加
{
    complex rt;
    rt.re = a.re + b.re;
    rt.im = a.im + b.im;
    return rt;
}

complex complexMult(complex a, complex b) //复数乘
{
    complex rt;
    rt.re = a.re * b.re - a.im * b.im;
    rt.im = a.im * b.re + a.re * b.im;
    return rt;
}

/// <summary>
/// 离散傅里叶变换
/// </summary>
/// <param name="f">时域数据</param>
/// <param name="F">频域数据</param>
/// <param name="N">数据长度</param>
void dft(complex f[], complex F[], int N)
{
    complex temp;
    int k, n;
    for (int k = 0; k < N; k++) {
        F[k].re = 0;
        F[k].im = 0;
        for (int n = 0; n < N; n++) {
            temp.re = (float)cos(2 * pi * k * n / N);
            temp.im = -(float)sin(2 * pi * k * n / N);
            F[k] = complexadd(F[k], complexMult(f[n], temp));
        }
    }
}

/// <summary>
/// 离散傅里叶逆变换
/// </summary>
/// <param name="F">频域数据</param>
/// <param name="f">时域数据</param>
/// <param name="N">数据长度</param>
void idft(complex F[], complex f[], int N) {
    complex temp;
    int k, n;
    for (int k = 0; k < N; k++) {
        f[k].re = 0;
        f[k].im = 0;
        for (int n = 0; n < N; n++) {
            temp.re = (float)cos(2 * pi * k * n / N);
            temp.im = (float)sin(2 * pi * k * n / N);
            f[k] = complexadd(f[k], complexMult(F[n], temp));
        }
        f[k].re /= N;
        f[k].im /= N;
    }
}

/// <summary>
/// 将零频分量移到频谱中心
/// </summary>
/// <param name="F">频域数据</param>
/// <param name="Shift">零频分量移到频谱中心的数据</param>
/// <param name="N">数据长度</param>
void fftshift(complex F[], complex Shift[], int N) {
    for (int i = 0; i < N; i++) {
        int shiftIndex = (i + (N - 1) / 2 + 1) % N;
        Shift[i] = F[shiftIndex];
    }
}

/// <summary>
/// 计算功率谱函数
/// </summary>
/// <param name="F">零频分量移到频谱中心的数据</param>
/// <param name="P">功率谱数据</param>
/// <param name="N">数据长度</param>
void powerSpectrum(complex F[], float P[], int N) {
    for (int i = 0; i < N; i++) {
        P[i] = F[i].re * F[i].re + F[i].im * F[i].im;
    }
}

int main() {
    complex samples[colums], X[colums], x[colums]; //samples[]示例

    for (int i = 0; i < colums; i++) {
        samples[i].re = i + 1;
        samples[i].im = 0;
    }
    dft(samples, X, colums);
    printf("一维傅里叶变换 DFI:\n");
    for (int i = 0; i < colums; i++) {
        printf("(%f,%f)\n", X[i].re, X[i].im);
    }
    printf("\n");

    complex ShiftX[colums];
    fftshift(X, ShiftX, colums);
    printf("零频分量移到频谱中心 FFTShift:\n");
    for (int i = 0; i < colums; i++) {
        printf("(%f,%f)\n", ShiftX[i].re, ShiftX[i].im);
    }
    printf("\n");

    float PowerX[colums];
    powerSpectrum(ShiftX, PowerX, colums);
    printf("功率谱 powerSpectrum:\n");
    for (int i = 0; i < colums; i++) {
        printf("(%f)\n", PowerX[i]);
    }
    printf("\n");
 
}

运行结果