sosyangliu 2021-03-16 13:06 采纳率: 100%
浏览 215
已采纳

快速傅里叶逆变换未还原数据

using System;

namespace ConsoleApp34
{
    class Program
    {
        static void Main(string[] args)
        {
            const int N = 8;
            //初始化时域数据
            Complex[] TD2FD = new Complex[N];
            for (int i = 0; i < N; i++)
            {
                Complex cpx = new Complex();
                cpx.re = i;
                cpx.im = 0;
                TD2FD[i] = cpx;
            }

            Console.WriteLine("------原始信号-----");
            Print(TD2FD);

            Console.WriteLine("----- FFT -----");
            FFT(TD2FD);
            Print(TD2FD);

            Console.WriteLine("----- IFFT -----");
            IFFT(TD2FD);
            Print(TD2FD);

            Console.Read();
        }

        //快速傅里叶变换
        public static void FFT(Complex[] TD2FD)
        {
            FFT_Core(TD2FD, WT_LUT(TD2FD.Length, 1));
        }

        //快速傅里叶逆变换
        public static void IFFT(Complex[] FD2TD)
        {
            //做FFT变换
            Complex[] WT = WT_LUT(FD2TD.Length, -1);
            FFT_Core(FD2TD, WT);
            //实部除以N
            for (int i = 0; i < FD2TD.Length; i++)
                FD2TD[i].re /= FD2TD.Length;
        }

        // 返回旋转因子查询表(twiddle factor lookup table)
        private static Complex[] WT_LUT(int N, int flag = 1)
        {
            Complex[] WT = new Complex[N];
            for (int i = 0; i < N; i++)
            {
                Complex cpx_wt = new Complex();
                float angle = (float)(-i * Math.PI * 2 / N);
                cpx_wt.re = (float)Math.Cos(angle);
                //IFFT flag=-1, FFT flag=1
                cpx_wt.im = flag * (float)Math.Sin(angle);
                WT[i] = cpx_wt;
            }
            return WT;
        }

        private static void FFT_Core(Complex[] TD2FD, Complex[] WT)
        {
            int power = (int)Math.Log(TD2FD.Length, 2);
            int butterfly;
            int p, s;
            Complex x1, x2, wt;

            BitReverse(TD2FD);

            //蝶形运算
            for (int k = 0; k < power; k++) //级数
            {
                //Console.WriteLine("第{0}级, 共{1}组", k, 1 << k);
                for (int j = 0; j < 1 << k; j++) //组数
                {
                    butterfly = 1 << (power - k);//每组有几个元素
                    //计算参与蝶形运算的两个元素的索引
                    p = j * butterfly;
                    s = p + butterfly / 2;
                    //Console.WriteLine("butterfly={0}, p={1}, s={2}", butterfly, p, s);
                    for (int i = 0; i < butterfly / 2; i++) //蝶形运算次数
                    {
                        //Console.Write("  ({0}x{1} wtIdx={2})", i + p, i + s, i * (1 << k));
                        x1 = TD2FD[i + p];
                        x2 = TD2FD[i + s];
                        wt = WT[i * (1 << k)];
                        TD2FD[i + p] = x1 + x2 * wt;
                        TD2FD[i + s] = x1 - x2 * wt;
                    }
                    //Console.WriteLine();
                }
            }
        }

        private static int BitReverse(int x)
        {
            //倒位排序
            //0   1   2   3   4   5   6   7   十进制
            //000 001 010 011 100 101 110 111 十进制对应的二进制
            //000 100 010 110 001 101 011 111 码位反转
            //0   4   2   6   1   5   3   7   码位反转后对应的十进制
            int[] table = new int[8] { 0, 4, 2, 6, 1, 5, 3, 7 };
            return table[x];
        }

        // 倒位排序——雷德算法
        private static void BitReverse(Complex[] array)
        {
            int i, j, k;
            int N = array.Length;
            Complex temp;
            j = 0;

            for (i = 0; i < N - 1; i++)
            {
                if (i < j)
                {
                    temp = array[i];
                    array[i] = array[j];
                    array[j] = temp;
                }

                // 求j的下一个倒序位
                // N/2的二进制最高位为1,其他位都为0
                // 判断最高位是否为1,可与N/2进行比较
                // 判断次高位是否为1,可与N/4进行比较
                k = N >> 1;

                //如果k<=j,表示j的最高位为1
                while (k <= j)
                {
                    //当k<=j时,需要将最高位变为0
                    j = j - k;
                    //判断次高位是否为1,依次类推,逐个比较,直到j某个位为0
                    k >>= 1;
                }

                j = j + k;//将0变为1
            }
        }

        // 打印
        private static void Print(Complex[] TD2FD)
        {
            for (int i = 0; i < TD2FD.Length; i++)
            {
                Console.WriteLine("(re={0}, im={1})", TD2FD[i].re, TD2FD[i].im);
            }
            Console.WriteLine();
        }
    }

    //定义复数
    public class Complex
    {
        public float re;//实数部
        public float im;//虚数部

        public static Complex operator +(Complex lhs, Complex rhs)
        {
            Complex result = new Complex();
            result.re = lhs.re + rhs.re;
            result.im = lhs.im + rhs.im;
            return result;
        }

        public static Complex operator -(Complex lhs, Complex rhs)
        {
            Complex result = new Complex();
            result.re = lhs.re - rhs.re;
            result.im = lhs.im - rhs.im;
            return result;
        }

        public static Complex operator *(Complex lhs, Complex rhs)
        {
            Complex result = new Complex();
            result.re = lhs.re * rhs.re;
            result.im = lhs.im * rhs.im;
            return result;
        }

        public static Complex operator *(float lhs, Complex rhs)
        {
            Complex result = new Complex();
            result.re = lhs * rhs.re;
            result.im = lhs * rhs.im;
            return result;
        }

        public static Complex operator *(Complex lhs, float rhs)
        {
            Complex result = new Complex();
            result.re = lhs.re * rhs;
            result.im = lhs.im * rhs;
            return result;
        }
    }
}

请问下我写的IFFT,为什么没还原成原始数据呢(re=0,1,2,3,4,5,6,7)?

  • 写回答

2条回答 默认 最新

  • 皮皮宽 2021-03-16 17:13
    关注

    赚个外块真不容易,看了一下午,

    你的问题如下:

    1,乘法错了, 复数的乘法,你写错了;(a + b j) * (c + d j) = ac-bd + (bc + ad) j 

    2.  你蝶形运算错了,你在好好对照蝶形图捋一捋,比如,第0级对应的旋转因子都是0,你的0123,第1级是0202,你的0202,第2级0123,你的0000,刚好反了,你的分组和每组的数,有问题,你再好好捋一捋吧

    附上修改后的程序和运行结果:

     

    using System;
    
    namespace ConsoleApp34
    {
        class Program
        {
            static void Main(string[] args)
            {
                const int N = 8;
                //初始化时域数据
                Complex[] TD2FD = new Complex[N];
                for (int i = 0; i < N; i++)
                {
                    Complex cpx = new Complex();
                    cpx.re = i;
                    cpx.im = 0;
                    TD2FD[i] = cpx;
                }
    
                Console.WriteLine("------原始信号-----");
                Print(TD2FD);
    
                Console.WriteLine("----- FFT -----");
                FFT(TD2FD);
                Print(TD2FD);
    
                Console.WriteLine("----- IFFT -----");
                IFFT(TD2FD);
                Print(TD2FD);
    
                Console.Read();
            }
    
            //快速傅里叶变换
            public static void FFT(Complex[] TD2FD)
            {
                FFT_Core(TD2FD, WT_LUT(TD2FD.Length, 1));
            }
    
            //快速傅里叶逆变换
            public static void IFFT(Complex[] FD2TD)
            {
                //做FFT变换
                Complex[] WT = WT_LUT(FD2TD.Length, -1);
                FFT_Core(FD2TD, WT);
                //实部除以N
                for (int i = 0; i < FD2TD.Length; i++)
                    FD2TD[i].re /= FD2TD.Length;
            }
    
            // 返回旋转因子查询表(twiddle factor lookup table)
            private static Complex[] WT_LUT(int N, int flag = 1)
            {
                Complex[] WT = new Complex[N];
                for (int i = 0; i < N; i++)
                {
                    Complex cpx_wt = new Complex();
                    float angle = (float)(-i * Math.PI * 2 / N);
                    cpx_wt.re = (float)Math.Cos(angle);
                    //IFFT flag=-1, FFT flag=1
                    cpx_wt.im = flag * (float)Math.Sin(angle);
                    WT[i] = cpx_wt;
                }
                return WT;
            }
    
            private static void FFT_Core(Complex[] TD2FD, Complex[] WT)
            {
                int power = (int)Math.Log(TD2FD.Length, 2);
                int butterfly;
                Complex x1, x2, wt;
                BitReverse(TD2FD);
    
                //蝶形运算
                for (int k = 0; k < power; k++) //级数
                {
                    //Console.WriteLine("第{0}级, 共{1}组", k, 1 << k);
                    for (int j = 0; j < 1 << (power - k - 1); j++) //组数
                    {
                        butterfly = 1 << k + 1;//每组有几个元素
                        //计算参与蝶形运算的两个元素的索引
                        for (int i = 0; i < butterfly / 2; i++) //蝶形运算次数
                        {
                            //Console.Write("  ({0}x{1} wtIdx={2})", i + p, i + s, i * (1 << k));
                            x1 = TD2FD[i + j * butterfly];
                            x2 = TD2FD[i + butterfly/2 + j * butterfly];
                            wt = WT[i * TD2FD.Length / butterfly];
                            TD2FD[i + j * butterfly] = x1 + x2 * wt;
                            TD2FD[i + butterfly / 2 + j * butterfly] = x1 - x2 * wt;
                        }
                        //Console.WriteLine();
                    }
                }
            }
    
            private static int BitReverse(int x)
            {
                //倒位排序
                //0   1   2   3   4   5   6   7   十进制
                //000 001 010 011 100 101 110 111 十进制对应的二进制
                //000 100 010 110 001 101 011 111 码位反转
                //0   4   2   6   1   5   3   7   码位反转后对应的十进制
                int[] table = new int[8] { 0, 4, 2, 6, 1, 5, 3, 7 };
                return table[x];
            }
    
            // 倒位排序——雷德算法
            private static void BitReverse(Complex[] array)
            {
                int i, j, k;
                int N = array.Length;
                Complex temp;
                j = 0;
    
                for (i = 0; i < N - 1; i++)
                {
                    if (i < j)
                    {
                        temp = array[i];
                        array[i] = array[j];
                        array[j] = temp;
                    }
    
                    // 求j的下一个倒序位
                    // N/2的二进制最高位为1,其他位都为0
                    // 判断最高位是否为1,可与N/2进行比较
                    // 判断次高位是否为1,可与N/4进行比较
                    k = N >> 1;
    
                    //如果k<=j,表示j的最高位为1
                    while (k <= j)
                    {
                        //当k<=j时,需要将最高位变为0
                        j = j - k;
                        //判断次高位是否为1,依次类推,逐个比较,直到j某个位为0
                        k >>= 1;
                    }
    
                    j = j + k;//将0变为1
                }
            }
    
            // 打印
            private static void Print(Complex[] TD2FD)
            {
                for (int i = 0; i < TD2FD.Length; i++)
                {
                    Console.WriteLine("(re={0}, im={1})", TD2FD[i].re, TD2FD[i].im);
                }
                Console.WriteLine();
            }
        }
    
        //定义复数
        public class Complex
        {
            public float re;//实数部
            public float im;//虚数部
    
            public static Complex operator +(Complex lhs, Complex rhs)
            {
                Complex result = new Complex();
                result.re = lhs.re + rhs.re;
                result.im = lhs.im + rhs.im;
                return result;
            }
    
            public static Complex operator -(Complex lhs, Complex rhs)
            {
                Complex result = new Complex();
                result.re = lhs.re - rhs.re;
                result.im = lhs.im - rhs.im;
                return result;
            }
    
            public static Complex operator *(Complex lhs, Complex rhs)
            {
                Complex result = new Complex();
                result.re = lhs.re * rhs.re - lhs.im * rhs.im;
                result.im = lhs.re * rhs.im + lhs.im * rhs.re;
                return result;
            }
    
            public static Complex operator *(float lhs, Complex rhs)
            {
                Complex result = new Complex();
                result.re = lhs * rhs.re;
                result.im = lhs * rhs.im;
                return result;
            }
    
            public static Complex operator *(Complex lhs, float rhs)
            {
                Complex result = new Complex();
                result.re = lhs.re * rhs;
                result.im = lhs.im * rhs;
                return result;
            }
        }
    }

    ok,存在可以忽略的误差,很正常,这个没毛病;

    本回答被题主选为最佳回答 , 对您是否有帮助呢?
    评论
查看更多回答(1条)

报告相同问题?

悬赏问题

  • ¥15 我需要全国每个城市的最新小区名字等数据。
  • ¥15 开发一个小区生态的小程序
  • ¥15 MddBootstrapInitialize2失败
  • ¥15 LCD Flicker
  • ¥15 Spring MVC项目,访问不到相应的控制器方法
  • ¥15 esp32在micropython环境下使用ssl/tls连接mqtt服务器出现以下报错Connected on 192.168.154.223发生意外错误: 5无法连接到 MQTT 代理,如何解决?
  • ¥15 关于#genesiscsheel#的问题,如何解决?
  • ¥15 Android aidl for hal
  • ¥15 STM32CubeIDE下载程序报错
  • ¥15 微信好友如何转变为会员系统?(相关搜索:小程序)