先说下问题:
我用8x8的数据测试FFT2和IFFT2两个算法,发现数据未还原,求帮忙检查下代码哪写错了?
(注:如果只对每一行或者每一列做FFT和IFFT,可以还原数据。)
打印的日志下如
封装的傅里叶变换类
using System;
/// <summary>
/// 傅里叶变换
/// </summary>
public sealed class Fourier
{
//快速傅里叶变换
public static void FFT(Complex[] TD2FD)
{
FFT_Core(TD2FD, WT_LUT(TD2FD.Length, 1));
}
//快速傅里叶变换 (二维)
public static void FFT2(Complex2D TD2FD)
{
//对每一行做FFT
for (int i = 0; i < TD2FD.Height; i++)
{
Complex[] row = TD2FD.GetRow(i);
FFT(row);
TD2FD.SetRow(i, row);
}
//对每一列做FFT
for (int i = 0; i < TD2FD.Width; i++)
{
Complex[] column = TD2FD.GetColumn(i);
FFT(column);
TD2FD.SetColumn(i, column);
}
}
//快速傅里叶逆变换
public static void IFFT(Complex[] FD2TD)
{
//做FFT变换
Complex[] WT = WT_LUT(FD2TD.Length, -1);
FFT_Core(FD2TD, WT);
//实部除以N
for (int i = 0; i < FD2TD.Length; i++)
FD2TD[i].re /= FD2TD.Length;
}
//快速傅里叶逆变换 (二维)
public static void IFFT2(Complex2D FD2TD)
{
//对每一行做IFFT
for (int i = 0; i < FD2TD.Height; i++)
{
Complex[] row = FD2TD.GetRow(i);
IFFT(row);
FD2TD.SetRow(i, row);
}
//对每一列做IFFT
for (int i = 0; i < FD2TD.Width; i++)
{
Complex[] column = FD2TD.GetColumn(i);
IFFT(column);
FD2TD.SetColumn(i, column);
}
}
// 将直流分量移到频谱图的中心
public static void FFT2Shift(Complex2D complex2D)
{
int halfH = complex2D.Height / 2;
int halfW = complex2D.Width / 2;
//将图像四个象限区域按对角线交换
for (int i=0; i<halfH; i++)
{
for (int j=0; j<complex2D.Width; j++)
{
if (j < halfW)
complex2D.SwapComplex(i, j, i + halfH, j + halfW);
else
complex2D.SwapComplex(i, j, i + halfH, j - halfW);
}
}
}
// 高通滤波
public static void HighPassFilting(Complex2D complex2D)
{
int halfH = complex2D.Height / 2;
int halfW = complex2D.Width / 2;
int H4 = complex2D.Height / 8;
int W4 = complex2D.Width / 8;
for (int i = halfH - H4; i < halfH + H4; i++)
{
for (int j = halfW - W4; j < halfW + W4; j++)
{
Complex cpx = complex2D.GetComplex(i, j);
cpx.re = 0;
cpx.im = 0;
complex2D.SetComplex(i, j, cpx);
}
}
}
// 低通滤波
public static void LowPassFilting(Complex2D complex2D)
{
int halfH = complex2D.Height / 2;
int halfW = complex2D.Width / 2;
int H4 = complex2D.Height / 8;
int W4 = complex2D.Width / 8;
for (int i=0; i < complex2D.Height; i++)
{
for (int j=0; j < complex2D.Width; j++)
{
if (i < halfH - H4 || i > halfH + H4 ||
j < halfW - W4 || j > halfW + W4)
{
Complex cpx = complex2D.GetComplex(i, j);
cpx.re = 0;
cpx.im = 0;
complex2D.SetComplex(i, j, cpx);
}
}
}
}
// 返回旋转因子查询表(twiddle factor lookup table)
private static Complex[] WT_LUT(int N, int flag = 1)
{
Complex[] WT = new Complex[N];
for (int i = 0; i < N; i++)
{
Complex cpx_wt = new Complex();
float angle = (float)(-i * Math.PI * 2 / N);
cpx_wt.re = (float)Math.Cos(angle);
//IFFT flag=-1, FFT flag=1
cpx_wt.im = flag * (float)Math.Sin(angle);
WT[i] = cpx_wt;
}
return WT;
}
private static void FFT_Core(Complex[] TD2FD, Complex[] WT)
{
int power = (int)Math.Log(TD2FD.Length, 2);
int butterfly;
int p, s;
Complex x1, x2, wt;
//倒位排序
BitReverse(TD2FD);
//蝶形运算
for (int k = 0; k < power; k++) //级数
{
for (int j = 0; j < 1 << (power - k - 1); j++) //组数
{
//每组有几个元素
butterfly = 1 << k + 1;
//计算参与蝶形运算的两个元素的索引
p = j * butterfly;
s = p + butterfly / 2;
for (int i = 0; i < butterfly / 2; i++) //蝶形运算次数
{
x1 = TD2FD[i + p];
x2 = TD2FD[i + s];
wt = WT[i * TD2FD.Length / butterfly];
TD2FD[i + p] = x1 + x2 * wt;
TD2FD[i + s] = x1 - x2 * wt;
}
}
}
}
// 倒位排序——雷德算法
private static void BitReverse(Complex[] array)
{
//倒位排序原理
//0 1 2 3 4 5 6 7 十进制
//000 001 010 011 100 101 110 111 十进制对应的二进制
//000 100 010 110 001 101 011 111 码位反转
//0 4 2 6 1 5 3 7 码位反转后对应的十进制
int i, j, k;
int N = array.Length;
Complex temp;
j = 0;
for (i = 0; i < N - 1; i++)
{
if (i < j)
{
temp = array[i];
array[i] = array[j];
array[j] = temp;
}
// 求j的下一个倒序位
// N/2的二进制最高位为1,其他位都为0
// 判断最高位是否为1,可与N/2进行比较
// 判断次高位是否为1,可与N/4进行比较
k = N >> 1;
//如果k<=j,表示j的最高位为1
while (k <= j)
{
//当k<=j时,需要将最高位变为0
j = j - k;
//判断次高位是否为1,依次类推,逐个比较,直到j某个位为0
k >>= 1;
}
j = j + k;//将0变为1
}
}
// 打印
public static void Print(Complex[] TD2FD)
{
for (int i = 0; i < TD2FD.Length; i++)
{
Console.WriteLine(TD2FD[i].ToString());
}
Console.WriteLine();
}
}
//定义复数
public class Complex
{
public float re;//实数部
public float im;//虚数部
// 幅值
public double Amplitude
{
get
{
//测试发现取值范围为
//min=0.0009918213, max=412.4615
return Math.Sqrt(re * re + im * im);
}
}
// 频谱图像素值
public double PixelAmplitude
{
get
{
//幅值范围很大,需要做以下处理:
//1. 将幅值范围调到 [1, ?]
//2. 利用Log函数压缩范围
//3. 将范围映射到颜色值[0,1]
double p = Math.Log(Amplitude * 10000) / 16f;
return p;
}
}
// 相位
public double Phase
{
get
{
return Math.Atan2(im, re);
}
}
public override string ToString()
{
return string.Format("re={0}, im={1}", re, im);
}
public static Complex operator +(Complex lhs, Complex rhs)
{
Complex result = new Complex();
result.re = lhs.re + rhs.re;
result.im = lhs.im + rhs.im;
return result;
}
public static Complex operator -(Complex lhs, Complex rhs)
{
Complex result = new Complex();
result.re = lhs.re - rhs.re;
result.im = lhs.im - rhs.im;
return result;
}
public static Complex operator *(Complex lhs, Complex rhs)
{
Complex result = new Complex();
result.re = lhs.re * rhs.re - lhs.im * rhs.im;
result.im = lhs.re * rhs.im + lhs.im * rhs.re;
return result;
}
public static Complex operator *(float lhs, Complex rhs)
{
Complex result = new Complex();
result.re = lhs * rhs.re;
result.im = lhs * rhs.im;
return result;
}
public static Complex operator *(Complex lhs, float rhs)
{
Complex result = new Complex();
result.re = lhs.re * rhs;
result.im = lhs.im * rhs;
return result;
}
}
public class Complex2D
{
private Complex[,] matrix;
private int m_width;
private int m_height;
// width:图像宽度 height:图像高度
public Complex2D(int width, int height)
{
m_width = width;
m_height = height;
matrix = new Complex[width, height];
}
public int Width { get { return m_width; } }
public int Height { get { return m_height; } }
public Complex[] GetRow(int i)
{
Complex[] row = new Complex[m_width];
for (int j = 0; j < m_width; j++)
row[j] = matrix[j,i];
return row;
}
public void SetRow(int i, Complex[] array)
{
for (int j = 0; j < m_width; j++)
matrix[j, i] = array[j];
}
public Complex[] GetColumn(int i)
{
Complex[] column = new Complex[m_height];
for (int j = 0; j < m_height; j++)
column[j] = matrix[i,j];
return column;
}
public void SetColumn(int i, Complex[] array)
{
for (int j = 0; j < m_width; j++)
matrix[i, j] = array[j];
}
//i: 第几行 j: 第几列
public Complex GetComplex(int i, int j)
{
return matrix[j,i];
}
//i: 第几行 j: 第几列
public void SetComplex(int i, int j, Complex src)
{
matrix[j, i] = src;
}
// 交换两个元素
// i: 第几行 j: 第几列
public void SwapComplex(int i0, int j0, int i1, int j1)
{
Complex tmp = matrix[j0,i0];
matrix[j0, i0] = matrix[j1, i1];
matrix[j1, i1] = tmp;
}
// 交换行
public void SwapRow(int i, int j)
{
for (int k=0; k<m_width; k++)
{
Complex cpx0 = matrix[k,i];
Complex cpx1 = matrix[k,j];
matrix[k,i] = cpx1;
matrix[k,j] = cpx0;
}
}
// 交换列
public void SwapColumn(int i, int j)
{
for (int k = 0; k < m_height; k++)
{
Complex cpx0 = matrix[i,k];
Complex cpx1 = matrix[j,k];
matrix[i,k] = cpx1;
matrix[j,k] = cpx0;
}
}
public void Print(string fileName)
{
System.Text.StringBuilder sb = new System.Text.StringBuilder();
for (int i = 0; i < m_height; i++)
{
for (int j = 0; j < m_width; j++)
sb.AppendFormat("{0:G} ", matrix[j,i].re.ToString().PadRight(5));
sb.AppendLine();
}
System.IO.File.WriteAllText(string.Format("D://{0}.txt", fileName), sb.ToString());
}
}
测试代码
//测试
Complex2D complex2D = new Complex2D(8, 8);
for (int i = 0; i < complex2D.Height; i++)
{
for (int j = 0; j < complex2D.Width; j++)
{
Complex cpx = new Complex();
cpx.re = i * complex2D.Width + j;
cpx.im = 0;
complex2D.SetComplex(i, j, cpx);
}
}
complex2D.Print("8x8");
Fourier.FFT2(complex2D);
complex2D.Print("fft2");
Fourier.IFFT2(complex2D);
complex2D.Print("ifft2");
ok,搞定了;
引发错误的地方在于:傅里叶逆变换,实数部分和虚数部分都需要除N,你只给实数部分除了N,一维数据没有问题是因为,你的一维数据全是实数,逆变换后虚数部分为0,所以不影响;而二维傅里叶变换需要两次逆变换(行和列的逆变换),因此就有了较大影响;
public static void IFFT(Complex[] FD2TD)
{
//做FFT变换
Complex[] WT = WT_LUT(FD2TD.Length, -1);
FFT_Core(FD2TD, WT);
//实部除以N
for (int i = 0; i < FD2TD.Length; i++)
{
FD2TD[i].re /= FD2TD.Length;
FD2TD[i].im /= FD2TD.Length;
}
}
这样改一下就好了,我不太会输出txt,我用控制台输出试了一下,这是逆变换以后的数据:
ok,误差在可忽略的范围内;
还有一个规范性的问题,一般矩阵(row,columns)的行是row(Height),列是columns(Width),你刚开始初始化时确实是这么定义的,但后面你就全部反着来了,因为这次输入的是8*8的数据,行列数相等,没有报错,如果下个数据行列数不相等,就有错了,还是规范起来写的好。(这个问题不影响你的这次的结果,但还是改掉的好)
谢谢,已经将矩阵规范化了,是我修改时改漏了一句。
public class Complex2D
{
private Complex[,] matrix;
private int m_width;
private int m_height;
// width:图像宽度 height:图像高度
public Complex2D(int width, int height)
{
m_width = width;
m_height = height;
matrix = new Complex[height,width];
}
public int Width { get { return m_width; } }
public int Height { get { return m_height; } }
public Complex[] GetRow(int i)
{
Complex[] row = new Complex[m_width];
for (int j = 0; j < m_width; j++)
row[j] = matrix[i,j];
return row;
}
public void SetRow(int i, Complex[] array)
{
for (int j = 0; j < m_width; j++)
matrix[i,j] = array[j];
}
public Complex[] GetColumn(int i)
{
Complex[] column = new Complex[m_height];
for (int j = 0; j < m_height; j++)
column[j] = matrix[j,i];
return column;
}
public void SetColumn(int i, Complex[] array)
{
for (int j = 0; j < m_height; j++)
matrix[j,i] = array[j];
}
//i: 第几行 j: 第几列
public Complex GetComplex(int i, int j)
{
return matrix[i,j];
}
//i: 第几行 j: 第几列
public void SetComplex(int i, int j, Complex src)
{
matrix[i,j] = src;
}
// 交换两个元素
// i: 第几行 j: 第几列
public void SwapComplex(int i0, int j0, int i1, int j1)
{
Complex tmp = matrix[i0,j0];
matrix[i0, j0] = matrix[i1, j1];
matrix[i1, j1] = tmp;
}
// 交换行
public void SwapRow(int i, int j)
{
for (int k=0; k<m_width; k++)
{
Complex cpx0 = matrix[i,k];
Complex cpx1 = matrix[j,k];
matrix[i,k] = cpx1;
matrix[j,k] = cpx0;
}
}
// 交换列
public void SwapColumn(int i, int j)
{
for (int k = 0; k < m_height; k++)
{
Complex cpx0 = matrix[k,i];
Complex cpx1 = matrix[k,j];
matrix[k,i] = cpx1;
matrix[k,j] = cpx0;
}
}
public void Print(string fileName)
{
System.Text.StringBuilder sb = new System.Text.StringBuilder();
for (int i = 0; i < m_height; i++)
{
for (int j = 0; j < m_width; j++)
sb.AppendFormat("{0:G} ", Math.Floor(matrix[i,j].re).ToString().PadRight(5));
sb.AppendLine();
}
System.IO.File.WriteAllText(string.Format("D://{0}.txt", fileName), sb.ToString());
}
}
额。。。咱俩对行和列的理解有偏差。。。。。
我本来想按照你的行和列来改程序。。我真的做不到,思维定势,弄到最后我自己混乱了。。
你这个问题就是行列有点乱导致的(当然了,也我不排除是我的锅,可能你本来是对的,我的理解和你是反的,把你带偏了)。。
我按照自己的行列从头到尾改了一遍,我按这个跑出来没问题,附上代码。。你自己好好检查下你的吧,
using System;
/// <summary>
/// 傅里叶变换
/// </summary>
public sealed class Fourier
{
//快速傅里叶变换
public static void FFT(Complex[] TD2FD)
{
FFT_Core(TD2FD, WT_LUT(TD2FD.Length, 1));
}
//快速傅里叶变换 (二维)
public static void FFT2(Complex2D TD2FD)
{
//对每一行做FFT
for (int i = 0; i < TD2FD.Height; i++)
{
Complex[] row = TD2FD.GetRow(i);
FFT(row);
TD2FD.SetRow(i, row);
}
//对每一列做FFT
for (int i = 0; i < TD2FD.Width; i++)
{
Complex[] column = TD2FD.GetColumn(i);
FFT(column);
TD2FD.SetColumn(i, column);
}
}
//快速傅里叶逆变换
public static void IFFT(Complex[] FD2TD)
{
//做FFT变换
Complex[] WT = WT_LUT(FD2TD.Length, -1);
FFT_Core(FD2TD, WT);
//实部除以N
for (int i = 0; i < FD2TD.Length; i++)
{
FD2TD[i].re /= FD2TD.Length;
FD2TD[i].im /= FD2TD.Length;
}
}
//快速傅里叶逆变换 (二维)
public static void IFFT2(Complex2D FD2TD)
{
//对每一行做IFFT
for (int i = 0; i < FD2TD.Height; i++)
{
Complex[] row = FD2TD.GetRow(i);
IFFT(row);
FD2TD.SetRow(i, row);
}
//对每一列做IFFT
for (int i = 0; i < FD2TD.Width; i++)
{
Complex[] column = FD2TD.GetColumn(i);
IFFT(column);
FD2TD.SetColumn(i, column);
}
}
// 返回旋转因子查询表(twiddle factor lookup table)
private static Complex[] WT_LUT(int N, int flag = 1)
{
Complex[] WT = new Complex[N];
for (int i = 0; i < N; i++)
{
Complex cpx_wt = new Complex();
float angle = (float)(-i * Math.PI * 2 / N);
cpx_wt.re = (float)Math.Cos(angle);
//IFFT flag=-1, FFT flag=1
cpx_wt.im = flag * (float)Math.Sin(angle);
WT[i] = cpx_wt;
}
return WT;
}
private static void FFT_Core(Complex[] TD2FD, Complex[] WT)
{
int power = (int)Math.Log(TD2FD.Length, 2);
int butterfly;
int p, s;
Complex x1, x2, wt;
//倒位排序
BitReverse(TD2FD);
//蝶形运算
for (int k = 0; k < power; k++) //级数
{
for (int j = 0; j < 1 << (power - k - 1); j++) //组数
{
//每组有几个元素
butterfly = 1 << k + 1;
//计算参与蝶形运算的两个元素的索引
p = j * butterfly;
s = p + butterfly / 2;
for (int i = 0; i < butterfly / 2; i++) //蝶形运算次数
{
x1 = TD2FD[i + p];
x2 = TD2FD[i + s];
wt = WT[i * TD2FD.Length / butterfly];
TD2FD[i + p] = x1 + x2 * wt;
TD2FD[i + s] = x1 - x2 * wt;
}
}
}
}
// 倒位排序——雷德算法
private static void BitReverse(Complex[] array)
{
//倒位排序原理
//0 1 2 3 4 5 6 7 十进制
//000 001 010 011 100 101 110 111 十进制对应的二进制
//000 100 010 110 001 101 011 111 码位反转
//0 4 2 6 1 5 3 7 码位反转后对应的十进制
int i, j, k;
int N = array.Length;
Complex temp;
j = 0;
for (i = 0; i < N - 1; i++)
{
if (i < j)
{
temp = array[i];
array[i] = array[j];
array[j] = temp;
}
// 求j的下一个倒序位
// N/2的二进制最高位为1,其他位都为0
// 判断最高位是否为1,可与N/2进行比较
// 判断次高位是否为1,可与N/4进行比较
k = N >> 1;
//如果k<=j,表示j的最高位为1
while (k <= j)
{
//当k<=j时,需要将最高位变为0
j = j - k;
//判断次高位是否为1,依次类推,逐个比较,直到j某个位为0
k >>= 1;
}
j = j + k;//将0变为1
}
}
// 打印
public static void Print(Complex[] TD2FD)
{
for (int i = 0; i < TD2FD.Length; i++)
{
Console.WriteLine(TD2FD[i].ToString());
}
Console.WriteLine();
}
}
//定义复数
public class Complex
{
public float re;//实数部
public float im;//虚数部
// 幅值
public double Amplitude
{
get
{
//测试发现取值范围为
//min=0.0009918213, max=412.4615
return Math.Sqrt(re * re + im * im);
}
}
// 频谱图像素值
public double PixelAmplitude
{
get
{
//幅值范围很大,需要做以下处理:
//1. 将幅值范围调到 [1, ?]
//2. 利用Log函数压缩范围
//3. 将范围映射到颜色值[0,1]
double p = Math.Log(Amplitude * 10000) / 16f;
return p;
}
}
// 相位
public double Phase
{
get
{
return Math.Atan2(im, re);
}
}
public override string ToString()
{
return string.Format("re={0}, im={1}", re, im);
}
public static Complex operator +(Complex lhs, Complex rhs)
{
Complex result = new Complex();
result.re = lhs.re + rhs.re;
result.im = lhs.im + rhs.im;
return result;
}
public static Complex operator -(Complex lhs, Complex rhs)
{
Complex result = new Complex();
result.re = lhs.re - rhs.re;
result.im = lhs.im - rhs.im;
return result;
}
public static Complex operator *(Complex lhs, Complex rhs)
{
Complex result = new Complex();
result.re = lhs.re * rhs.re - lhs.im * rhs.im;
result.im = lhs.re * rhs.im + lhs.im * rhs.re;
return result;
}
public static Complex operator *(float lhs, Complex rhs)
{
Complex result = new Complex();
result.re = lhs * rhs.re;
result.im = lhs * rhs.im;
return result;
}
public static Complex operator *(Complex lhs, float rhs)
{
Complex result = new Complex();
result.re = lhs.re * rhs;
result.im = lhs.im * rhs;
return result;
}
}
public class Complex2D
{
private Complex[,] matrix;
private int m_width;
private int m_height;
// width:图像宽度 height:图像高度
public Complex2D(int height,int width)
{
m_width = width;
m_height = height;
matrix = new Complex[height, width];
}
public int Width { get { return m_width; } }
public int Height { get { return m_height; } }
public Complex[] GetRow(int i)
{
Complex[] row = new Complex[m_width];
for (int j = 0; j < m_width; j++)
row[j] = matrix[i, j];
return row;
}
public void SetRow(int i, Complex[] array)
{
for (int j = 0; j < m_width; j++)
matrix[i, j] = array[j];
}
public Complex[] GetColumn(int i)
{
Complex[] column = new Complex[m_height];
for (int j = 0; j < m_height; j++)
column[j] = matrix[j, i];
return column;
}
public void SetColumn(int i, Complex[] array)
{
for (int j = 0; j < m_height; j++)
matrix[j,i] = array[j];
}
//i: 第几行 j: 第几列
public Complex GetComplex(int i, int j)
{
return matrix[i, j];
}
public void SetComplex(int i, int j, Complex src)
{
matrix[i, j] = src;
}
//i: 第几行 j: 第几列
public void Print()
{
for (int i = 0; i < m_height; i++)
{
for (int j = 0; j < m_width; j++)
{
Console.WriteLine("(re={0}, im={1})", matrix[i,j].re, matrix[i, j].im);
Console.WriteLine();
}
}
}
}
namespace HelloWorld
{
class Program //Main要包含在类当中
{
static void Main(string[] args)
{//测试
Complex2D complex2D = new Complex2D(16, 8);
for (int i = 0; i < complex2D.Height; i++)
{
for (int j = 0; j < complex2D.Width; j++)
{
Complex cpx = new Complex();
cpx.re = i * complex2D.Width + j;
cpx.im = 0;
complex2D.SetComplex(i, j, cpx);
}
}
complex2D.Print();
Fourier.FFT2(complex2D);
complex2D.Print();
Fourier.IFFT2(complex2D);
complex2D.Print();
}
}
}
我的二维数据是这样的 [row,columns] , 第一位维是行,第二维是列,一共有Height行,width列,定义二维数据是【Height,Width】
正确的代码已整理到我的博客上 http://www.devacg.com/?post=1358