下面是一个short time Fourier transform程序. 但许多地方我理解不了,而且执行结果好象还有错. 哪位"数学专业的研究生" 或 "信号处理专业的高手" 帮我讲解并注释一下.本人要用它来对麦克风采集的声音进行滤波去噪.提供相关资料也行. (关键是有怎样用代码实现Fourier transform的相关理论)#include <stdio.h>
#include <math.h>
#include <graphics.h>#define pi 3.1425926/* define the functions used throughout program */void fft(float* p1, float* p2, int n); /* function for computing FFT */
void curve(int i, float* p, int j); /* function for display signal */
void curve1(float* p, int s); /* function for display the SFTF results int time domain */int m, n;
float c[256][16], d[256][16];
float *a, *b, *w;main()
{
int i, j, k, window, s;
char name[30], wind[30];
int driver, mode;
float* p;
unsigned char p1[300];
FILE* fp; /* Read the input signal data file */ do
{
printf("please input the name of data file:");
scanf("%s", name);
}while((fp = fopen(name, "rb")) == NULL);
printf("please input the length of the data:");
scanf("%d", &m);
if((p = (float*)calloc(2 * m, sizeof(float))) == NULL)
{
printf("Can' t allocate memory to data n ");
exit(1);
} while((fscanf(fp, "%*f%f", p)) != EOF)
{++p;}
fclose(fp); /* Get the type of window to use on the time domain */ printf("please chose the style of window(1,2,3,4,5): n ");
printf(" 1. suare window n");
printf(" 2. Bartlett window n ");
printf(" 3. Hanning window n ");
printf(" 4. Hamming window n ");
printf(" 5. Blackman window n ");
printf(" 6. Gaussian window n ");
scanf("%d", &window);
if(window == 1) strcpy(wind, "squre window");
if(window == 2) strcpy(wind, "Bartlet window");
if(window == 3) strcpy(wind, "Hanning window");
if(window == 4) strcpy(wind, "Hamming window");
if(window == 5) strcpy(wind, "Blackman window");
if(window == 6) strcpy(wind, "Gaussian window");
scanf("%d", &n);
printf("please input the distance of the window shift:");
scanf("%d", &s);
if((a = (float*)calloc(2 * n, sizeof(float))) == NULL)
{
printf("Can't allocate memory to data n ");
exit(1);
}
if((b = (float*)calloc(2 * n, sizeof(float))) == NULL)
{
printf("Can't allocate memory to data n ");
exit(1);
}
if((w = (float*) calloc(m, sizeof(float))) == NULL)
{
printf("Can't allocate memory to data n ");
exit(1);
} /* perform the window function requested by the user */ switch(window)
{
case 1:
for(i=0; i <= n-1; i++)
w[i] = 1;
break;
case 2:
for(i=0; i <= (n-1)/2; i++)
w[i] = 2.0 * i /(n-1);
for(i=(n-1)/2+1;i <= n-1; i++)
w[i] = 2.0-2.0*i/(n-1);
break;
case 3:
for(i=0; i <= n-1; i++)
w[i] = (1.0 - cos(2.0 * pi * i/(n-1)))/2.0;
break;
case 4:
for(i=0; i <= n-1; i++)
w[i] = 0.54 - 0.46 * cos(2 * pi * i /(n-1));
break;
case 5:
for(i=0; i <= n-1; i++)
w[i] = 0.42 - 0.5 * cos(2 * pi * i / (n-1)) + 0.08 * cos(4 * pi * i / (n-1));
break;
case 6:
for(i=0; i <= n-1; i++)
w[i] = 1.0 / exp(18 * ((n-1)/2.0 - i) * ((n-1)/2.0 - i) / (n*n));
break;
default:
printf("error option!");
break;
} /* Initializing graphics */
driver = DETECT;
initgraph(&driver, &mode, ""); /* display the signal in graphics */
settextstyle(0, 0, 1);
moveto(1, 10);
outtext("original signal");
curve(m, p, 200);
moveto(1, 280);
outtext("window function");
for(i=0; i<n; i++)
*(w+i) = 100 * (*(w+i));
curve(m, w, 400);
getch();
cleardevice();
for(i=0; i<n; i++) *(w+i) = 0.01 * (*(w+i));
for(i=m; i<m+n; i++)
{
p[i] = p[i-m];
}
for(i=0; i<= n/2-1; i++)
{
a[i] = cos(2 * pi * i/n);
b[i] = sin(2 * pi * i/n);
} /* Find the spectrum of the windowed data */
for(i=0; i*s<m; i++)
{
for(j=0; j<n; j++)
c[i][j] = (*(p+s*i+j))*w[j];
fft(c[i],d[i], n);
}
if((fp = fopen("stft.data", "w")) == NULL) exit(1);
fprintf(fp, "the name of of data file: %s n ", name);
fprintf(fp, "the length of data: %sd n", m);
fprintf(fp, "the style of window: %s n", wind);
fprintf(fp, "the width of window: %d n", n);
fprintf(fp, "the distance of window shifting: %d n", s);
fprintf(fp, "STFT resultsq ");
fprintf(fp, "time frequency stft modular n q q ");
for(i=0; i*s<m; i++)
{
for(j=0; j<n/2; j++)
{
c[i][j] = sqrt(c[i][j] * c[i][j] + d[i][j] * d[i][j]);
fprintf(fp, "%d %d %f n ", i, j, c[i][j]);
}
} cleardevice();
curve1(p, s);
cleardevice();
closegraph();
}void fft(float* p1, float* p2, int n)
{
int i, j, k, power, w, cc;
int l;
float aa;
float t1, t2;
power = (int) (log10(n)/log10(2));
j = 1;
for(i=1; i<= n-1; i++)
{
if(i<j)
{
aa = *(p1+j-1); *(p1+j-1) = *(p1+i-1); *(p1+i-1) = aa;
aa = *(p2+j-1); *(p2+j-1) = *(p2+i-1); *(p2+i-1) = aa;
}
k = n/2;
while(1)
{
if(k>=j) break;
j = j-k;
k = k/2;
}
j = j+k;
}
l = 1;
for(i=1; i <= power; i++)
{
l = 2 * l;
for(j=1; j <= l/2; j++)
{
for(k = j-1; k <= n-1; k=k+1)
{
w = k+l/2;
cc = (j-1)*n/l;
t1 = (*(p1+w))*a[cc]-(*(p2+w))*b[cc];
t2 = (*(p1+w))*b[cc]+(*(p2+w))*a[cc];
*(p1+w) = *(p1+k) - t1; *(p2+w) = *(p2+k) - t2;
*(p1+k) = *(p1+k) + t2; *(p2+k) = *(p2+k) + t2;
}
}
}
}void curve(int i, float* p, int j)
{
int k, m1;
setcolor(RED);
line(10, j, 522, j);
for(k=0; k<i; k++)
{
m1 = 512/i;
setcolor(GREEN);
line(10+k*m1, j-*(p+k), 10+(k+1)*m1, j-*(p+k+1));
}
}void curve1(float*p, int s)
{
int i, j, k, kk;
float cc;
line(50, 1, 562, 1);
line(50, 450, 562, 450);
line(50, 1, 50, 450);
line(562, 1, 562, 450);
outtextxy(145, 455, "50");
outtextxy(245, 455, "100");
outtextxy(345, 455, "150");
outtextxy(445, 455, "200");
outtextxy(555, 455, "250");
k = n / 8;
for(kk=0; kk<k; kk++)
{
for(i=0; i<m; i++)
line(50+i, 100-p[i], 50+i+1, 100-p[i+1]);
for(i=0; i<4; i++)
{
line(50, 140+70*(i+1), 562, 140+70*(i+1));
for(j=0; j<m; j++)
line(50+s*j, 140+70*(i+1)-c[j][kk*4+i], 50 + s*(j+1), 140+70*(i+1)-c[j+1][kk*4+i]);
}
getch();
cleardevice();
}
}
#include <math.h>
#include <graphics.h>#define pi 3.1425926/* define the functions used throughout program */void fft(float* p1, float* p2, int n); /* function for computing FFT */
void curve(int i, float* p, int j); /* function for display signal */
void curve1(float* p, int s); /* function for display the SFTF results int time domain */int m, n;
float c[256][16], d[256][16];
float *a, *b, *w;main()
{
int i, j, k, window, s;
char name[30], wind[30];
int driver, mode;
float* p;
unsigned char p1[300];
FILE* fp; /* Read the input signal data file */ do
{
printf("please input the name of data file:");
scanf("%s", name);
}while((fp = fopen(name, "rb")) == NULL);
printf("please input the length of the data:");
scanf("%d", &m);
if((p = (float*)calloc(2 * m, sizeof(float))) == NULL)
{
printf("Can' t allocate memory to data n ");
exit(1);
} while((fscanf(fp, "%*f%f", p)) != EOF)
{++p;}
fclose(fp); /* Get the type of window to use on the time domain */ printf("please chose the style of window(1,2,3,4,5): n ");
printf(" 1. suare window n");
printf(" 2. Bartlett window n ");
printf(" 3. Hanning window n ");
printf(" 4. Hamming window n ");
printf(" 5. Blackman window n ");
printf(" 6. Gaussian window n ");
scanf("%d", &window);
if(window == 1) strcpy(wind, "squre window");
if(window == 2) strcpy(wind, "Bartlet window");
if(window == 3) strcpy(wind, "Hanning window");
if(window == 4) strcpy(wind, "Hamming window");
if(window == 5) strcpy(wind, "Blackman window");
if(window == 6) strcpy(wind, "Gaussian window");
scanf("%d", &n);
printf("please input the distance of the window shift:");
scanf("%d", &s);
if((a = (float*)calloc(2 * n, sizeof(float))) == NULL)
{
printf("Can't allocate memory to data n ");
exit(1);
}
if((b = (float*)calloc(2 * n, sizeof(float))) == NULL)
{
printf("Can't allocate memory to data n ");
exit(1);
}
if((w = (float*) calloc(m, sizeof(float))) == NULL)
{
printf("Can't allocate memory to data n ");
exit(1);
} /* perform the window function requested by the user */ switch(window)
{
case 1:
for(i=0; i <= n-1; i++)
w[i] = 1;
break;
case 2:
for(i=0; i <= (n-1)/2; i++)
w[i] = 2.0 * i /(n-1);
for(i=(n-1)/2+1;i <= n-1; i++)
w[i] = 2.0-2.0*i/(n-1);
break;
case 3:
for(i=0; i <= n-1; i++)
w[i] = (1.0 - cos(2.0 * pi * i/(n-1)))/2.0;
break;
case 4:
for(i=0; i <= n-1; i++)
w[i] = 0.54 - 0.46 * cos(2 * pi * i /(n-1));
break;
case 5:
for(i=0; i <= n-1; i++)
w[i] = 0.42 - 0.5 * cos(2 * pi * i / (n-1)) + 0.08 * cos(4 * pi * i / (n-1));
break;
case 6:
for(i=0; i <= n-1; i++)
w[i] = 1.0 / exp(18 * ((n-1)/2.0 - i) * ((n-1)/2.0 - i) / (n*n));
break;
default:
printf("error option!");
break;
} /* Initializing graphics */
driver = DETECT;
initgraph(&driver, &mode, ""); /* display the signal in graphics */
settextstyle(0, 0, 1);
moveto(1, 10);
outtext("original signal");
curve(m, p, 200);
moveto(1, 280);
outtext("window function");
for(i=0; i<n; i++)
*(w+i) = 100 * (*(w+i));
curve(m, w, 400);
getch();
cleardevice();
for(i=0; i<n; i++) *(w+i) = 0.01 * (*(w+i));
for(i=m; i<m+n; i++)
{
p[i] = p[i-m];
}
for(i=0; i<= n/2-1; i++)
{
a[i] = cos(2 * pi * i/n);
b[i] = sin(2 * pi * i/n);
} /* Find the spectrum of the windowed data */
for(i=0; i*s<m; i++)
{
for(j=0; j<n; j++)
c[i][j] = (*(p+s*i+j))*w[j];
fft(c[i],d[i], n);
}
if((fp = fopen("stft.data", "w")) == NULL) exit(1);
fprintf(fp, "the name of of data file: %s n ", name);
fprintf(fp, "the length of data: %sd n", m);
fprintf(fp, "the style of window: %s n", wind);
fprintf(fp, "the width of window: %d n", n);
fprintf(fp, "the distance of window shifting: %d n", s);
fprintf(fp, "STFT resultsq ");
fprintf(fp, "time frequency stft modular n q q ");
for(i=0; i*s<m; i++)
{
for(j=0; j<n/2; j++)
{
c[i][j] = sqrt(c[i][j] * c[i][j] + d[i][j] * d[i][j]);
fprintf(fp, "%d %d %f n ", i, j, c[i][j]);
}
} cleardevice();
curve1(p, s);
cleardevice();
closegraph();
}void fft(float* p1, float* p2, int n)
{
int i, j, k, power, w, cc;
int l;
float aa;
float t1, t2;
power = (int) (log10(n)/log10(2));
j = 1;
for(i=1; i<= n-1; i++)
{
if(i<j)
{
aa = *(p1+j-1); *(p1+j-1) = *(p1+i-1); *(p1+i-1) = aa;
aa = *(p2+j-1); *(p2+j-1) = *(p2+i-1); *(p2+i-1) = aa;
}
k = n/2;
while(1)
{
if(k>=j) break;
j = j-k;
k = k/2;
}
j = j+k;
}
l = 1;
for(i=1; i <= power; i++)
{
l = 2 * l;
for(j=1; j <= l/2; j++)
{
for(k = j-1; k <= n-1; k=k+1)
{
w = k+l/2;
cc = (j-1)*n/l;
t1 = (*(p1+w))*a[cc]-(*(p2+w))*b[cc];
t2 = (*(p1+w))*b[cc]+(*(p2+w))*a[cc];
*(p1+w) = *(p1+k) - t1; *(p2+w) = *(p2+k) - t2;
*(p1+k) = *(p1+k) + t2; *(p2+k) = *(p2+k) + t2;
}
}
}
}void curve(int i, float* p, int j)
{
int k, m1;
setcolor(RED);
line(10, j, 522, j);
for(k=0; k<i; k++)
{
m1 = 512/i;
setcolor(GREEN);
line(10+k*m1, j-*(p+k), 10+(k+1)*m1, j-*(p+k+1));
}
}void curve1(float*p, int s)
{
int i, j, k, kk;
float cc;
line(50, 1, 562, 1);
line(50, 450, 562, 450);
line(50, 1, 50, 450);
line(562, 1, 562, 450);
outtextxy(145, 455, "50");
outtextxy(245, 455, "100");
outtextxy(345, 455, "150");
outtextxy(445, 455, "200");
outtextxy(555, 455, "250");
k = n / 8;
for(kk=0; kk<k; kk++)
{
for(i=0; i<m; i++)
line(50+i, 100-p[i], 50+i+1, 100-p[i+1]);
for(i=0; i<4; i++)
{
line(50, 140+70*(i+1), 562, 140+70*(i+1));
for(j=0; j<m; j++)
line(50+s*j, 140+70*(i+1)-c[j][kk*4+i], 50 + s*(j+1), 140+70*(i+1)-c[j+1][kk*4+i]);
}
getch();
cleardevice();
}
}
complicated for me to understand.
Your words lost logic
不一定要看完上面的代码, 关键是能提供怎样编程实现Fourier transform的相关理论.对每一位热心肠的朋友都表示非常感谢!
建议你做如下工作：
1、用C实现FFT,用Matlab检查
2、用Matlab实现STFT
3、再用C实现STFT另外，你要给出一些描述啊，你用的是什么窗函数？如果你的要求不高，不如直接用Matlab吧，简单些。
--------------------------------------------------------------------------------
[本篇全文] [回复本文] [本篇作者: datou ] [本篇人气: 154]
发信人: datou (小七), 信区: DSP
标 题: 那里有fft快速傅立叶变换算法的c/c++源代码？
发信站: 浙江大学BBS (Sat Apr 20 14:39:16 2002), 转信谢谢--
※ 来源:·浙江大学BBS bbs.zju.edu.cn·[FROM: datou] --------------------------------------------------------------------------------
[本篇全文] [回复本文] [本篇作者: ALOHA ] [本篇人气: 75]
发信人: ALOHA (ALOHA), 信区: DSP
标 题: Re: 那里有fft快速傅立叶变换算法的c/c++源代?…
发信站: 浙江大学BBS (Sat Apr 20 21:08:14 2002), 站内信件
/*===================
fft中的整序函数，其中T为s的位数
===================*/
int sort(int s,int T)
{
int t,d;
for (d=0,t=0;t<T;t++)
{
d<<=1;
d+=s&SIEVE;
s>>=1;轻松的心情), 信区: DSP理 讨论区 [DSP
}
return d;
}/*==================
cbuff: 为复数的输入输出缓冲地址
iN : 序列长度，应为2的T次幂
Return: 1: Success
0: Out of Memory
2: iN不是T的整数次幂
==================*/
int fft(std::complex<double> cbuff[],int iN)
{
/*---------------------
Declaration & Allocate
---------------------*/
int it,ig,ii,T;
int *iSNo = new int[iN];
std::complex<double> ctemp;
std::complex<double> *cWN = new std::complex<double>[iN/2];
if (!(iSNo&&cWN))
return 0;傻男那?, 信区: DSP理 讨论区 [DSP
/*---------------------
Calculate T (仅限于iN为2的整数次幂)
---------------------*/
ii=iN; T=0;
do {
if (ii&1==1)
break;
ii>>=1;
T++;
}while(1);
if (1<<T!=iN)
return 2;
/*---------------------
Sorting
---------------------*/
for (ii=0;ii<iN;ii++)
{
iSNo[ii] = sort(ii,T);
if (iSNo[ii]>ii)
{
ctemp=cbuff[ii];
cbuff[ii =cbuff[iSNo[ii]]; 讨论区 [DSP
cbuff[iSNo[ii]]=ctemp;
}
}
/*---------------------
Prepare the cWN
---------------------*/
for(ii=0;ii<iN/2;ii++)
{
cWN[ii] = std::exp(std::complex<double>(0,-2*PAI*ii/iN));
}
/*---------------------
Main Procedure
---------------------*/
for (it=0;it<T;it++)
{
for (ig=0;ig<(1<<(T-it-1));ig++)
{
for (ii=0;ii<(1<<it);ii++)
{
ctemp = cWN[(1<<T-it-1)*ii] * cbuff[ig*(1<<it+
1)+ii+(1<<it)];
cbuff[ig*(1<<it+1)+ii+(1<<it)] = 讨论区 [DSP
cbuff[ig*(1<<it+1)+ii] - ctemp;
cbuff[ig*(1<<it+1)+ii] = cbuff[ig*(1<<it+1)+ii
] + ctemp;
}
}
}
/*---------------------
Free Memory
---------------------*/
delete iSNo;
delete cWN;
/*---------------------
Success
---------------------*/
return 1;
}--
DSP means Do Something Phenomenal
--
静女其姝 俟我于城隅 爱而不见 骚首趾蹰
静女其娈 贻我彤管 彤管有炜 说怿女美
自牧归荑 洵美且异 匪女之为美 美人之贻※ 来源:·浙江大学BBS bbs.zju.edu.cn·[FROM: ALOHA] --------------------------------------------------------------------------------
[本篇全文] [回复本文] [本篇作者: ALOHA ] [本篇人气: 46]
发信人: ALOHA (ALOHA), 信区: DSP
标 题: Re: 那里有fft快速傅立叶变换算法的c/c++源代?…
发信站: 浙江大学BBS (Sat Apr 20 21:08:58 2002), 站内信件 #include "math.h"
void kkfft(pr,pi,n,k,fr,fi,l,il)
int n,k,l,il;
double pr[],pi[],fr[],fi[];
{ int it,m,is,i,j,nv,l0;
double p,q,s,vr,vi,poddr,poddi;
for (it=0; it<=n-1; it++)
{ m=it; is=0;
for (i=0; i<=k-1; i++)
{ j=m/2; is=2*is+(m-2*j); m=j;}
fr[it]=pr[is]; fi[it]=pi[is];
}
pr[0]=1.0; pi[0]=0.0;
p=6.283185306/(1.0*n);
pr[1]=cos(p); pi[1]=-sin(p);
if (l!=0) pi[1]=-pi[1];
for (i=2; i<=n-1; i++)
{ p=pr[i-1]*pr[1]; q=pi[i-1]*pi[1];
s=(pr[i-1]+pi[i-1])*(pr[1]+pi[1]);
pr[i]=p-q; pi[i]=s-p-q;
}
for (it=0; it<=n-2; it=it+2)
{ vr=fr[it]; vi=fi[it];
fr[it]=vr+fr[it+1]; fi[it]=vi+fi[it+1];
fr[it+1]=vr-fr[it+1]; fi[it+1]=vi-fi[it+1];
}
m=n/2; nv=2;
for (l0=k-2; l0>=0; l0--)
{ m=m/2; nv=2*nv;
for (it=0; it<=(m-1)*nv; it=it+nv)
for (j=0; j<=(nv/2)-1; j++)
{ p=pr[m*j]*fr[it+j+nv/2];
q=pi[m*j]*fi[it+j+nv/2];
s=pr[m*j]+pi[m*j];
s=s*(fr[it+j+nv/2]+fi[it+j+nv/2]);
poddr=p-q; poddi=s-p-q;
fr[it+j+nv/2]=fr[it+j]-poddr;
fi[it+j+nv/2]=fi[it+j]-poddi;
fr[it+j]=fr[it+j]+poddr;
fi[it+j]=fi[it+j]+poddi;
}
}
if (l!=0)
for (i=0; i<=n-1; i++)
{ fr[i]=fr[i]/(1.0*n);
fi[i]=fi[i]/(1.0*n);
}
if (il!=0)
for (i=0; i<=n-1; i++)
{ pr[i]=sqrt(fr[i]*fr[i]+fi[i]*fi[i]);
if (fabs(fr[i])<0.000001*fabs(fi[i]))
{ if ((fi[i]*fr[i])>0) pi[i]=90.0;
else pi[i]=-90.0;
}
else
pi[i]=atan(fi[i]/fr[i])*360.0/6.283185306;
}
return;
}
--
静女其姝 俟我于城隅 爱而不见 骚首趾蹰
静女其娈 贻我彤管 彤管有炜 说怿女美
自牧归荑 洵美且异 匪女之为美 美人之贻※ 来源:·浙江大学BBS bbs.zju.edu.cn·[FROM: ALOHA] --------------------------------------------------------------------------------
[返回上一页][本讨论区]
另外毫不客气地指出，用“本姑娘有难”等字眼吸引眼球有时候真的很引人反感！试想一下，如果是个男的写“本少爷有难”等字句，你又会有何感想？
这样骗mm的qq倒是不错的嘛:))