2018-10-29 16:16:19 hhaowang 阅读数 8069

目录


1、图像变换
2、离散傅里叶变换(Discrete Fourier Transform)
3、DFT性质
4、DFT与数字图像处理
5、FFT-快速傅里叶变换
6、DFT与FFT的算法实现


1. 图像变换

— —数学领域中有很多种变换,如傅里叶变换、拉普拉斯变换、Z变换等。但是这些数学变换的用途大致相同,即利用某种变换使得遇到的问题更方便的观测或解决。数字图像处理技术(Digital Image Processing)是一门应用学科,它是建立在一定的数学理论基础之上的科学,因此,在解决数字图像处理的具体问题时,作为解决方案或技巧,就必须要用到各种数学变换,称为图像变换。

– 在数字图像处理中,图像增强、图像复原、图像编码和压缩、图像分析与描述等每一种处理手段和方法都要用到图像变换,例如:在对图像进行去噪或锐化操作时,low-pass filter、high-pass filter 可以借助傅里叶变换把要在空间域中解决的卷积运算问题转换并映射到频率域中的乘积运算。

2. 离散傅里叶变换(DFT)

回顾:连续信号的傅里叶变换
在这里插入图片描述
在这里插入图片描述

2.1 一维离散傅里叶变换

为了在计算机上实现傅里叶变换计算,必须把连续函数离散化,离散函数的额傅里叶变换称为离散傅里叶变换。

连续信号离散化过程如下:

  • 建立连续时间信号x(t)
  • 设定采样间隔△t
  • 对连续信号做时域等间隔采样,得到离散序列x(k△t), (k=0,1,2,N-1)
  • 可得到有N个元素的离散函数序列x(n)

在这里插入图片描述

因此,一维离散傅里叶变换(DFT)和一维离散傅里叶反变换
(Inverse DFT)可表示为:
在这里插入图片描述
一维离散傅里叶变换具有周期性,证明如下:
在这里插入图片描述
欧拉方程(Euler theorem):
在这里插入图片描述

2.2 二维离散傅里叶变换

对于一个具有MxN尺度的二维离散函数f(x,y), (x=0,1,2,M-1; y=0,1,2,N-1),其离散傅里叶变换为:
在这里插入图片描述
在运用计算机进行二维离散傅里叶变换时,可先对每一行或者每一列进行一维的离散傅里叶变换,然后对上一步的计算结果再进行列或者行的一维离散傅里叶变换,即可得到二维的离散傅里叶变换结果。具体步骤如下图所示:
在这里插入图片描述

3. DFT Properties

二维离散傅里叶变换具有许多数字图像处理中非常实用的性质,主要包括:可分离性平移性线性特性比例特性周期性和共轭对称性微分特性旋转特性等。
3-1. 可分离性
由二维傅里叶变换不难看出可以把变换方程写成如下形式:
在这里插入图片描述
这个性质称为二维傅里叶变换的可分离性,利用这个性质,可以通过进行两次一维离散傅里叶变换(或反变换)来实现一个二维傅里叶变换(反变换)。以正变换为例,具体算法步骤如下:
在这里插入图片描述
3-2. 平移性
在这里插入图片描述
这一性质表明,当用exp[j2π(ux/M+vy/N)]乘以f(x,y),求乘积的傅里叶变换,可以使得空间频率域u-v平面坐标原点从(0,0)平移到(u,v)的位置。同样,当用exp[-j2π(ux/M+vy/N)]乘以F(x,y),并求此乘积的离散傅里叶反变换,可以使空间x-y平面坐标系的原点从(0,0)平移到(x,y)的位置。

3-3. 线性特性

如果二维离散函数f1(x,y)和f2(x,y)的傅里叶变换分别为F1(x,y)和F1(x,y),则存在以下线性特性。
f1(x,y)+ f2(x,y)<=>F1(x,y)+ F1(x,y)
这一性质可使得我们节约求傅里叶变换的时间,如果已经得到了F1(x,y)和F1(x,y),则f1(x,y)+ f2(x,y)的傅里叶变换结果只需将F1(x,y)和F1(x,y)相加求和即可。

需要注意的是,如果在显示器上显示|F1(x,y)+ F1(x,y)|的结果,则往往由于其最大值超过了数字图像处理设备所允许的最大灰度造成信息丢失现象。因此,在显示|F1(x,y)+ F1(x,y)|的结果时应先检查其最大值,若超过了设备所允许的最大灰度值,则应当先做线性对比度压缩处理。

3-4. 比例特性
在这里插入图片描述

3-5. 周期性和共轭对称性

从二维离散傅里叶变换对的表达式可以看出,离散傅里叶变换和它的反变换具有周期性,即:
在这里插入图片描述

3- 6. 旋转特性

二维傅里叶变换具有如下旋转特性:

f(r,a+b) <=> F(w,φ+φ0)
该性质表明,如果f(x,y)在空间域旋转一个角度a,则对应的傅里叶变换也在频率域旋转同样的角度。反之,亦然。

3-7. 微分特性
在这里插入图片描述
在这里插入图片描述
在模式识别技术中,经常用到Laplace算子。

3- 8. 平均值特性
在这里插入图片描述

4. DFT与数字图像处理

接下来用Lena.jpg标准图来做二维傅里叶变换:
Mac OS 下使用OpenCV和C++集成开发环境做傅里叶变换
代码如下:

//
//  dft_t.cpp
//  OpenCV_demo1
//
//  Created by Hao Wang on 2018/11/2.
//  Copyright © 2018年 Hao Wang. All rights reserved.
//

#include <iostream>
#include <opencv2/core/core.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/opencv.hpp>

using namespace std;
using namespace cv;


int main()
{
    Mat Im_read = imread("/Users/haowang/Desktop/Workspace_C++/Opencv_demo/OpenCV_demo1/Standard Pictures/Lena.jpg", IMREAD_GRAYSCALE);       //读入图像灰度图

    //判断图像是否加载成功
    if (Im_read.empty())
    {
        cout << "图像加载失败!" << endl;
        return -1;    // read error,return -1
    }
    else
        cout << "图像加载成功!"  << endl;

    Mat padded;                 //以0填充输入图像矩阵
    int m = getOptimalDFTSize(Im_read.rows);
    int n = getOptimalDFTSize(Im_read.cols);

    //填充输入图像I,输入矩阵为padded,上方和左方不做填充处理
    copyMakeBorder(Im_read, padded, 0, m - Im_read.rows, 0, n - Im_read.cols, BORDER_CONSTANT, Scalar::all(0));

    Mat planes[] = { Mat_<float>(padded), Mat::zeros(padded.size(),CV_32F) };
    Mat complexI;
    merge(planes, 2, complexI);     //将planes融合合并成一个多通道数组complexI

    dft(complexI, complexI);        //进行傅里叶变换

    //计算幅值,转换到对数尺度(logarithmic scale)
    //=> log(1 + sqrt(Re(DFT(I))^2 + Im(DFT(I))^2))
    split(complexI, planes);        //planes[0] = Re(DFT(I),planes[1] = Im(DFT(I))
    //即planes[0]为实部,planes[1]为虚部
    magnitude(planes[0], planes[1], planes[0]);     //planes[0] = magnitude
    Mat magI = planes[0];

    magI += Scalar::all(1);
    log(magI, magI);                //转换到对数尺度(logarithmic scale)

    //如果有奇数行或列,则对频谱进行裁剪
    magI = magI(Rect(0, 0, magI.cols&-2, magI.rows&-2));

    //重新排列傅里叶图像中的象限,使得原点位于图像中心
    int cx = magI.cols / 2;
    int cy = magI.rows / 2;

    Mat q0(magI, Rect(0, 0, cx, cy));       //左上角图像划定ROI区域
    Mat q1(magI, Rect(cx, 0, cx, cy));      //右上角图像
    Mat q2(magI, Rect(0, cy, cx, cy));      //左下角图像
    Mat q3(magI, Rect(cx, cy, cx, cy));     //右下角图像

    //变换左上角和右下角象限
    Mat tmp;
    q0.copyTo(tmp);
    q3.copyTo(q0);
    tmp.copyTo(q3);

    //变换右上角和左下角象限
    q1.copyTo(tmp);
    q2.copyTo(q1);
    tmp.copyTo(q2);

    //归一化处理,用0-1之间的浮点数将矩阵变换为可视的图像格式
    normalize(magI, magI, 0, 1, CV_MINMAX);

    imshow("输入图像", Im_read);
    imshow("频谱图", magI);
    waitKey(0);
    return 0;

}

计算后得到Lena图的傅里叶变换结果如下:
在这里插入图片描述

5. FFT算法实现

5.1 一维离散傅里叶变换MATLAB实现

%FFT Discrete Fourier transform.
%   FFT(X) is the discrete Fourier transform (DFT) of vector X.  For
%   matrices, the FFT operation is applied to each column. For N-D
%   arrays, the FFT operation operates on the first non-singleton
%   dimension.
%
%   FFT(X,N) is the N-point FFT, padded with zeros if X has less
%   than N points and truncated if it has more.
%
%   FFT(X,[],DIM) or FFT(X,N,DIM) applies the FFT operation across the
%   dimension DIM.
%   
%   For length N input vector x, the DFT is a length N vector X,
%   with elements
%                    N
%      X(k) =       sum  x(n)*exp(-j*2*pi*(k-1)*(n-1)/N), 1 <= k <= N.
%                   n=1
%   The inverse DFT (computed by IFFT) is given by
%                    N
%      x(n) = (1/N) sum  X(k)*exp( j*2*pi*(k-1)*(n-1)/N), 1 <= n <= N.
%                   k=1
%
%   See also FFT2, FFTN, FFTSHIFT, FFTW, IFFT, IFFT2, IFFTN.

%   Copyright 1984-2005 The MathWorks, Inc.

%   Built-in function.

5.2 二维离散傅里叶变换MATLAB实现

function f = fft2(x, mrows, ncols)
%FFT2 Two-dimensional discrete Fourier Transform.
%   FFT2(X) returns the two-dimensional Fourier transform of matrix X.
%   If X is a vector, the result will have the same orientation.
%
%   FFT2(X,MROWS,NCOLS) pads matrix X with zeros to size MROWS-by-NCOLS
%   before transforming.
%
%   Class support for input X: 
%      float: double, single
%
%   See also FFT, FFTN, FFTSHIFT, FFTW, IFFT, IFFT2, IFFTN.

%   Copyright 1984-2010 The MathWorks, Inc. 

if ismatrix(x)
    if nargin==1
        f = fftn(x);
    else
        f = fftn(x,[mrows ncols]);
    end
else
    if nargin==1
        f = fft(fft(x,[],2),[],1);
    else
        f = fft(fft(x,ncols,2),mrows,1);
    end
end   

5.3 OpenCV(C++)实现

#include <cxcore.h>
#include <cv.h>
#include <highgui.h>
 
// Rearrange the quadrants of Fourier image so that the origin is at
// the image center
// src & dst arrays of equal size & type
void cvShiftDFT(CvArr * src_arr, CvArr * dst_arr )
{
    CvMat * tmp;
    CvMat q1stub, q2stub;
    CvMat q3stub, q4stub;
    CvMat d1stub, d2stub;
    CvMat d3stub, d4stub;
    CvMat * q1, * q2, * q3, * q4;
    CvMat * d1, * d2, * d3, * d4;
 
    CvSize size = cvGetSize(src_arr);
    CvSize dst_size = cvGetSize(dst_arr);
    int cx, cy;
 
    if(dst_size.width != size.width || 
       dst_size.height != size.height){
        cvError( CV_StsUnmatchedSizes, "cvShiftDFT", "Source and Destination arrays must have equal sizes", __FILE__, __LINE__ );   
    }
 
    if(src_arr==dst_arr){
        tmp = cvCreateMat(size.height/2, size.width/2, cvGetElemType(src_arr));
    }
 
    cx = size.width/2;
    cy = size.height/2; // image center
 
    q1 = cvGetSubRect( src_arr, &q1stub, cvRect(0,0,cx, cy) );
    q2 = cvGetSubRect( src_arr, &q2stub, cvRect(cx,0,cx,cy) );
    q3 = cvGetSubRect( src_arr, &q3stub, cvRect(cx,cy,cx,cy) );
    q4 = cvGetSubRect( src_arr, &q4stub, cvRect(0,cy,cx,cy) );
    d1 = cvGetSubRect( dst_arr, &d1stub, cvRect(0,0,cx,cy) );
    d2 = cvGetSubRect( dst_arr, &d2stub, cvRect(cx,0,cx,cy) );
    d3 = cvGetSubRect( dst_arr, &d3stub, cvRect(cx,cy,cx,cy) );
    d4 = cvGetSubRect( dst_arr, &d4stub, cvRect(0,cy,cx,cy) );
 
    if(src_arr!=dst_arr){
        if( !CV_ARE_TYPES_EQ( q1, d1 )){
            cvError( CV_StsUnmatchedFormats, "cvShiftDFT", "Source and Destination arrays must have the same format", __FILE__, __LINE__ ); 
        }
        cvCopy(q3, d1, 0);
        cvCopy(q4, d2, 0);
        cvCopy(q1, d3, 0);
        cvCopy(q2, d4, 0);
    }
    else{
        cvCopy(q3, tmp, 0);
        cvCopy(q1, q3, 0);
        cvCopy(tmp, q1, 0);
        cvCopy(q4, tmp, 0);
        cvCopy(q2, q4, 0);
        cvCopy(tmp, q2, 0);
    }
}
 
int main(int argc, char ** argv)
{
    const char* filename = argc >=2 ? argv[1] : "lena.jpg";
    IplImage * im;
 
    IplImage * realInput;
    IplImage * imaginaryInput;
    IplImage * complexInput;
    int dft_M, dft_N;
    CvMat* dft_A, tmp;
    IplImage * image_Re;
    IplImage * image_Im;
    double m, M;
 
    im = cvLoadImage( filename, CV_LOAD_IMAGE_GRAYSCALE );
    if( !im )
        return -1;
 
    realInput = cvCreateImage( cvGetSize(im), IPL_DEPTH_64F, 1);
    imaginaryInput = cvCreateImage( cvGetSize(im), IPL_DEPTH_64F, 1);
    complexInput = cvCreateImage( cvGetSize(im), IPL_DEPTH_64F, 2);
 
    cvScale(im, realInput, 1.0, 0.0);
    cvZero(imaginaryInput);
    cvMerge(realInput, imaginaryInput, NULL, NULL, complexInput);
 
    dft_M = cvGetOptimalDFTSize( im->height - 1 );
    dft_N = cvGetOptimalDFTSize( im->width - 1 );
 
    dft_A = cvCreateMat( dft_M, dft_N, CV_64FC2 );
    image_Re = cvCreateImage( cvSize(dft_N, dft_M), IPL_DEPTH_64F, 1);
    image_Im = cvCreateImage( cvSize(dft_N, dft_M), IPL_DEPTH_64F, 1);
 
    // copy A to dft_A and pad dft_A with zeros
    cvGetSubRect( dft_A, &tmp, cvRect(0,0, im->width, im->height));
    cvCopy( complexInput, &tmp, NULL );
    if( dft_A->cols > im->width )
    {
        cvGetSubRect( dft_A, &tmp, cvRect(im->width,0, dft_A->cols - im->width, im->height));
        cvZero( &tmp );
    }
 
    // no need to pad bottom part of dft_A with zeros because of
    // use nonzero_rows parameter in cvDFT() call below
 
    cvDFT( dft_A, dft_A, CV_DXT_FORWARD, complexInput->height );
 
    cvNamedWindow("win", 0);
    cvNamedWindow("magnitude", 0);
    cvShowImage("win", im);
 
    // Split Fourier in real and imaginary parts
    cvSplit( dft_A, image_Re, image_Im, 0, 0 );
 
    // Compute the magnitude of the spectrum Mag = sqrt(Re^2 + Im^2)
    cvPow( image_Re, image_Re, 2.0);
    cvPow( image_Im, image_Im, 2.0);
    cvAdd( image_Re, image_Im, image_Re, NULL);
    cvPow( image_Re, image_Re, 0.5 );
 
    // Compute log(1 + Mag)
    cvAddS( image_Re, cvScalarAll(1.0), image_Re, NULL ); // 1 + Mag
    cvLog( image_Re, image_Re ); // log(1 + Mag)
 
 
    // Rearrange the quadrants of Fourier image so that the origin is at
    // the image center
    cvShiftDFT( image_Re, image_Re );
 
    cvMinMaxLoc(image_Re, &m, &M, NULL, NULL, NULL);
    cvScale(image_Re, image_Re, 1.0/(M-m), 1.0*(-m)/(M-m));
    cvShowImage("magnitude", image_Re);
 
    cvWaitKey(-1);
    return 0;
}
2017-03-28 21:30:37 penkee 阅读数 2523

《数字图像处理》第三版 Rafael C. Gonzalez等著  P155

花了1,2个月断续的看书,终于有些进展了,一直到DFT这里,偏理论些。


到现在我都不清楚弄出图片的DFT的振幅图像,相位图来干啥的。不过没关系,按照作者写的来,先搞出来再说。


接下要做的是如何用程序画出变换的振幅图

按照公式F(u,v),可以看出来每个点的值是一个复数,假设F(u0,v0)=a + b j;书上定义,振幅图的值是  a的平方+b的平方的和开根号。


e到三角函数用的是欧拉公式




这样公式推导后,程序就非常容易计算了

发现每次计算振幅图的每一个点都需要遍历原图所有的点,所以性能这块需要考虑下


很奇怪,第一次搞这个,感觉和书上差很大,公式没有错啊

原图

     

图像的傅立叶变换振幅图

 


原图


图像的傅立叶变换振幅图



https://github.com/penkee/imagecal/blob/master/app-dao/src/main/java/com/dcloud/process/FourierService.java


按照书上说的,把灰度范围压缩到【0-255】,而不是直接去掉,效果果然一样了。









2016-11-09 12:07:55 BaiYH1994 阅读数 2655

  图像傅里叶变换方法有很多,可以通过空间光调制器输入图像后在通过平行光照明经过傅里叶变换透镜进行傅里叶变换,另一个方法就是利用计算机进行傅里叶变换,其中傅里叶变换有两种算法一种是DFT还有一种是FFT(快速傅里叶变换)。

  首先我介绍一下图像的定义,图像是怎么去得到的呢?图像是物体与点扩散函数卷积的结果加上一个噪声项(具体想了解可以查看“傅里叶光学"-吕乃光)。

  在这里我谈谈离散傅里叶变换的一些问题,以及离散傅里叶变换的MATLAB代码。首先谈谈傅里叶变换,想必大家一听到傅里叶变换整个人都不好了,傅里叶变换是做什么的呢,个人的理解就是通过无数个正交的向量去描述一个任意的曲线或者信号,最简单的理解就是我们可以通过xy直角坐标系通过一个关系建立一个二维曲线,这样我们就可以通过一个正交的向量去描述这个线上的任意一个点。

  在傅里叶变换中所谓的空间频率的高频与低频分别代表什么呢?高频代表的就是图像中灰度突变速度快的地方,说白了就是看起来密集的地方就是高频,看起来平坦的地方就是低频。这是最简单的理解方式。离散傅里叶变换可以处理任意大小的图像,FFT只能处理2^n*2^n大小的图像,当图像大小不够的时候需要填补空缺的地方。

一下为DFT函数的matlab代码,

<span style="font-size:18px;">function [out_dft_image,out_ab_dft_pic]=dft(input_image,mode)
%%  [out_dft_image,out_ab_dft_pic]=dft_test(input_image)     this function can make dft transform
% F(u-M/2,v-N/2)<=>f(x,y)*(-1).^(x+y)   or you can make dft transform like
% F(u,v)<=>f(x,y)
%     input_image:  had batter less than 150*150 large
%     mode:   if mode==1  this function is to do comfortable visable result   
%             else mode=~1 the result maybe make you uncomfortable 
%     out_dft_image: this is the complex matrix it can make the idft and
%                   get the surse image in spetial domain 
%     out_ab_dft_pic: this is the frequency spectrum picture matrix
%     design by baiyinhao 2015.9.28    20:07 Email:792499178@qq.com
%
%     作  者:光电科技协会     2015.9.28    20:07  白银浩
%


%空域变换到频率域
pic=double(input_image);
[u_p,v_p]=size(pic);%  F(u,v)中的u,v
%f(x,y)*(-1).^(x+y)   即对原图进行一次运算
if mode==1
    
    for x=1:u_p
        for y=1:v_p
            pic(x,y)=pic(x,y)*(-1)^(x+y);%在此处对输入的pic进行了改变  即变成f(x,y)*(-1).^(x+y) 白  添加 
        end
    end
    
    %居中傅里叶变换
    [x_p,y_p]=size(pic);
    dft_pic=zeros(u_p,v_p);
    flag=0;
    temp1=0;
    for u=1:u_p
        for v=1:v_p
            for x=1:x_p
                for y=1:y_p
                    temp1=temp1+pic(x,y)*exp(-1j*2*pi*(u*x/u_p+v*y/v_p));
                end
            end
            dft_pic(u,v)=temp1;
            temp1=0;
        end
        flag=flag+1;
        sprintf('运行到了图像的第:%d 行',flag)
    end
    %输出傅里叶变换 后的实部虚部
    R_dft_pic=real(dft_pic);%求出傅里叶变换后的实部
    I_dft_pic=imag(dft_pic);%求出傅里叶变换后的虚部

    %傅里叶谱
    ab_dft_pic(u_p,v_p)=0;
    for u=1:u_p
        for v=1:v_p
            ab_dft_pic(u,v)=sqrt(R_dft_pic(u,v).^2+I_dft_pic(u,v).^2);
       end
    end
    %最后输出值 
    m_dft_pic=max(max(ab_dft_pic));
    ab_dft_pic1=ab_dft_pic/m_dft_pic*255;



else
   %正常傅里叶变换
    [x_p,y_p]=size(pic);
    dft_pic=zeros(u_p,v_p);
    flag=0;
    temp1=0;
    for u=1:u_p
        for v=1:v_p
            for x=1:x_p
                for y=1:y_p
                    temp1=temp1+pic(x,y)*exp(-1j*2*pi*(u*x/u_p+v*y/v_p));
                end
            end
            dft_pic(u,v)=temp1;
            temp1=0;
        end
        flag=flag+1;
        sprintf('运行到了图像的第:%d 行',flag)
    end
    %输出傅里叶变换 后的实部虚部 白银浩添加
    R_dft_pic=real(dft_pic);%求出傅里叶变换后的实部
    I_dft_pic=imag(dft_pic);%求出傅里叶变换后的虚部

    %傅里叶谱
    ab_dft_pic(u_p,v_p)=0;
    for u=1:u_p
        for v=1:v_p
            ab_dft_pic(u,v)=sqrt(R_dft_pic(u,v).^2+I_dft_pic(u,v).^2);
        end
    end
    %最后输出值 
    m_dft_pic=max(max(ab_dft_pic));
    ab_dft_pic1=ab_dft_pic/m_dft_pic*255;
 end


    out_dft_image=dft_pic;
    out_ab_dft_pic=uint8(ab_dft_pic1);


end</span>

由于运行速度慢,有助于理解,如果有实际用途建议使用FFT。


2016.11.9   12:09

白银浩    

E-mail:BaiYH1994@163.com

Q    Q:792499178

2018-11-26 16:52:44 shanwenkang 阅读数 504

酉变换

酉变换可以由如下方式定义,其中输入和输出之间的关系可以写成矩阵相乘的形式,矩阵A称为酉矩阵,A满足A的逆矩阵等于A的共轭对称矩阵

DFT变换就是一个酉变换,系数矩阵A满足每一列的模是1并且由于不同频率正弦信号之间的正交性,列之间是相互正交,因此A也是一个酉矩阵

对于二维DFT我们可以看做两次一维的DFT,因此我们也可以写成矩阵相乘的形式

 

我们表达一个二维图像或者是一个一维向量,我们都是用基的形式来表示

空域基

我们在空域上表示一个图像实际上就是把每个像素看做一个基,然后每个点的像素合起来形成一整幅图像

傅里叶基

我们在频域上表示图像实际上是把图像中不同频率成分的向量看做不同的基来表示图像

如下图就是频域的基在空域上的体现,也就是说每幅图像都是由这些图像加权得到的

二维DFT有它的好处,但是它也有坏处,比如傅里叶变换后得到的频谱值是一个复数,但是我们实际上只需要一个实数,正因为如此我们很多时候不得不用abs来对系数取模;第二个缺点就是每个基实际上都包含了整幅图的信息,这使得有的时候我们表达一幅图像的开销更大。比如我们有一幅图只有一个像素点是1,那么我们用空域的表达只需要一个系数,但是如果我们要在频域上表示,那么我们将用很多个255来表示它(类比冲激响应的频谱是一条直线)

 

其他常用变换

离散余弦变换

为了解决以上第一个问题,我们可以改变图像的基,既然我们的基带有虚数,那我们只需要采用实数基就可以了,这也就是离散余弦变换(DCT)的由来。离散余弦变换也满足酉变换的条件,并且也像DFT一样存在快速算法

如下图所示是离散余弦变换的基在空域上的表示

DCT是JPEG标准非常重要的算法

与离散余弦变换类似的是还存在离散正弦变换

Hadamard变换

以上提到的变换都是浮点型运算,为了使运算中只存在加减,我们引入了Hadamard变换

如下图所示是Hadamard变换的基在空域上的表示

Harr变换

为了解决傅里叶变换的第二个问题,我们引入了Harr变换,Harr变换是一种小波变换,有关小波介绍戳我

如下图所示是小波变换的基在空域上的表示,可以看到小波变换的基不是像其他变换一样覆盖整个图像的,而是在高频的部分将基函数本地化来覆盖图像的不同区域

以下是小波变换的优点,我们可以看出最重要的一点就是它拥有局部的基函数,从而在时域划分与频域划分上找到了一种平衡

 

 

2017-09-16 14:40:58 u011271038 阅读数 1265

这几天一直学习傅立叶变换,看了很多国内外资料,网上讲原理的很多,到了程序实现这块大多是Matlab,opencv等,这些软件的api对于我们理解DFT在计算机中的实现并没有多大帮助。于是想用C/C++实现DFT,经过不断的阅读与编程实验,最终程序有了还算满意的结果。

“关于傅立叶变换,无论是书本还是在网上可以很容易找到关于傅立叶变换的描述,但是大都是些故弄玄虚的文章,太过抽象,尽是一些让人看了就望而生畏的公式的罗列,让人很难能够从感性上得到理解”—dznlong

这是在网上看到的一句话,也是我一直以来的心声,要想弄明白一些事情,终究是要自己动手。

至于原理性的内容,网上太多,我这里不多说了,仅列几个连接供参考。
https://en.wikipedia.org/wiki/Fourier_series
http://www.cnblogs.com/v-July-v/archive/2011/02/20/1983676.html
http://www.dspguide.com/pdfbook.htm

 for (int u = 0; u < width; u ++)
     for (int v =0; v < height; v ++){

         for (int x = 0; x < width; x++){
             for (int y = 0; y < height; y++){

                real += srcMat.at<Vec3b>(y, x)[0] * cos(- ((double)2 * PI * u * x / width + (double)2 * PI * v * y / height )) * 1;
                imaginary += srcMat.at<Vec3b>(y, x)[0] * sin(-((double)2 * PI * u * x / width + (double)2 * PI * v * y /height)) * 1;
            }
        }

u和v是变换后的变量,从程序可以看出,每一个F(u,v)的计算都需要遍历整副图像,确实是相当费时间,我只用了100 × 100 的图像,计算速度慢的惊人。

 real =  100 * real / (width * height);
 imaginary =100 *  imaginary / (width * height);

为什么要乘以100呢?因为此时的real和imaginary值很小,如果显示出来基本是黑色,所以为了显示清楚,我们随便乘以个值放大一下,这应该就是Opencv例子里说的归一化。

“`
tmpMat.at(u, v) = sqrt(real * real + imaginary * imaginary);

“`这里写图片描述
上面那个小方块窗口是原图像变换后的图像,用的是opencv的方法。
这里写图片描述
这个图是用C语言实现的DFT,为了看得清晰点,我把它放大了,相比与opencv的图这个没有做坐标的变换。

这个程序不是完美的,只是一种尝试,通过这个过程,对DFT有进一步了解,养成这种自己动手的习惯后,对学习其他的算法有一定帮助,接下来我该转移到其他问题学习,当知识丰富后还会在回过头来重新审视今天这些内容,并加以完善。

没有更多推荐了,返回首页