FFTW图像处理之频域滤波和频域分析

1. 频域滤波实例

低通滤波（模糊图像/去噪）

c

void lowpass_filter(const char* input_path, const char* output_path, int cutoff_radius) {int width, height;unsigned char *img = stbi_load(input_path, &width, &height, NULL, 1); // 灰度图像fftw_complex *in = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * width * height);fftw_complex *out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * width * height);// 准备输入数据for (int y = 0; y < height; y++) {for (int x = 0; x < width; x++) {in[y*width + x][0] = img[y*width + x] / 255.0;in[y*width + x][1] = 0;}}// 正向FFTfftw_plan plan_forward = fftw_plan_dft_2d(height, width, in, out, FFTW_FORWARD, FFTW_ESTIMATE);fftw_execute(plan_forward);// 频域中心坐标int center_x = width / 2;int center_y = height / 2;// 低通滤波：保留中心区域for (int y = 0; y < height; y++) {for (int x = 0; x < width; x++) {int dx = x - center_x;int dy = y - center_y;if (dx*dx + dy*dy > cutoff_radius*cutoff_radius) {out[y*width + x][0] = 0;out[y*width + x][1] = 0;}}}// 逆FFTfftw_plan plan_backward = fftw_plan_dft_2d(height, width, out, in, FFTW_BACKWARD, FFTW_ESTIMATE);fftw_execute(plan_backward);// 保存结果unsigned char *result = (unsigned char*)malloc(width * height);for (int i = 0; i < width*height; i++) {double val = in[i][0] / (width * height); // 归一化result[i] = (unsigned char)(fmin(fmax(val * 255, 0), 255));}stbi_write_png(output_path, width, height, 1, result, width);// 清理资源fftw_destroy_plan(plan_forward);fftw_destroy_plan(plan_backward);fftw_free(in);fftw_free(out);free(result);stbi_image_free(img);
}

高通滤波（边缘增强）

c

void highpass_filter(const char* input_path, const char* output_path, int cutoff_radius) {// ... (前面部分与低通滤波相同，直到频域处理部分)// 高通滤波：去除中心区域for (int y = 0; y < height; y++) {for (int x = 0; x < width; x++) {int dx = x - center_x;int dy = y - center_y;if (dx*dx + dy*dy <= cutoff_radius*cutoff_radius) {out[y*width + x][0] = 0;out[y*width + x][1] = 0;}}}// ... (剩余部分与低通滤波相同)
}

带通滤波（保留特定频率）

c

void bandpass_filter(const char* input_path, const char* output_path, int inner_radius, int outer_radius) {// ... (前面部分相同)// 带通滤波：保留环形区域for (int y = 0; y < height; y++) {for (int x = 0; x < width; x++) {int dx = x - center_x;int dy = y - center_y;int dist_sq = dx*dx + dy*dy;if (dist_sq < inner_radius*inner_radius || dist_sq > outer_radius*outer_radius) {out[y*width + x][0] = 0;out[y*width + x][1] = 0;}}}// ... (剩余部分相同)
}

2. 频域分析实例

幅度谱可视化

c

void visualize_magnitude_spectrum(const char* input_path, const char* output_path) {int width, height;unsigned char *img = stbi_load(input_path, &width, &height, NULL, 1);fftw_complex *in = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * width * height);fftw_complex *out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * width * height);// 准备输入数据for (int y = 0; y < height; y++) {for (int x = 0; x < width; x++) {in[y*width + x][0] = img[y*width + x] / 255.0;in[y*width + x][1] = 0;}}// 正向FFTfftw_plan plan_forward = fftw_plan_dft_2d(height, width, in, out, FFTW_FORWARD, FFTW_ESTIMATE);fftw_execute(plan_forward);// 计算对数幅度谱并归一化double max_log_mag = -INFINITY;double *log_magnitude = (double*)malloc(width * height * sizeof(double));for (int i = 0; i < width*height; i++) {double re = out[i][0];double im = out[i][1];log_magnitude[i] = log(1.0 + sqrt(re*re + im*im));if (log_magnitude[i] > max_log_mag) max_log_mag = log_magnitude[i];}// 转换为图像unsigned char *spectrum_img = (unsigned char*)malloc(width * height);for (int i = 0; i < width*height; i++) {spectrum_img[i] = (unsigned char)(255 * log_magnitude[i] / max_log_mag);}stbi_write_png(output_path, width, height, 1, spectrum_img, width);// 清理资源fftw_destroy_plan(plan_forward);fftw_free(in);fftw_free(out);free(log_magnitude);free(spectrum_img);stbi_image_free(img);
}

//幅度谱可视化
void visualize_magnitude_spectrum(const QString& input_path, const QString& output_path) {// 加载图像QImage inputImg(input_path);if (inputImg.isNull()) {qDebug() << "Error loading image";return;}// 转换为灰度图像QImage grayImg = inputImg.convertToFormat(QImage::Format_Grayscale8);int width = grayImg.width();int height = grayImg.height();// 分配FFTW数组fftw_complex *in = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * width * height);fftw_complex *out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * width * height);// 准备输入数据for (int y = 0; y < height; y++) {const uchar* scanLine = grayImg.scanLine(y);for (int x = 0; x < width; x++) {in[y*width + x][0] = scanLine[x] / 255.0;in[y*width + x][1] = 0;}}// 正向FFTfftw_plan plan_forward = fftw_plan_dft_2d(height, width, in, out, FFTW_FORWARD, FFTW_ESTIMATE);fftw_execute(plan_forward);// 计算对数幅度谱并归一化double max_log_mag = -INFINITY;double *log_magnitude = (double*)malloc(width * height * sizeof(double));for (int i = 0; i < width*height; i++) {double re = out[i][0];double im = out[i][1];log_magnitude[i] = log(1.0 + sqrt(re*re + im*im));if (log_magnitude[i] > max_log_mag) max_log_mag = log_magnitude[i];}// 转换为图像QImage spectrumImg(width, height, QImage::Format_Grayscale8);for (int y = 0; y < height; y++) {uchar* scanLine = spectrumImg.scanLine(y);for (int x = 0; x < width; x++) {int i = y * width + x;scanLine[x] = static_cast<uchar>(255 * log_magnitude[i] / max_log_mag);}}// 保存图像spectrumImg.save(output_path);// 清理资源fftw_destroy_plan(plan_forward);fftw_free(in);fftw_free(out);free(log_magnitude);
}

相位谱可视化

c

void visualize_phase_spectrum(const char* input_path, const char* output_path) {// ... (前面部分与幅度谱相同)// 计算相位谱并归一化到[0,255]unsigned char *phase_img = (unsigned char*)malloc(width * height);for (int i = 0; i < width*height; i++) {double phase = atan2(out[i][1], out[i][0]); // [-π, π]phase_img[i] = (unsigned char)(255 * (phase + M_PI) / (2 * M_PI));}stbi_write_png(output_path, width, height, 1, phase_img, width);// ... (清理资源)
}

void visualize_phase_spectrum(const QString& input_path, const QString& output_path) {// 加载图像QImage inputImg(input_path);if (inputImg.isNull()) {qDebug() << "Error loading image";return;}// 转换为灰度图像QImage grayImg = inputImg.convertToFormat(QImage::Format_Grayscale8);int width = grayImg.width();int height = grayImg.height();// 分配FFTW数组fftw_complex *in = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * width * height);fftw_complex *out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * width * height);// 准备输入数据for (int y = 0; y < height; y++) {const uchar* scanLine = grayImg.scanLine(y);for (int x = 0; x < width; x++) {in[y*width + x][0] = scanLine[x] / 255.0;in[y*width + x][1] = 0;}}// 正向FFTfftw_plan plan_forward = fftw_plan_dft_2d(height, width, in, out, FFTW_FORWARD, FFTW_ESTIMATE);fftw_execute(plan_forward);// 计算相位谱并归一化到[0,255]QImage phaseImg(width, height, QImage::Format_Grayscale8);for (int y = 0; y < height; y++) {uchar* scanLine = phaseImg.scanLine(y);for (int x = 0; x < width; x++) {int i = y * width + x;double phase = atan2(out[i][1], out[i][0]); // [-π, π]scanLine[x] = static_cast<uchar>(255 * (phase + M_PI) / (2 * M_PI));}}// 保存图像phaseImg.save(output_path);// 清理资源fftw_destroy_plan(plan_forward);fftw_free(in);fftw_free(out);
}

功率谱分析

c

void analyze_power_spectrum(const char* input_path) {// ... (加载图像和FFT部分与前面相同)// 计算径向功率谱分布int max_radius = (int)sqrt(width*width + height*height)/2;double *power_dist = (double*)calloc(max_radius, sizeof(double));int *count = (int*)calloc(max_radius, sizeof(int));int center_x = width / 2;int center_y = height / 2;for (int y = 0; y < height; y++) {for (int x = 0; x < width; x++) {int dx = x - center_x;int dy = y - center_y;int r = (int)sqrt(dx*dx + dy*dy);if (r < max_radius) {double power = out[y*width + x][0]*out[y*width + x][0] + out[y*width + x][1]*out[y*width + x][1];power_dist[r] += power;count[r]++;}}}// 输出径向功率分布for (int r = 0; r < max_radius; r++) {if (count[r] > 0) {printf("半径 %d: 平均功率 %f\n", r, power_dist[r]/count[r]);}}// ... (清理资源)free(power_dist);free(count);
}

3. 使用示例

c

int main() {// 低通滤波示例lowpass_filter("input.png", "output_lowpass.png", 30);// 高通滤波示例highpass_filter("input.png", "output_highpass.png", 20);// 带通滤波示例bandpass_filter("input.png", "output_bandpass.png", 10, 50);// 频域分析示例visualize_magnitude_spectrum("input.png", "magnitude_spectrum.png");visualize_phase_spectrum("input.png", "phase_spectrum.png");analyze_power_spectrum("input.png");return 0;
}

4. 实际应用技巧

1. 频域中心化：在显示频谱前，可以使用fftshift将低频移到图像中心

   c

   void fftshift(fftw_complex *data, int width, int height) {for (int y = 0; y < height/2; y++) {for (int x = 0; x < width/2; x++) {// 交换四个象限swap(&data[y*width + x], &data[(y+height/2)*width + (x+width/2)]);swap(&data[y*width + (x+width/2)], &data[(y+height/2)*width + x]);}}
   }

2. 优化技巧：

   • 对于实数输入图像，使用fftw_plan_dft_r2c_2d和fftw_plan_dft_c2r_2d更高效

   • 重用FFTW计划(plan)避免重复创建开销

   • 对于大图像，考虑使用多线程FFTW

3. 窗口函数：应用汉宁窗等可以减少频谱泄漏

   c

   void apply_hann_window(fftw_complex *in, int width, int height) {for (int y = 0; y < height; y++) {double wy = 0.5 * (1 - cos(2*M_PI*y / (height-1)));for (int x = 0; x < width; x++) {double wx = 0.5 * (1 - cos(2*M_PI*x / (width-1)));in[y*width + x][0] *= wx * wy;in[y*width + x][1] *= wx * wy;}}
   }