数字视频处理（五）——频率域陷波滤波
【数字视频处理（五）——频率域陷波滤波】
- 编程实现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);
intbitrev(int inp, int numbits);
intlog_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
#include
#include
#include
#include
#include
#include
#include
#include "FFT.h"#defineM_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];
}
else
FFT2D_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
#include
#include
#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;
} 实验结果
文章图片
??至此，实验结束。
推荐阅读
- Java|Java OpenCV图像处理之SIFT角点检测详解
- 事件处理程序
- 4月23日海军节，我在青岛等你，一起看强大的中国海军。（如图如视频）
- 视频转换器哪种好用（）
- 最喜6.8.9
- 不懂法，害人终害己
- 腾讯视频（我有一段rap想给你说）
- 爬虫数据处理HTML转义字符
- 百度云极速下载，来体验飞的感觉，还可以看最新动漫的百度云视频哦
- Android|Android BLE蓝牙连接异常处理