用 nlinfit()拟合函数,来拟合出系数a,b,c,w.。
创新互联建站长期为近千家客户提供的网站建设服务,团队从业经验10年,关注不同地域、不同群体,并针对不同对象提供差异化的产品和服务;打造开放共赢平台,与合作伙伴共同营造健康的互联网生态环境。为蓬安企业提供专业的网站建设、网站制作,蓬安网站改版等技术服务。拥有10余年丰富建站经验和众多成功案例,为您定制开发。
x=[]
y=[]
myfun=inline('a(1)+a(2)*sin(a(3)*x)+a(4)*cos(a(3)*x)','a''x')
[a,r,J]=nlinfit(x,y,myfun,rand(1,4));
float ar[1024],ai[1024];/* 原始数据实部,虚部 */
float a[2050];
void fft(int nn) /* nn数据长度 */
{
int n1,n2,i,j,k,l,m,s,l1;
float t1,t2,x,y;
float w1,w2,u1,u2,z;
float fsin[10]={0.000000,1.000000,0.707107,0.3826834,0.1950903,0.09801713,0.04906767,0.02454123,0.01227154,0.00613588,};
float fcos[10]={-1.000000,0.000000,0.7071068,0.9238796,0.9807853,0.99518472,0.99879545,0.9996988,0.9999247,0.9999812,};
switch(nn)
{
case 1024: s=10; break;
case 512: s=9; break;
case 256: s=8; break;
}
n1=nn/2; n2=nn-1;
j=1;
for(i=1;i=nn;i++)
{
a[2*i]=ar[i-1];
a[2*i+1]=ai[i-1];
}
for(l=1;ln2;l++)
{
if(lj)
{
t1=a[2*j];
t2=a[2*j+1];
a[2*j]=a[2*l];
a[2*j+1]=a[2*l+1];
a[2*l]=t1;
a[2*l+1]=t2;
}
k=n1;
while (kj)
{
j=j-k;
k=k/2;
}
j=j+k;
}
for(i=1;i=s;i++)
{
u1=1;
u2=0;
m=(1i);
k=m1;
w1=fcos[i-1];
w2=-fsin[i-1];
for(j=1;j=k;j++)
{
for(l=j;lnn;l=l+m)
{
l1=l+k;
t1=a[2*l1]*u1-a[2*l1+1]*u2;
t2=a[2*l1]*u2+a[2*l1+1]*u1;
a[2*l1]=a[2*l]-t1;
a[2*l1+1]=a[2*l+1]-t2;
a[2*l]=a[2*l]+t1;
a[2*l+1]=a[2*l+1]+t2;
}
z=u1*w1-u2*w2;
u2=u1*w2+u2*w1;
u1=z;
}
}
for(i=1;i=nn/2;i++)
{
ar[i]=4*a[2*i+2]/nn; /* 实部 */
ai[i]=-4*a[2*i+3]/nn; /* 虚部 */
a[i]=4*sqrt(ar[i]*ar[i]+ai[i]*ai[i]); /* 幅值 */
}
}
给你改了下,自己看看吧,以后注意各种运算符的应用,及函数的使用
#includestdio.h
#includemath.h
main()
{
float x[1000],s[100],t;
double pi=3.1415926;
int n;
int E=2,w=1;
scanf("%f",t);
s[0]=0;
for(n=1;n=100;n++)
{
x[n]=2*E/pi*(1/(2*n-1)*cos((2*n-1)*w*t)*pow(-1.,n-1));
s[n]=x[n]+s[n-1];
}
printf("%f",s[n]);
}
#includestdio.h
#include math.h
class complex //定义一个类,实现复数的所有操作
{
double Real,Image; //实部与虚部
public:
complex(double r="0",double i="0"){Real=r;Image=i;}
double GetR(){return Real;} //取出实部
double GetI(){return Image;} //取出虚部
complex operator + (complex ); //复数加法
complex operator - (complex ); //复数减法
complex operator * (complex ); //复数乘法
void operator =(complex ); //复数 赋值
};
complex complex::operator + (complex c) //复数加法
{
complex t;
t.Real=Real+c.Real;
t.Image=Image+c.Image;
return t;
}
complex complex::operator - (complex c) //复数减法
{
complex t;
t.Real=Real-c.Real;
t.Image=Image-c.Image;
return t;
}
complex complex::operator * (complex c) //复数乘法
{
complex t;
t.Real=Real*c.Real-Image*c.Image;
t.Image=Real*c.Image+Image*c.Real;
return t;
}
void complex::operator = (complex c) //复数 赋值
{
Real=c.Real;
Image=c.Image;
}
void fft(complex a[],int length,int jishu) //实现fft的函数
{
const double PI="3".141592653589793;
complex u,Wn,t;
int i,j,k,m,kind,distance,other;
double tmp;
for(i=0;ilength;i++) //实现倒叙排列
{
k="i";
j=0;
for(m=0;mjishu;m++)
{
j="j"*2+k%2;
k/=2;
}
if(ij)
{
t="a";
a=a[j];
a[j]=t;
}
}
for(m=1;m=jishu;m++) //第m级蝶形运算,总级数为jishu
{
kind = (int)pow(2,m-1); //第m级有2^(m-1)种蝶形运算
distance = 2*kind; //同种蝶形结相邻距离为2^m
u=complex(1,0); //旋转因子初始值为 1
tmp=PI/kind;
Wn=complex(cos(tmp),-sin(tmp));//旋转因子Wn
for(j=0;jkind;j++) //每种蝶形运算的起始点为j,共有kind种
{
for(i=j;ilength;i+=distance) //同种蝶形运算
{
other=i+kind;//蝶形运算的两个因子对应单元下标的距离为2^(m-1)
t=a[other]*u; // 蝶形运算的乘积项
a[other]=a-t; //蝶形运算
a=a+t; //蝶形运算
}
u="u"*Wn; //修改旋转因子,多乘一个基本DFT因子WN
}
}
}
void main(void)
{
double a,b;
complex x[8]; //此程序以8点序列测试
printf("8点序列:\n");
for(int i="0";i8;i++) //初始化并输出原始序列
{
x=complex(i,i+1);
printf("x(%d) = %lf + %lf i\n",i+1,x.GetR(),x.GetI());
}
fft(x,8,3); //调用fft函数
printf("fft变换的结果为:\n");
for(i=0;i8;i++) //输出结果
printf("X(%d)= %lf + %lf i\n",i+1,x.GetR(),x.GetI());
}