FFT算法实现

最新推荐文章于 2024-08-21 11:22:17 发布

Awanyyy

最新推荐文章于 2024-08-21 11:22:17 发布

阅读量106

点赞数 3

文章标签：算法

本文链接：https://blog.csdn.net/Awanyyy/article/details/141167765

版权

使用Rader算法进行二进制倒序

function x = rader(x)
    lengthx = length(x);
    j = lengthx/2;
    
    for i=1 : lengthx/2-1
                      
        if i<j
            t = x(j+1);
            x(j+1) = x(i+1);
            x(i+1) = t;
        end
       
            k = lengthx / 2;
            while(j >= k)
                j = j - k;
                k = k/ 2;
            end
            if j < k
                j = j + k;
            end
    end

FFT算法实现

function H = newfft(x)
    y = rader(x);
    n = length(x);
    m =  log2(n);
    wn = exp(-j * (2 * pi / n));
    H = zeros(1, n);
    for nnum = 1 : m
        delta = 2 ^ (nnum -1);
        L = 2 ^ nnum;
        for i = 1 : n * (1/2) ^ nnum
            for k2 = 1 : L
                k1 = mod(k2 - 1, delta) + 1;
                alpha = 2 * (k1 -1);
                if k2 <= (L / 2)
                    H(L * (i - 1) + k2) = y(L/2 * (2 * i - 1 - 1) + k1) + y(L/2 * (2 * i - 1) + k1) * (wn ^ alpha);
                else
                    H(L * (i - 1) + k2) = y(L/2 * (2 * i - 1 - 1) + k1) - y(L/2 * (2 * i - 1) + k1) * (wn ^ alpha);
                end
            end
        end
        y = H;
    end
end
% H(L * (i - 1) + k2)把H序列拆分成几个长为L的子序列，然后选择第i个子序列的第k2个元素

结果分析，与官方FFT对比

clear;clc;
x = [1, 2, 3, 1, 4, 1, 10, 1];
H1 = newfft(x);
H2 = fft(x);
n = 0 : length(x) - 1;
error = sum((abs(H1) - abs(H2)).^2)/length(x);
mean = sum(abs(H2))/length(x);
com_e = error/mean;
subplot(2, 1, 1);stem(n, abs(H1));
subplot(2, 1, 2);stem(n, abs(H2));

频谱采样点均方误差为1.1022

相对误差为0.1067

FFT频谱结果（只考虑幅度谱）