FFT变换的C++代码-MATLAB代码

1FFT变换C++代码

public static double[] FFT(int N, double[] Data, int NFFT)

 存储返回结果前面存储实部 后面存储虚部

        int n = NFFT / 2;                   //FFT变换所需要的点数目的一半

        int m=0;                              //蝶形计算的级数

        double[] x = new double[n];         //存储采样的偶数数据

        double[] y = new double[n];         //存储采样的奇数数据

        int n1;                             //倒序运算

        int k;                              //倒序运算所需变量

        double tr, ti;                       //倒序运算所需中间变量

 数据长度不够，末尾填零

 进行奇偶分解，分解为另个序列

 计算蝶形运算的级数

 倒序运算

 蝶型运算

 计算旋转因子，何不直接套用公式 =cos(2πP/N)-sin(2πP/N)，P=J* ,会更简单

 提取运算结果 将x[i],y[i]的FFT变换的结果的实部与虚部分别存储在fr gr fi gi中

2FFT变换C++代码在MATLAB下的实现

function Y = FFT_WYCFun( N,X ,NFFT)

n = NFFT/2; %FFT变换点数的一半

m = log2(n);%蝶型算法的迭代次数

EvenData = zeros(n,1);%用于存储原始数据的计数

OddData = zeros(n,1); %用于存储原始数据的偶数

 % 汉明窗

u=0;

twopi = 8.0*atan(1.0);

f1 = n-1;

win = zeros(n,1);

for i=1:n

    win(i) = 0.54-0.46*cos(twopi*(i-1)/f1);

    u = u + win(i)*win(i);

end

 %数据不够，补0

for i=N+1:NFFT

    X(i) = 0;

End

%分奇数与偶数

for a =1:n

    EvenData(a) = X(2*a-1);%*win(a);

    OddData(a) = X(2*a);%*win(a);

end

%倒序运算

n1 = n-1;

j = 1;

for i=1:n1

    if(i        tr = EvenData(j);

        ti = OddData(j);

        EvenData(j) = EvenData(i);

        OddData(j) = OddData(i);

        EvenData(i) = tr;

        OddData(i) = ti;

    end

    k = n/2;

    while(k <(j+1))

        j = j-k;

        k = k/2;

    end

     j = j+k;

end

%蝶型算法

n1 = 1;

for  l=1:m

    n1 = 2*n1;

    n2 = n1 /2;

    e = 3.14159265359/n2;

    c = 1;

    s = 0;

    c1 = cos(e);

    s1 = -1*sin(e);

    for j = 1:n2

        for i = j:n1:n

            k = i +n2;

            tr = c*OddData(k) - s*EvenData(k);

            ti = c*EvenData(k) + s*OddData(k);

            OddData(k) = OddData(i) - tr;

            EvenData(k) = EvenData(i) - ti;

            OddData(i) = OddData(i) + tr;

            EvenData(i) = EvenData(i) + ti;

        end

        t = c;

        c = c*c1 - s*s1;

        s = t*s1 + s*c1;

    end

end

er = 3.14159265359/n1;

%结果提取

Y=zeros(NFFT,1);

Y(1) = OddData(1) + EvenData(1);

Y(n1+1) = OddData(1) - EvenData(1);

for i=1:n1-1

    fr = (OddData(i) + OddData(n1-i))/2;

    fi = (EvenData(i) - EvenData(n1-i))/2;

    gr = (EvenData(n1-i) + EvenData(i))/2;

    gi = (OddData(n1-i) - OddData(i))/2;

    a = i*er;

    p = cos(a);

    q = sin(a);

    Y(i) = fr + p*gr +q*gi;

    Y(NFFT - i) = fi + p*gi -q*gr;

end

% norm = 2.0*u;

% for i=1:n

%     Y(i) = Y(i)/norm;

% end

3数据长度及FFT点数对FFT变换结果的影响

根据之前讨论过的对原始接地电流信号进行重采样，不同的重采样率对FFT变换（MATLAB自带算法）的结果将会产生影响；最终经过对中重采样率的对比，确定一下两者重采样方式：

1)每20个点进行一次求平均，重采样后的采样率为20K，数据长度为2048；

2)每40个点进行一次球平均，重采样后的采样率为10K，数据长度为1024。

下面就对这两种重采样了进行探讨：

3.1原始文件及MATLAB下FFT变换频谱

可以看出当前信号中包含50Hz（对应幅值18mA）、150Hz（对应幅值7mA）以及幅值很小的350Hz、450Hz

3.2重采样率为20K（采样点数为2048）

3.2.1FFT点数为1024

当前FFT变换后基本看不出特征频率。

3.2.2FFT点数为1024*2

与上一中情况类似。

3.2.3FFT点数为1024*4

当前FFT变换结果中可以看到两个明显的峰值点，分别对应40Hz、140Hz左右，幅值分别为7mA、4mA；且峰值点附近的波动很大，及变换后的方差很大（特征频率附近的旁瓣很大）。

这两个峰值点分别对应的是原始信号中50Hz、150Hz，但其幅值缩小了两倍左右。因此当前FFT变换：主频（次频）提取不精确，幅值提取不准确，且旁瓣较大。

3.2.4FFT点数为1024*8

当前FFT变换结果可明显看到两峰值点，且其对应的频率更接近于本质频率（分别为50Hz、150Hz），其幅值分别8mA、3mA；且峰值点附近波动较小，旁瓣小。

因此当前FFT变化：主频（次频）提取精确，但幅值不准，旁瓣较小。

3.2.5FFT点数为1024*16

当前FFT变换结果可明显看到两峰值点，且其对应的频率接近于本质频率（分别为50Hz、150Hz），其幅值分别11mA、6mA；且峰值点附近波动较大，旁瓣较大。

因此当前FFT变化：主频（次频）提取精确，但幅值不够准确，旁瓣较大。

3.2.6FFT点数为1024*32

当前FFT变换结果可明显看到两峰值点，且其对应的频率接近于本质频率（分别为50Hz、150Hz），其幅值分别14mA、6mA；且峰值点附近波动大，旁瓣大。

因此当前FFT变化：主频（次频）提取精确，但幅值基本准确，但旁瓣较大。

3.2.7FFT点数为1024*

与上一中情况下的结果基本一致。

3.3重采样率为10K（采样点数为1024）

3.3.1FFT点数为1024

当前FFT变换后基本看不出特征频率。

3.3.2FFT点数为1024*2

当前FFT变换结果中可以看到两个明显的峰值点，分别对应40Hz、140Hz左右，幅值分别为7mA、4mA；且峰值点附近的波动很大，及变换后的方差很大（特征频率附近的旁瓣很大）。

这两个峰值点分别对应的是原始信号中50Hz、150Hz，但其幅值缩小了两倍左右。因此当前FFT变换：主频（次频）提取不精确，幅值提取不准确，且旁瓣较大。

3.3.3FFT点数为1024*4

当前FFT变换结果可明显看到两峰值点，且其对应的频率更接近于本质频率（分别为50Hz、150Hz），其幅值分别8mA、3mA；且峰值点附近波动较小，旁瓣小。

因此当前FFT变化：主频（次频）提取精确，但幅值不准，旁瓣较小。

3.3.4FFT点数为1024*8

当前FFT变换结果可明显看到两峰值点，且其对应的频率接近于本质频率（分别为50Hz、150Hz），其幅值分别11mA、6mA；且峰值点附近波动较大，旁瓣较大。

因此当前FFT变化：主频（次频）提取精确，但幅值不够准确，旁瓣较大。

3.3.5FFT点数为1024*16

当前FFT变换结果可明显看到两峰值点，且其对应的频率接近于本质频率（分别为50Hz、150Hz），其幅值分别14mA、6mA；且峰值点附近波动大，旁瓣大。

因此当前FFT变化：主频（次频）提取精确，但幅值基本准确，但旁瓣较大。

3.3.6FFT点数为1024*32

与上一中情况下的结果基本一致。

3.3.7FFT点数为1024*

与上一中情况下的结果基本一致。

3.4测试结果统计对比

当FFT变换点数时采样点数的16倍时，可以精确的提取出特征频率，且特征频率对应的幅值接近实际幅值。

4结论

当前算法能够较好的提取出特征频率，但提取特征频率对应幅值精度还不够。且最佳的FFT变换方式是，FFT变换点数为采样点数的16倍或以上，为了计算效率，建议选择16倍。

因此，结合现场实际运算及效果，选择一下参数：

重采样点数间隔：40点；

重采样后的频率：10KHz；

重采样后的点数：1024；

FFT变换的点数：1024*16.