C# 实现 FFT 正反变换 和 频域滤波

C# 实现 FFT 正反变换 和 频域滤波

要进行FFT运算首先要构造复数类,参考

http://blog.csdn.net/iamoyjj/archive/2009/05/15/4190089.aspx

 

下面的程序在依赖上述复数类的基础上实现了FFT正反变换算法和频域滤波算法,另外由于一般如果是对实数进行FFT的话,要将FFT得到的复数数组转为实数数组,下面类中的Cmp2Mdl方法的作用就是这个。这个FFT算法是基-2FFT算法,因此,如入的序列必须是2n次方个点长。

频域滤波的基本原理是:

1、  对输入序列进行FFT

2、  得到的频谱乘以一个权函数(滤波器,系统的传递函数)

3、  得到的结果进行IFFT

4、  如果是实数运算的话用Cmp2Mdl方法转为实数

代码如下:

    /// <summary>

    /// 频率分析器

    /// </summary>

    public static class FreqAnalyzer

    {

 

        /// <summary>

        /// 求复数complex数组的模modul数组

        /// </summary>

        /// <param name=”input”>复数数组</param>

        /// <returns>模数组</returns>

        public static double[] Cmp2Mdl(Complex[] input)

        {

            ///有输入数组的长度确定输出数组的长度

            double[] output = new double[input.Length];

 

            ///对所有输入复数求模

            for (int i = 0; i < input.Length; i++)

            {

                if (input[i].Real > 0)

                {

                    output[i] = -input[i].ToModul();

                }

                else

                {

                    output[i] = input[i].ToModul();

                }

            }

 

            ///返回模数组

            return output;

        }

 

        /// <summary>

        /// 傅立叶变换或反变换,递归实现多级蝶形运算

        /// 作为反变换输出需要再除以序列的长度

        /// !注意:输入此类的序列长度必须是2^n

        /// </summary>

        /// <param name=”input”>输入实序列</param>

        /// <param name=”invert”>false=正变换,true=反变换</param>

        /// <returns>傅立叶变换或反变换后的序列</returns>

        public static Complex[] FFT(double[] input, bool invert)

        {

            ///由输入序列确定输出序列的长度

            Complex[] output = new Complex[input.Length];

 

            ///将输入的实数转为复数

            for (int i = 0; i < input.Length; i++)

            {

                output[i] = new Complex(input[i]);

            }

 

            ///返回FFTIFFT后的序列

            return output = FFT(output, invert);

        }

 

        /// <summary>

        /// 傅立叶变换或反变换,递归实现多级蝶形运算

        /// 作为反变换输出需要再除以序列的长度

        /// !注意:输入此类的序列长度必须是2^n

        /// </summary>

        /// <param name=”input”>复数输入序列</param>

        /// <param name=”invert”>false=正变换,true=反变换</param>

        /// <returns>傅立叶变换或反变换后的序列</returns>

        private static Complex[] FFT(Complex[] input, bool invert)

        {

            ///输入序列只有一个元素,输出这个元素并返回

            if (input.Length == 1)

            {

                return new Complex[] { input[0] };

            }

 

            ///输入序列的长度

            int length = input.Length;

 

            ///输入序列的长度的一半

            int half = length / 2;

 

            ///有输入序列的长度确定输出序列的长度

            Complex[] output = new Complex[length];

 

            ///正变换旋转因子的基数

            double fac = -2.0 * Math.PI / length;

 

            ///反变换旋转因子的基数是正变换的相反数

            if (invert)

            {

                fac = -fac;

            }

 

            ///序列中下标为偶数的点

            Complex[] evens = new Complex[half];

 

            for (int i = 0; i < half; i++)

            {

                evens[i] = input[2 * i];

            }

 

            ///求偶数点FFTIFFT的结果,递归实现多级蝶形运算

            Complex[] evenResult = FFT(evens, invert);

 

            ///序列中下标为奇数的点

            Complex[] odds = new Complex[half];

 

            for (int i = 0; i < half; i++)

            {

                odds[i] = input[2 * i + 1];

            }

 

            ///求偶数点FFTIFFT的结果,递归实现多级蝶形运算

            Complex[] oddResult = FFT(odds, invert);

 

            for (int k = 0; k < half; k++)

            {

                ///旋转因子

                double fack = fac * k;

 

                ///进行蝶形运算

                Complex oddPart = oddResult[k] * new Complex(Math.Cos(fack), Math.Sin(fack));

                output[k] = evenResult[k] + oddPart;

                output[k + half] = evenResult[k] – oddPart;

            }

 

            ///返回FFTIFFT的结果

            return output;

        }

 

        /// <summary>

        /// 频域滤波器

        /// </summary>

        /// <param name=”data”>待滤波的数据</param>

        /// <param name=”Nc”>滤波器截止波数</param>

        /// <param name=”Hd”>滤波器的权函数</param>

        /// <returns>滤波后的数据</returns>

        private static double[] FreqFilter(double[] data, int Nc, Complex[] Hd)

        {

            ///最高波数==数据长度

            int N = data.Length;

 

            ///输入数据进行FFT

            Complex[] fData = FreqAnalyzer.FFT(data, false);

 

            ///频域滤波

            for (int i = 0; i < N; i++)

            {

                fData[i] = Hd[i] * fData[i];

            }

 

            ///滤波后数据计算IFFT

            ///!未除以数据长度

            fData = FreqAnalyzer.FFT(fData, true);

 

            ///复数转成模

            data = FreqAnalyzer.Cmp2Mdl(fData);

 

            ///除以数据长度

            for (int i = 0; i < N; i++)

            {

                data[i] = -data[i] / N;

            }

 

            ///返回滤波后的数据

            return data;

        }

    }

 

版权声明:本文内容由互联网用户自发贡献,该文观点仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 举报,一经查实,本站将立刻删除。

发布者:全栈程序员-用户IM,转载请注明出处:https://javaforall.cn/110913.html原文链接:https://javaforall.cn

【正版授权,激活自己账号】: Jetbrains全家桶Ide使用,1年售后保障,每天仅需1毛

【官方授权 正版激活】: 官方授权 正版激活 支持Jetbrains家族下所有IDE 使用个人JB账号...

(0)


相关推荐

发表回复

您的电子邮箱地址不会被公开。

关注全栈程序员社区公众号