傅里叶级数系数 的计算

http://blog.163.com/echo_8404@126/blog/#m=0&t=1&c=fks_087069080094086074080081085095085084088065083095085071086

傅里叶级数系数 的计算

2009-10-30 20:56:03| 分类： 数值计算方法|字号 订阅

/*

f(t) = a0+E(an*cos(n*w1*t)+bnsin(n*w1*t))

a0 = 1/T*积分（f(t)*dt）

an = 2/T*积分（f(t)*cos(n*w1)*dt)?

an = 2/T*积分（f(t)*sin(n*w1)*dt）

//********************************************************************************

文件名 : Fourier.c

作者 ： Signal And System 马

初次时间 : 09.10.30 晚上 20：43

描述 :

a,b ---- 1个周期内信号有值区间的下、上限。

t ------ 信号周期

l--------谐波次数

m--------子区间组数

cn-------n次谐波振幅

w ------基波角频率

an，bn---三角式系数

h ------步长

nw ------n次谐波

效果 ： 成功！！！

屏幕上打印出三角式系数an，bn，cn

和离散的幅度频谱图

//********************************************************************************

*/

//周期方波信号 傅里叶基数系数 计算

#include"stdio.h"

#include"stdlib.h"

#include"math.h"

main()

{

float ta[20][3];

float a = 0;

float b = 1;

float t = 4;

float pi = 3.1415926;

float x,a0,c0,pc,ps,an,bn,cn,h,w,t0,nw;

int i,j;

int n1,n;

int l=16,m = 16;

float san();

float sbn();

void gw12();

h = (b - a)/(2*m);

w = 2*pi/t;

t0 = a;

//求取直流分量幅度

nw = 0;

pc = san(t0,h,m,nw,pc);

a0 = pc/t;

c0 = a0;

ta[1][3] = c0*2;

printf("\na0=%4.2f,c0=%4.2f\n",a0,c0);

getch();

//求N次谐波分量频谱幅度

for(n=1;n<=l+1;n++)

{

nw = n*w;

pc = san(t0,h,m,nw,pc);

an = 2*pc/t;

ps = sbn(t0,h,m,nw,ps);

bn = 2*ps/t;

x = an*an + bn*bn;

cn = sqrt(x);

n1 = n+1;

ta[n1][1] = an;

ta[n1][2] = bn;

ta[n1][3] = cn;

}

//打印计算结果

printf("\n N= AN= BN= CN=");

for(j=2;j<=l+1;j++)

{

n = j-1;

printf("\n%4d%7.4f%7.4f%7.4f",n,ta[j][1],ta[j][2],ta[j][3]);

}

getch();

gw12(ta,l);

getch();

}

//求三角式系数AN的函数

float san(t0,h,m,nw,pc)

float t0,h,nw,pc;

int m;

{

float y[4],sm,tc,t;

int i,j;

sm = 0;

tc = t0;

for(j=1;j<=m;j++)

{

t = tc;

for(i=1;i<=3;i++)

{

y[i] = cos(nw*t);

t = tc+i*h;

}

sm = sm+y[1]+4*y[2]+y[3];

tc=tc+2*h;

}

pc=sm*h/3;

return(pc);

}

//求三角式系数BN的函数

float sbn(t0,h,m,nw,ps)

float t0,h,nw,ps;

int m;

{

float y[4],sm,tc,t;

int i,j;

sm = 0;

tc = t0;

for(j=1;j<=m;j++)

{

t = tc;

for(i=1;i<=3;i++)

{

y[i] = sin(nw*t);

t = tc+i*h;

}

sm = sm+y[1] + 4*y[2]+y[3];

tc=tc+2*h;

}

ps=sm*h/3;

return(ps);

}

//打印幅度频谱的函数--------------------------------------------------------------

void gw12(ta,l)

float ta[20][3];

int l;

{

float gm,sm,bg,ccn[20];

int i,j,n1,nx;

char ch[45];

for(j=0;j<=45;j++)

ch[j] = '*';

for(j=1;j<=l+1;j++)

{ccn[j] = ta[j][3]/2;}

//求最大最小值---------------------------------------------------------------

bg = ccn[1];

sm = ccn[1];

for(j=2;j<=l+1;j++)

{

if(ccn[j]>bg)

bg=ccn[j];

else if(ccn[j]<sm)

sm = ccn[j];

}

printf("\n\n min=%6.4fmax=%6.4f",sm,bg);//small,big

//--求出比例---------------------------------------------------------------------------

gm = bg -sm;

for(j=1;j<=l+1;j++)

ccn[j] = (ccn[j]-sm)/gm*44;

//打印负轴频谱-----------------------------------------------------------------------

printf("\n\nCN/2N");

for(j=1;j<=l+1;j++)

{

n1 = l+1-j+1;

nx = ccn[n1]+1;

printf("\n%6.4f%3d",ta[n1][3]/2,-n1+1);

for(i=1;i<=nx;i++)

putchar(ch[i]);

}

getch();

//--打印正轴频谱----------------------------------------------------------------

for(j=2;j<=l+1;j++)

{

n1 = j;

nx = ccn[n1]+1;

printf("\n%6.4f%3d",ta[n1][3]/2,n1-1);

for(i=1;i<=nx;i++)

putchar(ch[i]);

}

}