信号处理模块Word下载.docx
《信号处理模块Word下载.docx》由会员分享，可在线阅读，更多相关《信号处理模块Word下载.docx（74页珍藏版）》请在冰点文库上搜索。
f1=2*pi*f0/fs;
f2=F/T;
f3=pi*f2*（1/fs）*（1/fs）;
y=a*cos（f1*k+f3*k*k）;
returny;
intmain（）
doublea,f0,T,F,fs;
//n=50;
a=1.5;
f0=10.0;
T=50.0;
F=200.0;
FILE*fp;
=n;
{
x[i]=lfm（a,f0,T,F,fs,i）;
printf（"
%10.4f"
fp=fopen（"
lfm1.txt"
"
w+"
fprintf（fp,"
%f"
fclose（fp）;
3.方波
4.白噪声（高斯分布）
5.均匀分布随机信号
6.双曲调频（对数相位信号或者线性周期信号）
（HyperbolicFMsignal还没做，赶紧做啊）
7.CW脉冲（相当于取一小段0~T内的sin函数或者cos函数）
二．信号处理
1.首先是重头戏FFT:
#include<
stdio.h>
math.h>
doublefft（doublex[],doubley[],intn,intsign）
//x[]先x（i）放fft的实部,FFT后放X（k）的实部；
y[]先y（i）放fft的实部,FFT后放Y（k）的实部，n为2的整数次方，sign为1时做FFT,sign为-1时，做IFFT；
inti,j,k,l,m,n1,n2;
doublec,c1,e,s,s1,t,tr,ti;
j=1;
for（i=1;
16;
i++）//做n点的FFT，并且N<
2的16次方
m=i;
j=j*2;
if（j==n）break;
n1=n-1;
j=0;
n1;
i++）//进行码位倒置
if（i<
j）
tr=x[j];
ti=y[j];
x[j]=x[i];
y[j]=y[i];
x[i]=tr;
y[i]=ti;
}
k=n/2;
while（k<
j+1）
j=j-k;
k=k/2;
j=j+k;
n1=1;
for（l=1;
l<
=m;
l++）//FFT
n1=2*n1;
n2=n1/2;
e=3.14159265359/n2;
c=1.0;
s=0.0;
c1=cos（e）;
s1=-sign*sin（e）;
for（j=0;
j<
n2;
j++）
for（i=j;
n;
i+=n1）
k=i+n2;
tr=c*x[k]-s*y[k];
ti=c*y[k]+s*x[k];
x[k]=x[i]-tr;
y[k]=y[i]-ti;
x[i]=x[i]+tr;
y[i]=y[i]+ti;
t=c;
c=c*c1-s*s1;
s=t*s1+s*c1;
if（sign==-1）//IFFT的公式前面有1/n
x[i]/=n;
y[i]/=n;
#definen16
doublew[16]={1.0,0.0,-1.0,0.0,1.0,0.0,-1.0,0.0,1.0,0.0,-1.0,0.0,1.0,-1.0,0.0,1.0};
doublez[16]={0.0};
fft（w,z,n,1）;
printf（"
\nFFT\n"
for（i=0;
%10.7f+J%10.7f"
w[i],z[i]）;
如果序列x（n）是实数，那么其博立叶变换X（k）一般是复数，但其实部是偶对称，虚部是奇对称，即x（A）具有如下共扼对称性：
x（o）和x（付／2）都是实数，且有：
X（k）=X*（N-k）,1<
=k<
=（N/2-1）;
2.IIR滤波器：
（butterworth通带内平坦，但是过渡带衰减的慢，chebyshiev通带内有波动，但是过渡带衰减的快）
ifilt——整型变量。
滤波器的类型。
取值为1、2和3，分别对应切比雪夫、逆切比雪夫和巴持沃兹滤波器。
band——整型变量。
滤波器的通带形式。
取值为1、2、3和4，分别对应低遇、高通、带通和带阻滤波器。
ns——整型变量。
滤波器的n阶节数。
n——整型变量。
滤波器每节的阶数。
对于低通和高通滤波器，n＝2；
对于带通和带阻滤波器，n＝4。
对于低通和高通滤波器fc通带边界频率；
fr：
阻带边界频率。
对于带通和带阻滤波器flc：
通带下边界频率；
fhc：
通带上边界频率；
fls：
阻带下边界频率；
fhs：
阻带上边界频率。
stdlib.h>
voidiirbcf（intifilt,intband,intns,intn,doublef1,doublef2,doublef3,doublef4,doubledb,doubleb[],doublea[]）
intk;
doubleomega,lamda,epslon,fl,fh;
doubled[5],c[5];
voidchebyi（int,int,int,double,double*,double*）,chebyii（int,int,int,double,double,double*,double*）;
voidbwtf（int,int,int,double*,double*）;
doublecosh1（double）,warp（double）,bpsub（double,double,double）,omin（double,double）;
voidfblt（double*,double*,int,int,double,double,double*,double*）;
if（（band==1）||（band==4））fl=f1;
if（（band==2）||（band==3））fl=f2;
if（band<
=3）fh=f3;
if（band==4）fh=f4;
if（ifilt<
3）
switch（band）
{
case1:
case2:
{
omega=warp（f2）/warp（f1）;
break;
}
case3:
omega=omin（bpsub（warp（f1）,fh,fl）,bpsub（warp（f4）,fh,fl））;
case4:
omega=omin（1.0/bpsub（warp（f2）,fh,fl）,1.0/bpsub（warp（f3）,fh,fl））;
}
lamda=pow（10.0,（db/20.0））;
//?
?
epslon=lamda/cosh（2*ns*cosh1（omega））;
for（k=0;
k<
ns;
k++）
switch（ifilt）
chebyi（2*ns,k,4,epslon,d,c）;
chebyii（2*ns,k,4,omega,lamda,d,c）;
bwtf（2*ns,k,4,d,c）;
}
fblt（d,c,n,band,fl,fh,&
b[k*（n+1）+0],&
a[k*（n+1）+0]）;
staticdoublecosh1（doublex）
doublez;
z=log（x+sqrt（x*x-1.0））;
return（z）;
staticdoublewarp（doublef）
doublepi,z;
z=tan（pi*f）;
staticdoublebpsub（doubleom,doublefh,doublefl）
z=（om*om-warp（fh）*warp（fl））/（（warp（fh）-warp（fl））*om）;
staticdoubleomin（doubleom1,doubleom2）
doublez,z1,z2;
z1=fabs（om1）;
z2=fabs（om2）;
z=（z1<
z2）?
z1:
z2;
staticvoidbwtf（intln,intk,intn,doubled[],doublec[]）
doublepi,tmp;
d[0]=1.0;
c[0]=1.0;
d[i]=0.0;
c[i]=0.0;
tmp=（k+1）-（ln+1.0）/2.0;
if（tmp==0.0）//当n为奇数时，在n/2处系数为1/（p+1）
c[1]=1.0;
else
c[1]=-2.0*cos（（2*（k+1）+ln-1）*pi/（2*ln））;
c[2]=1.0;
staticvoidchebyi（intln,intk,intn,doubleep,doubled[],doublec[]）
doublepi,gam,omega,sigma;
gam=pow（（（1.0+sqrt（1.0+ep*ep））/ep）,1.0/ln）;
sigma=0.5*（1.0/gam-gam）*sin（（2*（k+1）-1）*pi/（2*ln））;
omega=0.5*（1.0/gam+gam）*cos（（2*（k+1）-1）*pi/（2*ln））;
if（（（ln%2）==1）&
&
（（k+1）==（ln+1）/2））//n为奇数，在n/2处的系数
d[0]=-sigma;
//分子
c[0]=d[0];
//分母
c[0]=sigma*sigma+omega*omega;
//分母
c[1]=-2.0*sigma;
d[0]=c[0];
//分子
if（（（ln%2）==0）&
（k==0））
d[0]=d[0]/sqrt（1.0+ep*ep）;
//前面的系数A
staticvoidchebyii（intln,intk,intn,doublews,doubleatt,doubled[],doublec[]）
doublepi,gam,alpha,beta,sigma,omega,scln,scld;
gam=pow（（att+sqrt（att*att-1.0））,1.0/ln）;
alpha=0.5*（1.0/gam-gam）*sin（（2*（k+1）-1）*pi/（2*ln））;
beta=0.5*（1.0/gam+gam）*cos（（2*（k+1）-1）*pi/（2*ln））;
sigma=ws*alpha/（alpha*alpha+beta*beta）;
omega=-1.0*ws*beta/（alpha*alpha+beta*beta）;
（（k+1）==（ln+1）/2））
d[0]=-1.0*sigma;
scln=sigma*sigma+omega*omega;
scld=pow（（ws/cos（（2*（k+1）-1）*pi/（2*ln）））,2）;
d[0]=scln*scld;
//分子常数项
d[2]=scln;
//分子平方项
c[1]=-2.0*sigma*scld;
c[2]=scld;
staticvoidfblt（doubled[],doublec[],intn,intband,doublefln,doublefhn,doubleb[],doublea[]）
inti,k,m,n1,n2,ls;
doublepi,w,w0,w1,w2,tmp,tmpd,tmpc,*work;
doublecombin（int,int）;
voidbilinear（double*,double*,double*,double*,int）;
w1=tan（pi*fln）;
//?
for（i=n;
i>
=0;
i--）
if（（c[i]!
=0.0）||（d[i]!
=0.0））
break;
switch（band）
case1:
case2:
n2=m;
n1=n2+1;
if（band==2）
for（i=0;
=m/2;
{
tmp=d[i];
d[i]=d[m-i];
d[m-i]=tmp;
tmp=c[i];
c[i]=c[m-i];
c[m-i]=tmp;
}
for（i=0;
d[i]=d[i]/pow（w1,i）;
c[i]=c[i]/pow（w1,i）;
case3:
case4:
n2=2*m;
work=（double*）malloc（n1*n1*sizeof（double））;
w2=tan（pi*fhn）;
w=w2-w1;
w0=w1*w2;
if（band==4）
=n2;
work[0*n1+i]=0.0;
work[1*n1+i]=0.0;
tmpd=d[i]*pow（w,（m-i））;
tmpc=c[i]*pow（w,（m-i））;
for（k=0;
=i;
ls=m+i-2*k;
tmp=combin（i,i）/（combin（k,k）*combin（i-k,i-k））;
work[0*n1+ls]+=tmpd*pow（w0,k）*tmp;
work[1*n1+ls]+=tmpc*pow（w0,k）*tmp;
d[i]=work[0*n1+i];
c[i]=work[1*n1+i];
free（work）;
bilinear（d,c,b,a,n）;
staticdoublecombin（inti1,inti2）
doubles;
s=1.0;
if（i2==0）
return（s）;
for（i=i1;
（i1-i2）;
s*=i;
return（s）;
staticvoidbilinear（doubled[],doublec[],doubleb[],doublea[],intn）
inti,j,n1;
doublesum,atmp,scale,*temp;
n1=n+1;
temp=（double*）malloc（n1*n1*sizeof（double））;
temp[j*n1+0]=1.0;
sum=1.0;
sum=sum*（double）（n-i+1）/（double）i;
temp[0*n1+i]=sum;
for（j=1;
temp[j*n1+i]=temp[（j-1）*n1+i]-temp[j*n1+i-1]-temp[（j-1）*n1+i-1];
b[i]=0.0;
atmp=0.0;
for（j=0;
b[i]=b[i]+temp[j*n1+i]*d[j];
atmp=atmp+temp[j*n1+i]*c[j];
scale=atmp;
if（i!
=0）a[i]=atmp;
b[i]=b[i]/scale;
a[i]=a[i]/scale;
a[0]=1.0;
free（temp）;
inti,k,n,ns,ifilt,band;
doublea[50],b[50];
doublef1,f2,f3,f4,fc,fr,fs,flc,fls,fhc,fhs,db;
enter1forchebyshevI,2forchebyshevII,3forbutterworth\n"
scanf（"
%d"
&
ifilt）;
enter1forlowpass,2forhighpass,3forbandpass,4forbandstop\n"
band）;
n=（band<
=2）?
2:
4;
enterthenumberoffiltersection\n"
ns）;
entersamplefrequencyfs\n"
%lf"
fs）;
=2）
printf（"
enterpassbandedgefrequencyfc\n"
scanf（"
fc）;
enterstopbandedgefrequencyfr\n"
fr）;
if（band=1）
f1=fc;
f2=fr;
else
f1=fr;
f2=fc;
f3=f4=0;
enterthelowerpassbandedgefrequencyflc\n"
flc）;
enterthehigherpassbandedgefrequencyfhc\n"
fhc）;
enter
升级会员 