图像处理 提取波形 在新窗口重绘波形 FFT变换 绘出变换结果

[vctuxiangchuli]

小辈正在进行图像处理，要求把示波器出图中的波形提取出来，并将波形进行FFT变换，我现在知道部分算法，但是还有很多问题，如怎样在一个新窗口中绘出波形、FFT变换后的频域图，求可供参考的相似的例子。若有源代码，感激不尽。（[e-mail]）

[liumenghappy]

lz，FFT变换的代码已发到你邮箱

[尹成]

#include "dip.h"

void TDipWindow::add(COMPLEX *a,COMPLEX *b,COMPLEX *ret)
{
ret->re = a->re + b->re;
ret->im = a->im + b->im;
}

void TDipWindow::sub(COMPLEX *a,COMPLEX *b,COMPLEX *ret)
{
ret->re = a->re - b->re;
ret->im = a->im - b->im;
}

void TDipWindow::mul(COMPLEX *m,COMPLEX *n,COMPLEX *ret)
{
float a,b,c;

a = ( m->re - m->im ) * n->im;
b = m->re * ( n->re - n->im );
c = m->im * ( n->re + n->im );
ret->re = a + b;
ret->im = a + c;
}

int TDipWindow::reverse(int t,int k)
{
int i,j,m,s;

m = t;
s = 0;
for ( i = 0 ; i < k ; i++ )
{
j = m / 2;
s = s * 2 + ( m - j * 2 );
m = j;
}

return s;
}

void TDipWindow::root(int n)
{
int j;

omega->re = 1;
omega->im = 0;
omega[1].re = (float)cos(2 * PI / n);
omega[1].im = -(float)sin(2 * PI / n);

for ( j = 2 ; j < n/2 ; j++ )
mul(&omega[1],&omega[j-1],&omega[j]);

for ( j = n/2 ; j < n ; j++ )
{
omega[j].re = -omega[j-n/2].re;
omega[j].im = -omega[j-n/2].im;
}
}

void TDipWindow::fft(int n)
{
int s,k,m,l,nv,t,j;
COMPLEX podd,ret;

k = (int)(log(n) / log(2) + 0.5);
nv = n;
m = 1;
for ( l = k-1 ; l >= 0 ; l-- )
{
for ( t = 0 ; t < m * nv ; t+=nv )
for ( j = 0 ; j < nv/2 ; j++ )
{
s = (t+j) / (int)pow(2,l);
s = reverse(s,k);

ret = omega[s];
mul(&ret,&p[t+j+nv/2],&podd);
sub(&p[t+j],&podd,&p[t+j+nv/2]);
add(&p[t+j],&podd,&p[t+j]);
}
m *= 2;
nv /= 2;
}

for ( t = 0 ; t < n ; t++ )
{
s = reverse(t,k);
f[t] = p[s];
}
}

void TDipWindow::nfft(int n)
{
int i;

for ( i = 0 ; i < n ; i++ )
p[i] = f[i];

fft(n);

p[0] = f[0];
for ( i = 1 ; i < n ; i++ )
p[i] = f[n-i];
for ( i = 0 ; i < n ; i++ )
{
p[i].re /= n;
p[i].im /= n;
}
}

void TDipWindow::tfft(int m,int n)
{
int i,j;

root(n);
p = (COMPLEX *) new COMPLEX[n];
f = (COMPLEX *) new COMPLEX[n];

for ( i = 0 ; i < m ; i++ )
{
for ( j = 0 ; j < n ; j++ )
p[j] = SpaceField[i][j];

fft(n);

for ( j = 0 ; j < n ; j++ )
FreqField[i][j] = f[j];
}

delete p;
delete f;

for ( i = 0 ; i < m ; i++ )
for ( j = 0 ; j < n ; j++ )
SpaceField[i][j] = FreqField[i][j];

root(m);
p = (COMPLEX *) new COMPLEX[m];
f = (COMPLEX *) new COMPLEX[m];

for ( j = 0 ; j < n ; j++ )
{
for ( i = 0 ; i < m ; i++ )
p[i] = SpaceField[i][j];

fft(m);

for ( i = 0 ; i < m ; i++ )
FreqField[i][j] = f[i];
}

delete p;
delete f;
}

void TDipWindow::ntfft(int m,int n)
{
int i,j;

root(n);
p = (COMPLEX *) new COMPLEX[n];
f = (COMPLEX *) new COMPLEX[n];

for ( i = 0 ; i < m ; i++ )
{
for ( j = 0 ; j < n ; j++ )
f[j] = FreqField[i][j];

nfft(n);

for ( j = 0 ; j < n ; j++ )
SpaceField[i][j] = p[j];
}

delete p;
delete f;

for ( i = 0 ; i < m ; i++ )
for ( j = 0 ; j < n ; j++ )
FreqField[i][j] = SpaceField[i][j];

root(m);
p = (COMPLEX *) new COMPLEX[m];
f = (COMPLEX *) new COMPLEX[m];

for ( j = 0 ; j < n ; j++ )
{
for ( i = 0 ; i < m ; i++ )
f[i] = FreqField[i][j];

nfft(m);

for ( i = 0 ; i < m ; i++ )
SpaceField[i][j] = p[i];
}

delete p;
delete f;
}

[尹成]

另外在FFT变换的时候：
X轴代表频率，它的分辨率与你进行FFT变换的点数和采样频率有关，设点数为N,采样频率为FS,则FFT变换后的结果（模）中，第0到N/2个数据代表了对应从0Hz到FS/2的功率谱，所以结果数据的分辨率为FS/N
Y轴代表功率，可用dB为单位，取个合适的比例系数表示高度就可以了。

[vctuxiangchuli]

谢谢楼主们

[xz801025]

不错啊，学习了

[wangyingyingqq]

友情帮顶

[vctuxiangchuli]

结贴，懂了

[lqa876270]

您好，我也在找FFT,也给我发一份。[e-mail]
谢谢！