/*
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]);
}
}
