`

图像的时频变换——离散傅里叶变换

阅读更多

 

一直很纳闷,几乎所有数字图像处理的书都会介绍数字时频变换,但是却很少书会讲时频变换的作用,这也是让我一直很疑惑的地方(不过也极有可能是本人愚钝)。频谱技术通常用于提高图像的处理操作速率,频谱相当于是图像的特征,时频变换虽然是一种数学技巧,但是运用到数字图像处理上会方便和简单。研究的图像变换基本上都是正交变换,正交变换可以减少图像数据的相关性,获取图像的整体特点,有利于用较少的数据量表示原始图像,这对图像的分析、存储以及图像的传输都是非常有意义的这里介绍了离散傅立叶变换、离散余弦变换、沃尔什-哈达玛变换及小波变换的基本理论和知识并进行图像时频变换实验。

一、二维离散傅里叶变换

一个图像尺寸为M×N的 函数f(x,y)的离散傅里叶变换由以下等式给出:


 其中 u=0,1,2,...,M-1和v=0,1,2,...,N-1。其中变量u和v用于确定它们的频率,频域系统是由F(u,v)所张成的坐标系,其中u和v用做(频率)变量。空间域是由f(x,y)所张成的坐标系。可以得到频谱系统在频谱图四角(0,0),(0,N-1),(N-1,0),(N-1,N-1)处沿u和v方向的频谱分量均为0。

离散傅里叶逆变换由下式给出:

令R和I分别表示F的实部和需部,则傅里叶频谱,相位角,功率谱(幅度)定义如下:



 在频谱的原点变换值称为傅里叶变换的直流分量,下面是傅里叶变换的周期公式:

 
 DFT实现仅计算一个周期,所以可以处理M X N的数组,由周期性可知,A,B,C,D是四个分别的四分之一周期。

  图像的频率是表征图像中灰度变化剧烈程度的指标,是灰度在平面空间上的梯。傅立叶频谱图上我们看到的明暗不一的亮点,实际上图像上某一点与邻域点差异的强弱,即梯度的大小,也即该点的频率的大小(可以这么理解,图像中的低频部分指低梯度的点,高频部分相反)。一般来讲,梯度大则该点的亮度强,否则该点亮度弱。这样通过观察傅立叶变换后的频谱图,也叫功率图,我们首先就可以看出,图像的能量分布,如果频谱图中暗的点数更多,那么实际图像是比较柔和的(因为各点与邻域差异都不大,梯度相对较小),反之,如果频谱图中亮的点数多,那么实际图像一定是尖锐的,边界分明且边界两边像素差异较大的。对频谱移频到原点以后,可以看出图像的频率分布是以原点为圆心,对称分布的。变换最慢的频率成分(u=v=0)对应一幅图像的平均灰度级。当从变换的原点移开时,低频对应着图像的慢变换分量,较高的频率开始对应图像中变化越来越快的灰度级。这些事物体的边缘和由灰度级的突发改变(如噪声)标志的图像成分。通常在进行傅里叶变换之前用(-1)^(x+y)乘以输入的图像函数,这样就可以将傅里叶变换的原点F(0,0)移到(M/2,N/2)上。

 

二、算法原理

 

离散快速傅里叶变换(FFT)是在离散傅立叶变换(DFT)的基础上改进的,为了更好的理解FFT算法,先对DFT进行梳理下,

1、DFT运算的特点:

首先分析有限长序列x(n)进行一次DFT运算所需的运算量。

一般,x(n)和都是复数,因此,每计算一个X(k)值,要进行N次复数相乘,和N-1次复数相加,X(k)一共有N个点,故完成全部DFT运算,需要N^2次复数相乘和N(N-1)次复数相加,在这些运算中,乘法比加法复杂,需要的运算时间多,尤其是复数相乘,每个复数相乘包括4个实数相乘和2个实数相加,例

又每个复数相加包括2个实数相加,所以,每计算一个X(k)要进行4^N次实数相乘和2N+2(N-1)=2(2N-1)次实数相加,因此,整个DFT运算需要4N^2实数相乘和2N(2N-1)次实数相加。

2、FFT算法的基本思想:

考察DFT与IDFT的运算发现,利用以下两个特性可减少运算量:

1)系数 是一个周期函数,它的周期性和对称性可用来改进运算,提高计算效率。

 

又如因此因此

利用这些周期性和对称性,使DFT运算中有些项可合并;

2)利用的 周期性和对称性,把长度为N点的大点数的DFT运算依次分解为若干个小点数的DFT。因为DFT的计算量正比于N二次方,N小,计算量也就小。

FFT算法正是基于这样的基本思想发展起来的。它有多种形式,但基本上可分为两类:时间抽取法和频率抽取法。

3、按时间抽取的FFT(N点DFT运算的分解)

假设采样序列点数为N=2^L,L为整数,如果不满足这个条件可以人为地添加若干个0以使采样序列点数满足这一要求。首先我们将序列x(n)按照奇偶分为两组如下: 

将DFT运算也相应分为两组:


注意到,H(k),G(k)有N/2个点,即k=0,1,…,N/2-1,还必须应用系数 的周期性和对称性表示X(k)的N/2 ~N-1点:

可见,一个N点的DFT被分解为两个N/2点的DFT,这两个N/2点的DFT再合成为一个N点DFT

 依此类推,G(k)和H(k)可以继续分下去,这种按时间抽取算法是在输入序列分成越来越小的子序列上执行DFT运算,最后再合成为N点的DFT。

蝶形信号流图

将G(k)和H(k) 合成X(k)运算可归结为:

蝶形运算的简化

图(a)为实现这一运算的一般方法,它需要两次乘法、两次加减法。考虑到-bW和bW两个乘法仅相差一负号,可将图(a)简化成图2.7(b),此时仅需一次乘法、两次加减法。图(b)的运算结构像一蝴蝶通常称作蝶形运算结构简称蝶形结,采用这种表示法,就可以将以上所讨论的分解过程用流图表示。

N=2^3=8的例子。

 按照这个办法,继续把N/2用2除,由于N=2M,仍然是偶数,可以被2整除,因此可以对两个N/2点的DFT再分别作进一步的分解。即对{G(k)}和{H(k)}的计算,又可以分别通过计算两个长度为N/4=2点的DFT,进一步节省计算量,见图。这样,一个8点的DFT就可以分解为四个2点的DFT。



 最后剩下的是2点DFT,它可以用一个蝶形结表示:


 

这样,一个8点的完整的按时间抽取运算的流图

由于这种方法每一步分解都是按输入时间序列是属于偶数还是奇数来抽取的,所以称为“按时间抽取法”或“时间抽取法”。

 4、按频率抽取的FFT(按输出X(k)在频域的顺序上属于偶数还是奇数分解为两组)

对于N=2^M情况下的另外一种普遍使用的FFT结构是频率抽取法。

对于频率抽取法,输入序列不是按偶数奇数,而是按前后对半分开,这样便将N点DFT写成前后两部分:


 
 令

这两个序列都是N/2点的序列,将其代入上两式,得


 
 这正是两个N/2点的DFT运算,即将一个N点的DFT分解为两个N/2点的DFT,上式的运算关系可用下图表示.

以N=8的频率抽取为例



  按频率抽取将8点DFT分解成四个2点DFT

与时间抽取法一样,由于N=2M,N/2仍是一个偶数,这样,一个N=2M点的DFT通过M次分解后,最后只剩下全部是2点的DFT,2点DFT实际上只有加减运算。但为了比较,也为了统一运算的结构,仍用一个系数为W0N的蝶形运算来表示。

频率抽取法的流图是时域抽取法流图的左右翻转。下图是N=8的频率抽取法FFT流图。

 N=8的按频率抽取FFT运算流图

5、实数序列的FFT

以上讨论的FFT算法都是复数运算,包括序列x(n)也认为是复数,但大多数场合,信号是实数序列,任何实数都可看成虚部为零的复数,例如,求某实信号x(n)的复谱,可认为是将实信号加上数值为零的虚部变成复信号(x(n)+j0),再用FFT求其离散傅里叶变换。这种作法很不经济,因为把实序列变成复序列,存储器要增加一倍,且计算机运行时,即使虚部为零,也要进行涉及虚部的运算,浪费了运算量。合理的解决方法是利用复数据FFT对实数据进行有效计算,下面介绍两种方法。

(1)用一个N点FFT同时计算两个N点实序列的DFT设x(n)、y(n)是彼此独立的两个N点实序列,且X(k)=DFT[x(n)],Y(k)=DFT[y(n)],则X(k)、Y(k)可通过一次FFT运算同时获得。

首先将x(n)、y(n)分别当作一复序列的实部及虚部,令g(n)=x(n)+jy(n)

通过FFT运算可获得g(n)的DFT值

 

利用离散傅里叶变换的共轭对称性


 

通过g(n)的FFT运算结果G(k),由上式可得到X(k)的值。

通过g(n)的FFT运算结果G(k),由上式也可得

到Y (k) 的值。

作一次N点复序列的FFT,再通过加、减法运算就可以

将X(k)与Y(k)分离出来。显然,这将使运算效率提高一倍。

6、用FFT计算二维离散的傅里叶变换

前面介绍的FFT原理都是关于一维傅里叶变换的,终于又能切回正题了,数字图像傅里叶变换明显是一个二维的离散傅里叶变换,二维信号有图象信号、时空信号、时频信号等。二维离散傅里叶变换可用于处理二维离散信号。

二维离散傅里叶变换的定义为:

 二维离散傅里叶变换可通过两次一维离散傅里叶变换来实现:

1)作一维N点DFT(对每个m做一次,共M次)

2)作M点的DFT(对每个k做一次,共N次)


这两次离散傅里叶变换都可以用快速算法求得,若M和N都是2的幂,则可使用基二FFT算法,所需要乘法次数为

而直接计算二维离散傅里叶变换所需的乘法次数为(M+N)MN,当M和N比较大时用用FFT运算,可节约很多运算量。

 

三、基于数字图像处理FFT编程实现

 

数字图像处理的快速傅里叶变换操作一共有下面3步

1)获取图片的像素的灰度数组

 

double[] iPix = new double[iw * ih];//输入灰度

               pixels = common.grabber(iImage, iw, ih); //获得像素数组

        for (int i = 0; i < iw*ih; i++)   //获得像素灰度值

                 iPix[i] = pixels[i]&0xff;   

 

2)对灰度数组进行傅里叶变换

 

               //FFT变换

               FFT2 fft2 = new FFT2();

               fft2.setData2(iw, ih, iPix);

               fftData = fft2.getFFT2();

首先是将实数转换成复数,然后根据欧拉公式  e^ix=cosx+isinx 计算蝶形因子(乘数),下面是一维FFT的核心程序,程序分为三层循环
  第一层循环:由于N=2^m需要m级计算,第一层循环对运算的级数进行控制。

   第二层循环:由于第L级有2^(L-1)个蝶形因子(乘数),第二层循环根据乘数进行控制,保证对于每一个蝶形因子第三层循环要执行一次,这样,第三层循环在第二层循环控制下,每一级要进行2^(L-1)次循环计算。(L=1,2.......m)

   第三层循环:由于第L级共有N/2^L个群,并且同一级内不同群的乘数分布相同,当第二层循环确定某一乘数后,第三层循环要将本级中每个群中具有这一乘数的蝶形计算一次,即第三层循环每执行完一次要进行N/2^L个碟形计算。

   可以得出结论:在每一级中,第三层循环完成N/2^L个碟形计算;第二层循环使第三层循环进行 2^(L-1)次,因此,第二层循环完成时,共进行2^(L-1) *N/2^L=N/2个碟形计算。实质是:第二、第三层循环完成了第L级的计算。

  几个要注意的数据:

   ① 在第L级中,每个碟形的两个输入端相距b=2^(L-1)个点(程序中的bfsize/2)。

   ② 同一乘数对应着相邻间隔为2^L=2*b个点的N/2^L个碟形。

   ③ L级的2^(L-1)个碟形因子W N[p]中的P,可表示为p = j*2^(m-L)

       其中j = 012...2^(L-1)-1


private Complex[] getData1()
    {
        //蝶形运算
        for (k = 0; k < power; k++)
        {
            for (j = 0; j < 1 << k; j++)
            {
                bfsize = 1 << (power - k);
                for (i = 0; i < bfsize / 2; i++)
                {
                    Complex temp1 = new Complex(0, 0);
                    Complex temp2 = new Complex(0, 0);

                    p = j * bfsize;
                    x2[i + p] = temp1.Add(x1[i + p], x1[i + p + bfsize / 2]);

                    temp2 = temp1.Sub(x1[i + p], x1[i + p + bfsize / 2]);

                    x2[i + p + bfsize / 2] = temp1.Mul(temp2, wc[i * (1 << k)]);                
                }
            }
            x = x1;
            x1 = x2;
            x2 = x;
        }

        //重新排序
        for (j = 0; j < count; j++)
        {
            p = 0;
            for (i = 0; i < power; i++)
                if ((j & (1 << i)) != 0)
                    p += 1 << (power - i - 1);

            fd1[j] = x1[p];
        }
        return fd1;
    }
 

 

3)FFT数据的可视化

使用傅里叶频谱来可视化FFT变换得到的数据

 

//FFT数据可视化
    public int[] toPix(Complex[] fftData, int iw, int ih)
    {
        int[] pix = new int[iw*ih];
	    
        int u, v, r;
        for (int j = 0; j < ih; j++)
        {
            for (int i = 0; i < iw; i++)
            {
                double tem = fftData[i + j * iw].re * fftData[i + j * iw].re
                           + fftData[i + j * iw].im * fftData[i + j * iw].im;
                r = (int)(Math.sqrt(tem) / 100);
                if (r > 255) r = 255;
                //移到中心
                if (i < iw / 2) u = i + iw / 2;
                else u = i - iw / 2;
                if (j < ih / 2) v = j + ih / 2;
                else v = j - ih / 2;
                
                pix[i + j * iw] = 255 << 24|r << 16|r << 8|r;
            }
        }
        return pix; 
    }         

  实验效果如下:


 

 

 

 

 

 

 

 

 

 

 

  • 描述: 二维离散傅里叶变换公式
  • 大小: 5.5 KB
  • 大小: 5.5 KB
  • 大小: 8.2 KB
  • 大小: 16.5 KB
  • 大小: 17.6 KB
  • 大小: 15.3 KB
  • 大小: 26.3 KB
  • 大小: 14.1 KB
  • 大小: 1.8 KB
  • 大小: 19 KB
  • 大小: 4.1 KB
  • 大小: 6.4 KB
  • 大小: 64.6 KB
  • 大小: 31 KB
  • 大小: 3 KB
  • 大小: 14 KB
  • 大小: 4.9 KB
  • 大小: 68.1 KB
  • 大小: 57.2 KB
  • 大小: 26.3 KB
  • 大小: 29.2 KB
  • 大小: 23.2 KB
  • 大小: 64.6 KB
  • 大小: 89.6 KB
  • 大小: 24.1 KB
  • 大小: 98.3 KB
  • 大小: 37.6 KB
  • 大小: 75 KB
  • 大小: 13.4 KB
  • 大小: 15 KB
  • 大小: 12 KB
  • 大小: 42.4 KB
  • 大小: 82.3 KB
  • 大小: 45.8 KB
  • 大小: 22.3 KB
  • 大小: 46.8 KB
  • 大小: 17 KB
  • 大小: 41.1 KB
  • 大小: 14.9 KB
  • 大小: 12.6 KB
  • 大小: 13.6 KB
  • 大小: 12.4 KB
  • 大小: 54.8 KB
  • 大小: 12.6 KB
分享到:
评论
1 楼 bao674785731 2016-04-01  
pix[i + j * iw] = 255 << 24|r << 16|r << 8|r;
能请教一下这个是什么意思吗

相关推荐

    傅立叶变换的原理、意义和应用 (2).docx

    傅立叶变换分为连续傅里叶变换(Continuous Fourier Transform)和离散傅里叶变换(Discrete Fourier Transform,简称DFT)。连续傅里叶变换用于分析连续的模拟信号,而离散傅里叶变换则是处理数字信号的关键。由于...

    jpeg中数据的时频转换.rar_jpeg_离散余弦变换

    标题中的“jpeg中数据的时频转换.rar_jpeg_离散余弦变换”指的是JPEG(Joint Photographic Experts Group)图像压缩标准中的一项关键技术——离散余弦变换(DCT,Discrete Cosine Transform)。JPEG是一种广泛使用的...

    基于的系统分析与设计——时频分析.rar

    时频分析方法如短时傅里叶变换(Short-Time Fourier Transform, STFT)、小波变换(Wavelet Transform)以及 chirplet 变换等,能够提供更为精细的时频分布图,帮助我们更好地解析信号的动态特征。 在MATLAB中,...

    【数字信号处理】数字信号离散序列短时傅里叶变换含Matlab源码.zip

    《数字信号处理——离散序列短时傅里叶变换与Matlab实现》 在现代通信、音频处理和图像分析等领域,数字信号处理是一项至关重要的技术。离散序列短时傅里叶变换(Short-Time Fourier Transform, STFT)是其中一种...

    数字信号处理 小波变换 课件

    1. 时频分析的重要性:解释了在处理非平稳信号时,传统傅立叶变换的局限性。 2. 小波函数的定义和性质:介绍小波函数的基本结构,如正交性、平移不变性和缩放特性。 3. 小波变换的类型:包括连续小波变换和离散小波...

    专题讲座——小波变换PPT学习教案.pptx

    时频分析包括了多种变换,如【短时傅里叶变换 (STFT)】、【Gabor变换】和【小波变换 (WT)】。 【短时傅里叶变换】通过在不同时间点对信号加窗然后做傅里叶变换,来获取信号的局部频率特性。但是,STFT的分辨率受到...

    小波与傅里叶分析基础课件

    常见的傅里叶变换有连续傅里叶变换和离散傅里叶变换(DFT),其中离散傅里叶变换的快速算法——快速傅里叶变换(FFT)在实际应用中尤为重要,因为它大大降低了计算复杂度。 **小波分析概述** 小波分析是对傅里叶...

    小波变换详解

    与传统的傅立叶变换相比,小波变换能够提供更为精细的时频局部化信息,这主要得益于其灵活的窗口大小调整机制。 ### 小波变换的起源与发展 小波变换的理论基础可追溯至傅立叶变换,由法国科学家Joseph Fourier于...

    tft.zip_fft_simulink 信号_傅里叶 simulink_时域频域转换_时频转换

    FFT是一种高效算法,用于计算离散傅里叶变换(DFT)和其逆变换。在信号处理中,FFT是将时域信号转化为频域信号的关键工具。时域信号是随时间变化的连续或离散函数,而频域信号则表示信号包含的不同频率成分及其幅度...

    基于FPGA和matlab的图像压缩算法的设计与实现

    接着,论文详细分析了四种不同的图像压缩算法:傅里叶变换、离散余弦变换(DCT)、小波变换和沃尔什——哈达玛变换,并探讨了它们在图像分析中的应用及各自的不足和性能分析。 傅里叶变换是一种将信号从时域转换到...

    The Fractional Fourier Transform and Its Applications

    分数傅里叶变换可以通过多种数学手段实现,比如基于离散傅里叶变换(DFT)的实现方法。由于分数傅里叶变换和傅里叶变换在本质上有相似之处,它们在理论上可能的应用范围也是相似的。这表明,分数傅里叶变换可以应用...

    《数字信号处理教程——MATLAB释义与实现》

    2. **傅里叶变换**:讲解傅里叶级数和傅里叶变换的理论,包括离散傅里叶变换(DFT)和快速傅里叶变换(FFT)。通过MATLAB实现信号的频域分析。 3. **滤波器设计**:介绍各种滤波器类型,如低通、高通、带通和带阻...

    Visual C++小波变换技术与工程实践pdf

    小波变换是信号处理和图像分析中的一个重要工具,它结合了时频分析的优点,能对非平稳信号进行精细分析。Visual C++作为Microsoft开发的面向对象的编程环境,为实现小波变换提供了强大的支持。 小波变换的基本概念...

    信号与系统分析实验(基于MATLAB)

    在MATLAB中,虽然无法直接计算DTFT,但可以通过DFT(离散傅里叶变换)和窗函数来近似。实验会讲解DTFT的性质,如周期性和复共轭对称性,并探讨采样率和窗函数对结果的影响。 在这些实验中,你将学习如何使用MATLAB...

    VC++,C++

    相比于傅立叶变换,小波变换能够提供更为精确的时频信息,特别适合于处理非平稳信号。 在VC++或C++中实现小波变换,通常会涉及到以下关键步骤和知识点: 1. **小波基的选择**:常见的小波基有Haar小波、Daubechies...

    WT.rar_小波变换讲义

    小波变换相较于传统的傅立叶变换,其独特之处在于具有时频局部化特性,即能够同时在时间和频率域上提供精细的信息。这使得小波变换在处理非平稳信号时表现出色。 **小波的基本概念** 小波(Wavelet)是一种具有...

Global site tag (gtag.js) - Google Analytics