数字视频处理（五）——频率域陷波滤波
2024-05-15 14:57:00
- 编程实现2-d DFT正变换和反变换
- 频率域陷波滤波
实验代码
解决方案资源管理器如下：
FFT.h
#pragma once
void compute_W(int n, double* W_re, double* W_im);
void permute_bitrev(int n, double* A_re, double* A_im);
int bitrev(int inp, int numbits);
int log_2(int n);
void fft_butterfly(int n, double* A_re, double* A_im, double* W_re, double* W_im);void fft_1d(int N, double* A_re, double* A_im);
void Notch(unsigned char* IYmoon, double* IYNotch,int moonheight, int moonwidth);
void Notch_2(double* IYmoon, double* IYNotch, int moonheight, int moonwidth);
void Notch_2(double* NY, double* IYNotch, int moonheight, int moonwidth);
void Change_Y(double* IYNotch, double* FFT2D_IM, double* FFT2D_RE, int moonheight, int moonwidth, int FFTheight,int FFTwidth);
void FFT_Column(double* FFT2D_RE, double* FFT2D_IM, int FFTheight, int FFTwidth);
void FFT_Row(double* FFT2D_RE, double* FFT2D_IM, int FFTheight, int FFTwidth);
void CreateH(double* H, int FFTheight, int FFTwidth);
void Filter(double* FFT2D_IM, double* FFT2D_RE, double* H, int FFTheight, int FFTwidth);
void Conjugate(double* FFT2D_IM, int FFTsize);
void Delete_zero(double* FFT2D_RE, double* NY, int FFTheight, int FFTwidth, int moonheight, int moonwidth);
fft.cpp
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>
#include <iostream>
#include <cstdio>
#include <fstream>
#include <cmath>
#include "FFT.h"#define M_PI 3.1415926/*频率陷波*/
void Notch(unsigned char* IYmoon, double* IYNotch, int moonheight, int moonwidth)
{for (int i = 0; i < moonheight; i++){for (int j = 0; j < moonwidth; j++){IYNotch[i * moonwidth + j] = pow(-1,i+j) * double(IYmoon[i * moonwidth + j]);}}
}
void Notch_2(double* NY, double* IYNotch, int moonheight, int moonwidth)
{for (int i = 0; i < moonheight; i++){for (int j = 0; j < moonwidth; j++){IYNotch[i * moonwidth + j] = pow(-1, i + j) * NY[i * moonwidth + j];}}
}
void Delete_zero(double* FFT2D_RE, double* NY, int FFTheight, int FFTwidth, int moonheight, int moonwidth)
{for (int i = 0; i < moonheight; i++){for (int j = 0; j < moonwidth; j++){NY[i * moonwidth + j] = FFT2D_RE[i * FFTwidth + j]/((double)FFTheight*FFTwidth);}}
}
void Change_Y(double* IYNotch, double* FFT2D_IM, double* FFT2D_RE, int moonheight, int moonwidth, int FFTheight, int FFTwidth)
{for (int i = 0; i < FFTheight; i++){for (int j = 0; j < FFTwidth; j++){if (i < moonheight && j < moonwidth){FFT2D_RE[i * FFTwidth + j] = IYNotch[i * moonwidth + j];}elseFFT2D_RE[i * FFTwidth + j] = 0;}}
}
void FFT_Column(double* FFT2D_RE, double* FFT2D_IM, int FFTheight, int FFTwidth)
{for (int i = 0; i < FFTwidth; i++)//对每列做FFT_1D变换{double* A_re = new double[sizeof(double) * FFTheight];double* A_im = new double[sizeof(double) * FFTheight];for (int k = 0; k < FFTheight; k++){A_re[k] = FFT2D_RE[k * FFTwidth + i];A_im[k] = FFT2D_IM[k * FFTwidth + i];}fft_1d(FFTheight, A_re, A_im);for (int k = 0; k < FFTheight; k++){FFT2D_RE[k * FFTwidth + i] = A_re[k];FFT2D_IM[k * FFTwidth + i] = A_im[k];}delete[]A_re;delete[]A_im;}
}void FFT_Row(double* FFT2D_RE, double* FFT2D_IM, int FFTheight, int FFTwidth)
{for (int i = 0; i < FFTheight; i++)//对每行做FFT_1D变换{double* A_re = new double[sizeof(double) * FFTwidth];double* A_im = new double[sizeof(double) * FFTwidth];for (int k = 0; k < FFTwidth; k++){A_re[k] = FFT2D_RE[i * FFTwidth + k];A_im[k] = FFT2D_IM[i * FFTwidth + k];}fft_1d(FFTwidth, A_re, A_im);for (int k = 0; k < FFTwidth; k++){FFT2D_RE[i * FFTwidth + k] = A_re[k];FFT2D_IM[i * FFTwidth + k] = A_im[k];}delete[]A_re;delete[]A_im;}
}void Filter(double* FFT2D_IM, double* FFT2D_RE, double* H, int FFTheight, int FFTwidth)
{for (int i = 0; i < FFTheight; i++){for (int j = 0; j < FFTwidth; j++){FFT2D_IM[i * FFTwidth + j] = FFT2D_IM[i * FFTwidth + j] * H[i * FFTwidth + j];}}
}void CreateH(double* H, int FFTheight, int FFTwidth)
{for (int i = 0; i < FFTheight; i++){for (int j = 0; j < FFTwidth; j++){H[i * FFTwidth + j] = 1.5;if ((i == FFTheight / 2) && (j == FFTwidth / 2)){H[i * FFTwidth + j] = 0.5;}}}
}void Conjugate(double* FFT2D_IM, int FFTsize)
{for (int i = 0; i < FFTsize; i++){FFT2D_IM[i] = -FFT2D_IM[i];}
}
/*********************************************************************************************
对N点序列进行fft变换：
1. A_re为该序列的实部，A_im为该序列的虚部；
2. 运算结果仍然存放在A_re和A_im
3. 输入输出都是自然顺序Note：
N必须为2的整数次幂
**********************************************************************************************/
void fft_1d(int N, double* A_re, double* A_im)
{double* W_re, * W_im;W_re = (double*)malloc(sizeof(double) * N / 2);W_im = (double*)malloc(sizeof(double) * N / 2);assert(W_re != NULL && W_im != NULL);compute_W(N, W_re, W_im);fft_butterfly(N, A_re, A_im, W_re, W_im);permute_bitrev(N, A_re, A_im);free(W_re);free(W_im);
}/* W will contain roots of unity so that W[bitrev(i,log2n-1)] = e^(2*pi*i/n)
* n should be a power of 2
* Note: W is bit-reversal permuted because fft(..) goes faster if this is done.
* see that function for more details on why we treat 'i' as a (log2n-1) bit number.
*/
void compute_W(int n, double* W_re, double* W_im)
{int i, br;int log2n = log_2(n);for (i = 0; i < (n / 2); i++){br = bitrev(i, log2n - 1);W_re[br] = cos(((double)i * 2.0 * M_PI) / ((double)n));W_im[br] = -sin(((double)i * 2.0 * M_PI) / ((double)n));}}/* permutes the array using a bit-reversal permutation */
void permute_bitrev(int n, double* A_re, double* A_im)
{int i, bri, log2n;double t_re, t_im;log2n = log_2(n);for (i = 0; i < n; i++){bri = bitrev(i, log2n);/* skip already swapped elements */if (bri <= i) continue;t_re = A_re[i];t_im = A_im[i];A_re[i] = A_re[bri];A_im[i] = A_im[bri];A_re[bri] = t_re;A_im[bri] = t_im;}
}/* treats inp as a numbits number and bitreverses it.
* inp < 2^(numbits) for meaningful bit-reversal
*/
int bitrev(int inp, int numbits)
{int i, rev = 0;for (i = 0; i < numbits; i++){rev = (rev << 1) | (inp & 1);inp >>= 1;}return rev;
}/* returns log n (to the base 2), if n is positive and power of 2 */
int log_2(int n)
{int res;for (res = 0; n >= 2; res++)n = n >> 1;return res;
}/* fft on a set of n points given by A_re and A_im. Bit-reversal permuted roots-of-unity lookup table
* is given by W_re and W_im. More specifically, W is the array of first n/2 nth roots of unity stored
* in a permuted bitreversal order.
*
* FFT - Decimation In Time FFT with input array in correct order and output array in bit-reversed order.
*
* REQ: n should be a power of 2 to work.
*
* Note: - See www.cs.berkeley.edu/~randit for her thesis on VIRAM FFTs and other details about VHALF section of the algo
* (thesis link - http://www.cs.berkeley.edu/~randit/papers/csd-00-1106.pdf)
* - See the foll. CS267 website for details of the Decimation In Time FFT implemented here.
* (www.cs.berkeley.edu/~demmel/cs267/lecture24/lecture24.html)
* - Also, look "Cormen Leicester Rivest [CLR] - Introduction to Algorithms" book for another variant of Iterative-FFT
*/
void fft_butterfly(int n, double* A_re, double* A_im, double* W_re, double* W_im)
{double w_re, w_im, u_re, u_im, t_re, t_im;int m, g, b;int mt, k;/* for each stage */for (m = n; m >= 2; m = m >> 1){/* m = n/2^s; mt = m/2; */mt = m >> 1;/* for each group of butterfly */for (g = 0, k = 0; g < n; g += m, k++){/* each butterfly group uses only one root of unity. actually, it is the bitrev of this group's number k.* BUT 'bitrev' it as a log2n-1 bit number because we are using a lookup array of nth root of unity and* using cancellation lemma to scale nth root to n/2, n/4,... th root.** It turns out like the foll.* w.re = W[bitrev(k, log2n-1)].re;* w.im = W[bitrev(k, log2n-1)].im;* Still, we just use k, because the lookup array itself is bit-reversal permuted.*/w_re = W_re[k];w_im = W_im[k];/* for each butterfly */for (b = g; b < (g + mt); b++){/* t = w * A[b+mt] */t_re = w_re * A_re[b + mt] - w_im * A_im[b + mt];t_im = w_re * A_im[b + mt] + w_im * A_re[b + mt];/* u = A[b]; in[b] = u + t; in[b+mt] = u - t; */u_re = A_re[b];u_im = A_im[b];A_re[b] = u_re + t_re;A_im[b] = u_im + t_im;A_re[b + mt] = u_re - t_re;A_im[b + mt] = u_im - t_im;}}}
}main.cpp
#include <iostream>
#include <cstdio>
#include <fstream>
#include "FFT.h"using namespace std;
int main()
{int moonwidth = 464;int moonheight = 538;int FFTwidth = 512;int FFTheight = 1024;int FFTsize = FFTheight * FFTwidth;int moonsize = moonheight * moonwidth * 3;int moonYsize = moonheight * moonwidth;int moonEsize = moonheight * moonwidth * 2;ifstream Imoonfile("moon.yuv", ios::binary);ofstream Omoonfile("outmoon.yuv", ios::binary);if (!Imoonfile) { cout << "error to open moonfile!" << endl; }if (!Omoonfile) { cout << "error to create moonfile!" << endl; }unsigned char* IYmoon = new unsigned char[moonYsize];//moon的Y分量unsigned char* IEmoon = new unsigned char[moonEsize];//moon的其他分量unsigned char* OYmoon = new unsigned char[moonYsize];double* IYNotch = new double[moonYsize];//中心化后的Y分量double* FFT2D_RE = new double[FFTsize];//FFT以后的矩阵实部double* FFT2D_IM = new double[FFTsize];//FFT以后的矩阵虚部double* NY = new double[moonYsize];//反变换以后的Y分量Imoonfile.read((char*)IYmoon, moonYsize);Imoonfile.read((char*)IEmoon, moonEsize);Notch(IYmoon,IYNotch, moonheight,moonwidth);//补零Change_Y(IYNotch, FFT2D_IM, FFT2D_RE,moonheight,moonwidth,FFTheight,FFTwidth);//进行每列的FFT_1D变换FFT_Column(FFT2D_RE,FFT2D_IM, FFTheight, FFTwidth);FFT_Row(FFT2D_RE, FFT2D_IM, FFTheight, FFTwidth);//用滤波器函数H(u, v)乘以F(u, v)double* H = new double[FFTsize];CreateH(H,FFTheight,FFTwidth);Filter(FFT2D_IM, FFT2D_RE, H, FFTheight, FFTwidth);//虚部取共轭Conjugate(FFT2D_IM, FFTsize);FFT_Column(FFT2D_RE, FFT2D_IM, FFTheight, FFTwidth);FFT_Row(FFT2D_RE, FFT2D_IM, FFTheight, FFTwidth);Conjugate(FFT2D_IM, FFTsize);//去零Delete_zero(FFT2D_RE, NY, FFTheight, FFTwidth, moonheight, moonwidth);//去中心化Notch_2(NY, IYNotch, moonheight, moonwidth);//强制类型转换for (int i = 0; i < moonheight; i++){for (int j = 0; j < moonwidth; j++){if (IYNotch[i * moonwidth + j] > 255){IYNotch[i * moonwidth + j] = 255;}else if (IYNotch[i * moonwidth + j] < 0){IYNotch[i * moonwidth + j] = 0;}OYmoon[i * moonwidth + j] = unsigned char(IYNotch[i * moonwidth + j]);}}//写文件Omoonfile.write((char*)OYmoon, moonYsize);Omoonfile.write((char*)IEmoon, moonEsize);delete[]IYmoon;delete[]IEmoon;delete[]IYNotch;delete[]FFT2D_IM;delete[]FFT2D_RE;delete[]NY;delete[]H;Imoonfile.close();Omoonfile.close();return 0;
}实验结果
至此，实验结束。
数字视频处理（五）——频率域陷波滤波相关推荐
- 【OpenCV 例程200篇】90. 频率域陷波滤波器
[OpenCV 例程200篇]90. 频率域陷波滤波器 欢迎关注 『OpenCV 例程200篇』 系列,持续更新中 欢迎关注 『Python小白的OpenCV学习课』 系列,持续更新中 5.2 陷波滤 ...
- 《数字图像处理》手动实现最佳陷波滤波
1 最佳陷波滤波实现 1.1 最佳陷波滤波原理及步骤 最佳陷波滤波的可以良好地处理一个以上的干扰分量或者多个周期性的噪声,相比于其他的滤波方法,最佳陷波滤波可以最小化复原的估计值f ̂(x,y)的局部 ...
- 【OpenCV 例程200篇】83. 频率域低通滤波：印刷文本字符修复
[OpenCV 例程200篇]83. 频率域低通滤波案例:印刷文本字符修复 欢迎关注 『OpenCV 例程200篇』 系列,持续更新中 欢迎关注 『Python小白的OpenCV学习课』 系列,持续更 ...
- 第5章 Python 数字图像处理(DIP) - 图像复原与重建12 - 空间滤波 - 使用频率域滤波降低周期噪声 - 陷波滤波、最优陷波滤波
标题 使用频率域滤波降低周期噪声 陷波滤波深入介绍 最优陷波滤波 本章陷波滤波器有部分得出的结果不佳,如果有更好的解决方案,请赐教,不胜感激. 使用频率域滤波降低周期噪声 陷波滤波深入介绍 零相移滤波 ...
- MATLAB高斯陷波滤波图像
clc,clear,close all % 清理命令区.清理工作区.关闭显示图形 warning off % 消除警告 feature jit off % 加速代码运行 D0 = 4; % 阻止的频率 ...
- MATLAB巴特沃斯陷波滤波图像
clc,clear,close all % 清理命令区.清理工作区.关闭显示图形 warning off % 消除警告 feature jit off % 加速代码运行 D0 = 4; % 阻止的频率 ...
- MATLAB设计指定频率的陷波器
下面来介绍一下如何在MATLAB中编程实现陷波器,具体如下: 1.首先在MATLAB的主界面编辑器中写入下列代码: clf;clear; %设置初值 f0=70; %频率为70HZ Ts=0.001; ...
- 使用陷波滤波减少莫尔 (波纹) 模式（C++）
处理图片中的莫尔波纹,需要用到陷波滤波器,本文使用巴特沃斯滤波器设计滤波器,并且可以通过opencv库调用鼠标callback操作,捕捉频域中的噪音点.在频域对图片进行处理,达到很好的滤波效果. #i ...
- 在频率域进行音频滤波（Python）
最近一个项目需要对声音文件进行一个滤波处理,由于并没有这种经验,所以借助图像的滤波来实现一下 首先导入需要的包和设置路径 import wave as we import numpy as np im ...
最新文章
- linux系统下redhat7之虚拟机控制
- 神经网络反向传播算法
- leader选举的源码分析-QuorumPeer.start
- linux安装python3.6 setuptools_linux下安装Python3.6.1
- Go程序：利用命令行参数做四则运算
- Python中利用*打印不同的三角形
- jar httpclient 少包，此处 区别 common-http包
- PAT 1018 锤子剪刀布
- shop--8.商品类别--批量操作--添加（前端）
- matlab遥感代码,遥感融合定量评价matlab程序代码
- 计算机软件在哪里建文本文档,电脑点击右键的新建文本文档不见了的解决方法 怎么解决电脑点击右键的新建文本文档不见了...
- Funcode学习笔记：完成Run、Jump、Idle等动作【后续更新Roll、Attack动作的实现】【By Myself】
- 【3D模型分享】柴油机MMZ D-260柴油发动机
- postgresql 解锁表
- 电商项目必备技能=＞放大镜
- 正则匹配里面的(.*?)
- 【ThreeJs】（2）照相机 | 正交投影照相机 | 透视投影照相机
- 基于CNN 对车牌数字进行识别，（二）
- ipad怎么压缩文件？教你一招快捷压缩图片
- Android学习之ImageView放置gif动态图