Fast Hartley Transform (FHT) is known as a real substitute of Fast Fourier Transform (FFT). My problems are described below:

I have an 640*480 image and an Gaussian Kernel of 81*81. I want to do filtering by doing convolution in frequency domain.

I got an one dimension 1024 points FHT routines, which was tested to be correct, and the result is not normalized by 2^n. So for image processing, I need to pass the image in two steps one for column transform, and then rows.

For my case, I think I need to pad my image to 1024*1024. And I also heard somebody saying that I need to transpose both of my image and kernel before doing transform.

My questions are that:

1.Why do I need to transpose my image matrix? Just because I want to place the zero frenquency in the centre of the transformed image?

2. When shall I pad my image? After the transposition or padding first and then transpose?

3. After the inverse transform, the result is 1024*1024, but I need a resulting image of original 640*480. How can I get the result by correctly indexing into the 1024*1024 matrix?

My code doing convolution is listed as following:

float* filter(int* input_img, int width, int height, int fiter_size){

Scale_space ss = new Scale_space();
int max_d = 1024;
int N = max_d*max_d;
double *temp_img = new double [N];
assert(temp_img != NULL)

// load the image
for (int i = 0; i < width; i++){
for (int j = 0; j < height; j++){
temp_img[j*max_d+i] = input_img[j*width+i];
}
}

// Transposition
// swap left right
for(int l=0;l < width/2;++l){
for(int m=0;m < height;++m) {
int temp = temp_img[m*max_d+l];
temp_img[m*max_d+l] = temp_img[m*max_d+width/2+l];
temp_img[m*max_d+l+width/2] = temp;
}
}

// swap upside down
for(int m=0;m < height/2;++m){
for(int l=0;l< width;++l) {
int temp = temp_img[m*max_d+l];
temp_img[m*max_d+l]=temp_img[(m+height/2)*max_d+l];
temp_img[(m+height/2)*max_d+l]=temp;
}
}

// padding image
for(int i=0;i<max_d;i++){
for(int j=0;j<max_d;j++){
if (i<width){
if (j<height){
temp_img[j*max_d+i] = temp_img[j*max_d+i];
}
else
// pad zero
temp_img[j*max_d+i] = 0;
}
else
temp_img[j*max_d+i] = 0;
}
}

// forward transform the image
ss->transform(temp_img);

double* Gx = new double[N];
assert(Gx !=NULL);

int s = 10;
// size of the kernel
kernel_size = 8*s+1;
for (int i = 0; i < kernel_size; i++){
for (int j = 0; j < kernel_size;j++){
//p = i - kernel_size/2;
//q = j - kernel_size/2;
p = i - ((i>=kernel_size/2)?kernel_size:0);
q = j - ((j>=kernel_size/2)?kernel_size:0);
e = exp(-(p*p+q*q)/(2.0*s*s))/(2*PI*s*s);
Gx[j*max_d+i] = e;
}
}

// padding zero
for(int i=0;i<max_d;i++){
for(int j=0;j<max_d;j++){
if (i<kernel_size){
if (j<kernel_size)
Gx[j*max_d+i] = Gx[j*max_d+i];
else
Gx[j*max_d+i] = 0;
}
else
Gx[j*max_d+i] = 0;
}
}

// Transform the mask
ss->transform(Gx);

// point-wise multiplication
for(int l=0;l<N;++l) Gx[l]=Gx[l]*temp_img[l];

// inverse transform
ss->transform(Gx);

float* result = new float [width*height];

// get the useful result
for(int l=0;l<width;++l){
for(int m=0;m<height;++m) {
int lp = (l + 4*smax + (width/2)+4*s) % (width/2+4*s);
int mp = (m + 4*smax + (height/2)+4*s) % (height/2+ 4*s);
int p = mp*max_d+lp;
result[m*width+l] = Gx[p];
}
}
return result;
}
-------------------------------------------------------------------

Am I doing things in correct order?

Please help, urgent for Project, thanks a lot!