- Forward transformation:
$$X_{k} = \sum_{n = 0} ^{N-1} x_n e^{ -2i \pi k n /N }, ~\quad k=0,1,..., N-1;$$ - Backward transformation:
$$x_n = \frac{1}{N}\sum_{k = 0} ^{N-1} X_k e^{ + 2i \pi k n /N }, ~\quad n=0,1,..., N-1;$$ - Butterworth Filter:
$$H(x, y) = \frac{1}{1+\left( \frac{x^2 +y^2}{H_c^2 }\right)^N} ,\quad N=\text{order},~H_c =\text{cut-off};$$ - 주어진 수가 $2^n$인가?(kipl.tistory.com/84). 입력 데이터 개수가 2의 거듭제곱이 아니면 부족한 부분은 0으로 채워서 $2^n \ge \text{데이터 개수}$ 크기가 되도록 만든다.
/* 1차원 FFT:: dir=1 for forward, -1 for backward;
** x = real-component, and y = imaginary componet.
** nn = length of data: 2의 거듭제곱이어야 한다.
*/
int fft1d(int dir, int nn, double x[], double y[]) {
/* Calculate the number of points */
int m = 0;
while (nn > 1) {
nn >>= 1; m++;
}
nn = 1 << m;
/* Do the bit reversal */
int i2 = nn >> 1;
int j = 0;
for (int i = 0; i < nn - 1; i++) {
if (i < j) { // swap(x[i], x[j]), swap(y[i], y[j]);
double tx = x[i], ty = y[i];
x[i] = x[j]; y[i] = y[j];
x[j] = tx; y[j] = ty;
}
int k = i2;
while (k <= j) {
j -= k;
k >>= 1;
}
j += k;
}
/* Compute the FFT */
double c1 = -1.0;
double c2 = 0.0;
int l2 = 1;
for (l = 0; l < m; l++) {
int l1 = l2;
l2 <<= 1;
double u1 = 1.0;
double u2 = 0.0;
for (int j = 0; j < l1; j++) {
for (int i = j; i < nn; i += l2) {
int i1 = i + l1;
double t1 = u1 * x[i1] - u2 * y[i1];
double t2 = u1 * y[i1] + u2 * x[i1];
x[i1] = x[i] - t1;
y[i1] = y[i] - t2;
x[i] += t1;
y[i] += t2;
}
double z = u1 * c1 - u2 * c2;
u2 = u1 * c2 + u2 * c1;
u1 = z;
}
c2 = sqrt((1.0 - c1) / 2.0);
if (dir == 1) c2 = -c2;
c1 = sqrt((1.0 + c1) / 2.0);
}
/* Scaling for forward transform */
if (dir == 1) {
for (i = 0; i < nn; i++) {
x[i] /= double(nn);
y[i] /= double(nn);
}
}
return 1;
};/* 2차원 FFT. data = (2 * nn) x mm 크기의 버퍼;
** 각 행은 nn개의 실수 성분 다음에 nn개의 허수 성분이 온다.
** nn, mm은 각각 행과 열은 나타내며, 2의 거듭제곱이어야 한다.
** isign = 1 for forward, -1 for backward.
** copy = buffer of length 2*mm ;
*/
void fft2d(double* data, int nn, int mm, int isign, double* copy) {
int stride = nn << 1; //row_stride;
/* Transform by ROWS for forward transform */
if (isign == 1) {
int index1 = 0;
for (int i = 0; i < mm; i++) {
// real = &data[index1], imag = &data[index1 + nn];
fft1d(isign, nn, &data[index1], &data[index1 + nn]);
index1 += stride;
}
}
/* Transform by COLUMNS */
for (int j = 0; j < nn; j++) {
/* Copy pixels into temp array */
int index1 = j ;
int index2 = 0;
for (i = 0; i < mm; i++) {
copy[index2] = data[index1];
copy[index2 + mm] = data[index1 + nn];
index2++;
index1 += stride;
}
/* Perform transform */
fft1d(isign, mm, ©[0], ©[mm]);
/* Copy pixels back into data array */
index1 = j;
index2 = 0;
for (i = 0; i < mm; i++) {
data[index1] = copy[index2];
data[index1 + mm] = copy[index2 + mm];
index2++;
index1 += stride;
}
}
/* Transform by ROWS for inverse transform */
if (isign == -1) {
int index1 = 0;
for (i = 0; i < mm; i++) {
fft1d(isign, nn, &data[index1], &data[index1 + nn]);
index1 += stride;
}
}
}void butterworth( double * data, int Xdim, int Ydim,
int Homomorph, int LowPass,
double Power, double CutOff, double Boost);
더보기
// H * F;
void butterworth(double * data, int Xdim, int Ydim,
int Homomorph, int LowPass,
double Power, double CutOff, double Boost)
{
double CutOff2 = CutOff * CutOff;
/* Prepare to filter */
int stride = Xdim << 1; // width of a sigle row(real + imag);
int halfx = Xdim >> 1;
int halfy = Ydim >> 1;
double *rdata = &data[0]; //real part;
double *idata = &data[Xdim];//imag part;
/* Loop over Y axis */
for (int y = 0; y < Ydim; y++) {
int y1 = (y < halfy) ? y: y - Ydim;
/* Loop over X axis */
for (int x = 0; x < Xdim; x++) {
int x1 = (x < halfx) ? x: x - Xdim;
/* Calculate value of Butterworth filter */
double filter;
if (LowPass)
Filter = (1 / (1 + pow((x1 * x1 + y1 * y1) / CutOff2, Power)));
else if ((x1 != 0) || (y1 != 0))
Filter = (1 / (1 + pow( CutOff2 / (x1 * x1 + y1 * y1), Power)));
else
Filter = 0.0;
if (Homomorph)
Filter = Boost + (1 - Boost) * Filter;
/* Do pointwise multiplication */
rdata[x] *= Filter;
idata[x] *= Filter;
};
rdata += stride; //go to next-line;
idata += stride;
}
};
Butterworth Low Pass filter 적용의 예(order=2, cutoff=20)
PowerSpectrum 변화
728x90
'Image Recognition > Fundamental' 카테고리의 다른 글
| Bright Preserving Histogram Equalization with Maximum Entropy (0) | 2008.07.31 |
|---|---|
| Adaptive Binarization (2) | 2008.07.14 |
| Histogram Equalization (0) | 2008.06.22 |
| Otsu Algorithm (6) | 2008.05.30 |
| Hough Transform (2) | 2008.05.22 |
