고속 퓨리에 변환 프로그램 구현

A Program demonstrating Fast Fourier Transform(FFT) Functions

 

개요

고속 퓨리에 변환은 공간 영역에서의 신호를 주파수 영역으로 고속 변환하여 영상 신호의 대역폭을 줄임으로써 변환 기반 영상 압축을 가능하게 한다.

 

■ 퓨리에 변환식은 실수의 범위을 넘어 복소수 범위까지 확장 (Complex 클래스 구현)

■ 구현은 Visual Studio 2005 & C#

■ FFT는 2차원 DFT를 1차원 DFT로 분할 적용(행/열)

■ FFT는 비트리버셜(입력 데이터 순서 변환)로 반복 분해 연산

■ 고속 퓨리에 변환(FFT)와 역변환(iFFT)

■ FFT 변환된 복소수 데이터로부터 원래의 영상으로 완전히 복원

 

주제	영상변환 : FFT(Fast Fourier Transform, 고속 퓨리에 변환)
관련기술	※ 구현하고자 하는 요소기술에서 사용되는 기술 내용을 간단히 입력한다. 1)퓨리에(Fourier) 변환 기술 2)순방향 FFT와 역방향 FFT 적용 구현 기술
구현내용 및 방법	※ 구현하고자 하는 요소 기술에 대한 계획을 기술한다. 1. 구현 언어 : Visual C# (Microsoft Visual Studio 2005) 2. 구현 내용 : - 영상에 대해 고속 퓨리에 변환(FFT)를 수행하여 그 결과를 주파수 스펙트럼 정보를 갖는 데이터로 변환 - 역방향 FFT 구현 3. 구현 방법 : - 입력 데이터 영상을 불러오기 - 순방향 FFT와 역방향 FFT를 적용한 데이터를“원본 영상/주파수 스펙트럼/재배치된 스펙트럼/순,역방향 적용 영상“으로 구분하여 출력하는 UI제공

 

시스템 구성 및 명세

순방향FFT()	영상에 대해 FFT를 수행하여 그 결과를 주파수 스펙트럼으로 나타내는 메인 루틴
역방향FFT()	역방향 FFT 메인루틴
FFT_2D()	입력 영상 데이터를 갖는 2차원 복소수형 데이터를 FFT에의해 주파수 스펙트럼 정보를 갖는 데이터로 변환하는 2차원 FFT
iFFT_2D()	2차원 복소수형 주파수 데이터를 역방향 FFT에 의해 영상 정보를 갖는 데이터로 변환해 주는 역방향 2차원 FFT
FFT()	분해된 복소수 1차원 배열 데이터를 받아서 주파수 정보로 변환한 다음 그 결과를 반환
iFFT()	분해된 복소수 1차원 배열의 주파수 데이터를 받아서 영상 정보로 변환한 다음 그 결과를 반환하는 루틴
scramble()	재귀적인 DFT 계산의 주기와 대칭을 이루도록 데이터를적절히 재배치
butterflies()	데이터를 점 집합으로 나누고 이웃한점에 대해 DFT 수행하는 나비 흐름 계산
iButterflies()	butterflies()의 역방향
viewFourier()	FFT나 iFFT로 추출한 2차원 복소수 데이터를 화면에 표시할 수 있는 0 ~ 255 값으로 변환
viewCenter()	변환 영상을 4등분 하여 4모서리로 분산된 저주파 영역을 가운데로 대칭이동하여 주파수 스펙트럼을 재배치

 

 

FFT 및 iFFT 변환 결과

 

 

개발환경

프로그램 소스

 

/* 퓨리에변환식은실수의범위를넘어서복소수범위까지확장되어있으므로 
이를위한프로그램을작성하기위해서는복소수를표현할수있어야한다
그런데C#의기본자료형에는복소수가없으므로Complex 클래스를정의한다*/
/// <summary>
/// 퓨리에변환을위한복소수클래스
/// </summary>
public class Complex
{
double re;       // 실수
double im;       // 허수

public double Real
{
get { return re;  }
set { re = value }
}
public double Imaginary
{
get { return im; }
set { im = value}
}
public Complex()
{
}
public Complex(double x, double y)
{
Real      = x;
Imaginary = y;
}

        /// <summary>
        ///두복소수곱셈연산을위한메소드
        /// </summary>
        public static Complex Multiply(Complex a, Complex b)
{
Complex c = new Complex();
c.Real = a.Real*b.Real-a.Imaginary*b.Imaginary;
c.Imaginary = a.Imaginary*b.Real + a.Real*b.Imaginary;
return c;
}
}


/// <summary>
/// 고속퓨리에변환Fourier Transform:FFT)
/// </summary>
public partial class Form1 : Form
{
Graphics grBm;
Image image;
Bitmap bitmap;
int[,] grayArray; //입력영상의픽셀값저장배열

public Form1()
{
InitializeComponent();
setShadowBitmap();
}
void setShadowBitmap()
{
bitmap = new Bitmap(ClientSize.Width, ClientSize.Height);
grBm = Graphics.FromImage(bitmap);
grBm.Clear(BackColor);
}
private void Form1_Paint(object sender, PaintEventArgs e)
{
Graphics gr = e.Graphics;
gr.DrawImage(bitmap, 0, 0);
}
private void openToolStripMenuItem_Click(object sender, EventArgs e)
{
ofdOpen.Title = "Open Image"
ofdOpen.Filter = "JPEG file(*.jpg)|*.jpg|Bitmap File(*.bmp)|*.bmp|GIF file(*.gif)|*.gif|PNG file(*.png)|*.png|TIFF(*.tif)|*.tif|All Files(*.*)|*.*"

if (ofdOpen.ShowDialog() == DialogResult.OK)
{                
string strFilename = ofdOpen.FileName;
image = Image.FromFile(strFilename);
this.ClientSize = new Systehttp://m.Drawing.Size(image.Width, image.Height);
setShadowBitmap();
grBm.DrawImage(image, 0, 0);

Wait wait = new Wait("이미지의픽셀값을저장하고있습니다);
wait.Show();
copyBitmap2Array(); //입력영상을2차원배열에저장
          wait.Close();
}
Invalidate();
}
void copyBitmap2Array()
{
int x, y, brightness;
Color color;

grayArray = new int[bitmap.Height, bitmap.Width];
for (y = 0; y < bitmap.Height; y++)
for (x = 0; x < bitmap.Width; x++)
{
color = bitmap.GetPixel(x, y);
brightness = color.R;
grayArray[y, x] = brightness;
}
}
void displayResultArray(int[,] resultArray, int Width, int Height)
{
int x, y;
Color color;
Bitmap gBitmap = new Bitmap(Width, Height);

for (y = 0; y < Height; y++)
for (x = 0; x < Width; x++)
{
color = Color.FromArgb(resultArray[y, x], resultArray[y, x], resultArray[y, x]);
gBitmap.SetPixel(x, y, color);
}
setShadowBitmap();
grBm.DrawImage(gBitmap, 0, 0, gBitmap.Width, gBitmap.Height);
Invalidate();
}

/// <summary>
/// 퓨리에변환이나역변환으로추출한2차원복소수데이터를화면에표시할수있는
/// 0 ~ 255 값으로변환시키기위한메소드
/// </summary>
int[,] viewFourier(Complex[,] G, int Width, int Height)
{
int[,] R = new int[Height, Width];
double[,] t = new double[Height, Width]; //디스플레이Log 함수
         double[,] m = new double[Height, Width]; //복소수의크기
         int x, y;

for (y = 0; y < Height; y++)
for (x = 0; x < Width; x++)
m[y, x] = Math.Sqrt(G[y, x].Real * G[y, x].Real + G[y, x].Imaginary * G[y, x].Imaginary);

for (y = 0; y < Height; y++)
for (x = 0; x < Width; x++)
t[y, x] = Math.Log10(1 + m[y, x]); //로그의가수가0이되지않기위해1을더한다
            
double dmax = double.MinValue; //디스플레이함수값중가장큰값데이터를0~255로맞추기위해필요한비율을구한다
         for (y = 0; y < Height; y++)
for (x = 0; x < Width; x++)
if (t[y, x] > dmax) dmax = t[y, x];

double dConst = 255.0 / dmax;
for (y = 0; y < Height; y++)
for (x = 0; x < Width; x++)
R[y, x] = (int)(dConst * t[y, x]);
return R;
}

/// <summary>
/// 퓨리에변환에의해디스플레이되는주파수스펙트럼을재배치하는메소드 주파수영역에서는중

/// 심좌표인(Width/2, Height/2)이(0,0)이되어야하므로네모서리점으로분산된것을 중심으로

/// 모아야제대로된주파수스펙트럼을볼수있다 변환된영상을기준으로4등분했을때각각좌상좌

/// 하우상우하부분별로각각의중심원점에대한대칭이동을 수행함으로써네모서리로분산된것을

/// 중심점으로재배치할수있다
/// </summary>
int[,] viewCenter(int[,] G, int Width, int Height)
{
int[,] R = new int[Height, Width];
int x, y;
int xMid = Width / 2;
int yMid = Height / 2;

// 좌상원점대칭이동
         for (y = 0; y < yMid; y++)
for (x = 0; x < xMid; x++)
R[yMid - 1 - y, xMid - 1 - x] = G[y, x];
// 좌하원점대칭이동
         for (y = yMid; y < Height; y++)
for (x = 0; x < xMid; x++)
R[Height - 1 - y + yMid, xMid - 1 - x] = G[y, x];
// 우상원점대칭이동
         for (y = 0; y < yMid; y++)
for (x = xMid; x < Width; x++)
R[yMid - 1 - y, Width - 1 - x + xMid] = G[y, x];
// 우하원점대칭이동
         for (y = yMid; y < Height; y++)
for (x = xMid; x < Width; x++)
R[Height - 1 - y + yMid, Width - 1 - x + xMid] = G[y, x];
return R;
}

/// <summary>
/// 영상에대해FFT를수행하여그결과를주파수스펙트럼으로나타내는메인루틴
/// 1.영상의폭과높이가2의n승인지제약검사
/// 2.입력영상을복소수형의실수부분에들어가게하고허수부분은0.0이되게복소수데이터로변환
/// 3.2차원FFT 수행하여주파수영역데이터로변환
/// 4.주파수스펙트럼을화면에표시
/// 5.네모서리로분산된것을중심으로재배치
/// </summary>
private void 순방향FFTToolStripMenuItem_Click(object sender, EventArgs e)
{
int M = 0, N = 0, tnum;

tnum = image.Width;
while (tnum > 1)
{
tnum >>= 1; M++;
}
if ((tnum << M) != image.Width)
{
MessageBox.Show("영상의폭이2의n승이아닙니다);
return
}

tnum = image.Height;
while (tnum > 1)
{
tnum >>= 1; N++;
}
if ((tnum << N) != image.Height)
{
MessageBox.Show("영상의높이가2의n승이아닙니다);
return
}

Wait wait = new Wait("FFT 변환작업을진행중입니다);
wait.Show();

Complex[,] grayFourArray = new Complex[image.Height, image.Width];
for (int y = 0; y < image.Height; y++)
for (int x = 0; x < image.Width; x++)
grayFourArray[y, x] = new Complex((double)grayArray[y, x], 0.0); //입력영상은복소수형실수부에허수부는0.0이되게복소수형데이터로변환
            
Complex[,] FourierArray = FFT_2D(grayFourArray, image.Width, image.Height, M, N); //2차원FFT수행주파수영역데이터로변환
         int[,] viewArray = viewFourier(FourierArray, image.Width, image.Height); //주파수스펙트럼을화면에표시
         displayResultArray(viewArray, image.Width, image.Height);            
wait.Close();
MessageBox.Show("확인버튼을누르면주파수스펙트럼을재배치합니다,"변환완료);
wait = new Wait("주파수스펙트럼을재배치중입니다);
wait.Show();

int[,] ResultArray = viewCenter(viewArray, image.Width, image.Height);
//네모서리로분산된것을중심으로재배치
         displayResultArray(ResultArray, image.Width, image.Height);
wait.Close();
}

/// <summary>
/// 입력영상데이터를갖는2차원복소수형데이터를FFT에의해주파수스펙트럼정보를갖는데이터로변환하는2차원FFT
/// 먼저행에대해서FFT 수행그결과를다시열에대해적용하여변환수행
/// </summary>
Complex[,] FFT_2D(Complex[,] G, int Width, int Height, int M, int N)
{
Complex[,] R = new Complex[image.Height, image.Width];
Complex[] OneDimColData = new Complex[image.Height]; // 열데이터는각각높이만큼의길이를갖음
         Complex[] OneDimRowData = new Complex[image.Width];  // 행데이터는각각폭만큼의길이를갖음

         //행에대한퓨리에변환
         for (int y = 0; y < Height; y++)
for (int x = 0; x < Width; x++)
OneDimRowData[x] = new Complex();

for (int y = 0; y < Height; y++)
{
for (int x = 0; x < Width; x++)
{
OneDimRowData[x].Real = G[y, x].Real;
OneDimRowData[x].Imaginary = G[y, x].Imaginary;
}
Complex[] FourierRow = FFT(OneDimRowData, M, Width);
for (int x = 0; x < Width; x++)
R[y, x] = new Complex(FourierRow[x].Real, FourierRow[x].Imaginary);
}
//

         //열에대한퓨리에변환
         for (int x = 0; x < Width; x++)
for (int y = 0; y < Height; y++)
OneDimColData[y] = new Complex();

for (int x = 0; x < Width; x++)
{
for (int y = 0; y < Height; y++)
{
OneDimColData[y].Real = R[y, x].Real;
OneDimColData[y].Imaginary = R[y, x].Imaginary;
}
Complex[] FourierCol = FFT(OneDimColData, N, Height);
for (int y = 0; y < Height; y++)
{
R[y, x].Real = FourierCol[y].Real;
R[y, x].Imaginary = FourierCol[y].Imaginary;
}
}
//
return R;
}

/// <summary>
/// 분해된복소수1차원배열데이터를받아서주파수정보로변환한다음그결과를반환하는루틴
/// scramble(), butterflies()
/// </summary>
Complex[] FFT(Complex[] G, int log2N, int length)
{
Complex[] scrambleG = scramble(G, length);
Complex[] R = butterflies(scrambleG, log2N, length);
return R;
}

/// <summary>
/// 재귀적인DFT 계산의주기와대칭을이루도록데이터를적절히재배치
/// 역순의인덱스비트를배열의성분으로계산하며데이터는역순에의해지정되는것과교환
/// </summary>
Complex[] scramble(Complex[] G, int length)
{
int i, j, m;
double temp;
Complex[] R = new Complex[length];

for (i = 0; i < length; i++)
R[i] = new Complex(G[i].Real, G[i].Imaginary);

j = 0;

for (i = 0; i < length; i++)
{
if (i > j)
{
temp = R[j].Real;
R[j].Real = R[i].Real;
R[i].Real = temp;
temp = R[j].Imaginary;
R[j].Imaginary = R[i].Imaginary;
R[i].Imaginary = temp;
}
m = length >> 1;
while ((j >= m) & (m > 1))
{
j -= m; m >>= 1;
}
j += m;
}
return R;
}

/// <summary>
/// 나비흐름계산
/// 데이터를점집합으로나누고이웃한점에대해DFT 수행
/// </summary>
Complex[] butterflies(Complex[] G, int log2N, int length)
{
int i, j, k, offset;
int n, halfN;
double theta, tdbl;
Complex exp = new Complex();
Complex iw = new Complex();
Complex w = new Complex();
Complex[] R = new Complex[length];

for (i = 0; i < length; i++)
R[i] = new Complex(G[i].Real, G[i].Imaginary);

n = 1;
for (i = 0; i < log2N; i++)
{
halfN = n;
n <<= 1;
theta = -2.0 * Math.PI / (double)n;
tdbl = Math.Sin(0.5 * theta);
exp.Real = -2.0 * tdbl * tdbl;
exp.Imaginary = Math.Sin(theta);
w.Real = 1.0; w.Imaginary = 0.0;
for (offset = 0; offset < halfN; offset++)
{
for (k = offset; k < length; k += n)
{
j = k + halfN;
iw = Complex.Multiply(w, R[j]);
R[j].Real = R[k].Real - iw.Real;
R[k].Real += iw.Real;
R[j].Imaginary = R[k].Imaginary - iw.Imaginary;
R[k].Imaginary += iw.Imaginary;
}
tdbl = w.Real;
w.Real = tdbl * exp.Real - w.Imaginary * exp.Imaginary + w.Real;
w.Imaginary = w.Imaginary * exp.Real + tdbl * exp.Imaginary + w.Imaginary;
}
}

for (i = 0; i < length; i++)
{
R[i].Real /= (double)length;
R[i].Imaginary /= (double)length;
}
return R;
}

/// <summary>
/// 역방향FFT 메인루틴
/// 1.영상의폭과높이가2의n승인지제약검사
/// 2.입력영상을복소수형의실수부분에들어가게하고허수부분은0.0이되게복소수데이터로변환       

/// 3.2차원FFT 수행하여주파수영역데이터로변환
/// 4.곧바로역방향FFT 수행
/// 5.화면에표시      
/// </summary>
private void 역방향FFTToolStripMenuItem_Click(object sender, EventArgs e)
{
int M = 0, N = 0, tnum;

tnum = image.Width;
while (tnum > 1)
{
tnum >>= 1; M++;
}
if ((tnum << M) != image.Width)
{
MessageBox.Show("영상의폭이2의n승이아닙니다);
return
}           

tnum = image.Height;
while (tnum > 1)
{
tnum >>= 1; N++;
}
if ((tnum << N) != image.Height)
{
MessageBox.Show("영상의높이가2의n승이아닙니다);
return
}

Wait wait = new Wait("FFT 역변환작업을진행중입니다);
wait.Show();

Complex[,] grayFourArray = new Complex[image.Height, image.Width];
for (int y = 0; y < image.Height; y++)
for (int x = 0; x < image.Width; x++)
   grayFourArray[y, x] = new Complex((double)grayArray[y, x], 0.0);

   Complex[,] FourierArray = FFT_2D(grayFourArray, image.Width, image.Height, M, N);
   Complex[,] iFourierArray = iFFT_2D(FourierArray, image.Width, image.Height, M, N);
int[,] viewArray = viewFourier(iFourierArray, image.Width, image.Height);
displayResultArray(viewArray, image.Width, image.Height);
wait.Close();
}

/// <summary>
/// 2차원복소수형주파수데이터를역방향FFT에의해영상정보를갖는데이터로변환해주는역방향2
/// 차원FFT먼저행에대해서역방향FFT 수행그결과를다시열에대해적용하여변환수행
/// </summary>
Complex[,] iFFT_2D(Complex[,] G, int Width, int Height, int M, int N)
{
Complex[,] R = new Complex[image.Height, image.Width];
Complex[] OneDimColData = new Complex[image.Height]; // 열데이터는각각높이만큼의길이를갖음
         Complex[] OneDimRowData = new Complex[image.Width];  // 행데이터는각각폭만큼의길이를갖음

         for (int y = 0; y < Height; y++)
for (int x = 0; x < Width; x++)
OneDimRowData[x] = new Complex();

for (int y = 0; y < Height; y++)
{
for (int x = 0; x < Width; x++)
{
OneDimRowData[x].Real = G[y, x].Real;
OneDimRowData[x].Imaginary = G[y, x].Imaginary;
}
Complex[] FourierRow = iFFT(OneDimRowData, M, Width);
for (int x = 0; x < Width; x++)
R[y, x] = new Complex(FourierRow[x].Real, FourierRow[x].Imaginary);
}

for (int x = 0; x < Width; x++)
for (int y = 0; y < Height; y++)
OneDimColData[y] = new Complex();

for (int x = 0; x < Width; x++)
{
for (int y = 0; y < Height; y++)
{
OneDimColData[y].Real = R[y, x].Real;
OneDimColData[y].Imaginary = R[y, x].Imaginary;
}
Complex[] FourierCol = iFFT(OneDimColData, N, Height);
for (int y = 0; y < Height; y++)
{
R[y, x].Real = FourierCol[y].Real;
R[y, x].Imaginary = FourierCol[y].Imaginary;
}
}
return R;
}

/// <summary>
/// 분해된복소수1차원배열의주파수데이터를받아서영상정보로변환한다음그결과를반환
/// scramble(), ibutterflies()
/// </summary>
Complex[] iFFT(Complex[] G, int log2N, int length)
{
Complex[] scrambleG = scramble(G, length);
Complex[] R = iButterflies(scrambleG, log2N, length);
return R;
}

/// <summary>
/// 역방향FFT 나비흐름계산
/// </summary>
Complex[] iButterflies(Complex[] G, int log2N, int length)
{
int i, j, k, offset;
int n, halfN;
double theta, tdbl;
Complex exp = new Complex();
Complex iw = new Complex();
Complex w = new Complex();
Complex[] R = new Complex[length];

for (i = 0; i < length; i++)
R[i] = new Complex(G[i].Real, G[i].Imaginary);

n = 1;
for (i = 0; i < log2N; i++)
{
halfN = n;
n <<= 1;
theta = 2.0 * Math.PI / (double)n; //theta의부호만순방향FFT와다르다
                  tdbl = Math.Sin(0.5 * theta);
exp.Real = -2.0 * tdbl * tdbl;
exp.Imaginary = Math.Sin(theta);
w.Real = 1.0; w.Imaginary = 0.0;
for (offset = 0; offset < halfN; offset++)
{
for (k = offset; k < length; k += n)
{
j = k + halfN;
iw = Complex.Multiply(w, R[j]);
R[j].Real = R[k].Real - iw.Real;
R[k].Real += iw.Real;
R[j].Imaginary = R[k].Imaginary - iw.Imaginary;
R[k].Imaginary += iw.Imaginary;
}
tdbl = w.Real;
w.Real = tdbl * exp.Real - w.Imaginary * exp.Imaginary + w.Real;
w.Imaginary = w.Imaginary * exp.Real + tdbl * exp.Imaginary + w.Imaginary;
}
}

for (i = 0; i < length; i++)
{
R[i].Real /= (double)length;
R[i].Imaginary /= (double)length;
}
return R;
}
}


[참고문헌 및 사이트]
[1]정민영, 이칠우, “C# 디지털 영상처리” 2005년 9월 도서출판 미래컴
[2]정성태, “Visual C++를 이용한 실용 영상 처리” 2007년 2월 생능출판사
[3]http://mecha.tu.ac.kr/common/download.asp?bcode=pds&bseq=413 

[4]http://mbm.konkuk.ac.kr/%EB%94%94%EC%A7%80%ED%84%B8%EC%98%81%EC%83%81%EC%B2%98%EB%A6%AC10.ppt
[5]http://logos.mokwon.ac.kr/lect/DSP/DSP07.ppt
[6]http://home.megapass.co.kr/~londonfo/down/FT.pdf