如何提高频谱分辨率
对密集信号在不增加采用长度如何提高频谱分辨率！找到一篇文章，仿真了下没效果，大家帮忙看看
function = hrft(N, q, Omega)
CO=cos(Omega);
SO=sin(Omega);
T0=1.0;
T1=CO;
U0=1.0;
U1=2.0*CO;
Qr=q(1)+q(2)*T1;
Qi=q(2)*SO*U0;
CO = CO*2.0;
for n=3:N
T2=CO*T1-T0;
T0=T1;
T1=T2;
U2=CO*U1-U0;
U0=U1;
U1=U2;
Cn=T2;
Sn=SO*U0;
Qr = Qr+q(n)*Cn;
Qi = Qi+q(n)*Sn;
end C版本#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#define MAXLINESIZE 128
/*Calculates the Fourier transform over Nf points from frequency f1 to f2,of N(N>1) samples of
a signal taken at sps sample per second.The values of cos(n*Omega)&sin(n*Omega) are calculated
using Chebyshev polynomials.*/
void hrft(int N,double *q,double Omega,double *Qr,double *Qi)
{
int n;
double C0,S0;
double T0,T1,T2;
double U0,U1,U2;
double Cn,Sn;
C0=cos(Omega);
S0=sin(Omega);
T0=1.0;
T1=C0;
U0=1.0;
U1=2.0*C0;
*Qr=q+q*T1;
*Qi=q*S0*U0;
C0 *=2.0;
for(n=2;n<N;++n)
{
T2=C0*T1-T0;
T0=T1;
T1=T2;
U2=C0*U1-U0;
U0=U1;
U1=U2;
Cn=T2;
Sn=S0*U0;
*Qr +=q*Cn;
*Qi +=q*Sn;
}
}
int main(int argc,char *argv[])
{
double f1,f2;
int Nf;
double df;
int N;
int i,n;
double sps;
double Omega;
double dOmega;
FILE *fp;
char linebuff;
double *q;
double Qr,Qi;
if (argc !=8)
{
printf ("\nhrft calculates the Fourier transform at Nf equally\n");
printf ("spaced points in the frequency interval \n");
printf ("usage: hrft f1 f2 Nf N sps infile outfile\n");
printf (" f1=starting frequency \n");
printf (" f2=ending frequency \n");
printf (" Nf=number of frequengcy points to calculate FT\n");
printf (" N=number of samples to read from the input file\n");
printf (" sps=samples per second\n");
printf (" infile=input file\n");
printf ("outfile=output file\n");
exit (-1);
}
f1=atof(argv);
f2=atof(argv);
Nf=atoi(argv);
N=atoi(argv);
sps=atof(argv);
if (N<2)
{
printf("N must be greater than 1\n");
exit (-1);
}
if ((fp=fopen(argv,"r"))==NULL)
{
perror("Error opening input file");
exit (-1);
}
printf("Reading data file: %s\n\n",argv);
for(n=0;n<N;++n)
{
fgets (linebuff,MAXLINESIZE,fp);
if (linebuff != '#') break;
printf("%s",linebuff);
}
q=(double *)malloc(N*sizeof(double));
sscanf (linebuff,"%lf",&q);
for (n=0;n<N;++n)
{
fscanf(fp,"%lf",&q);
}
fclose(fp);
if ((fp=fopen(argv,"w"))==NULL)
{
perror("Error opening output file");
exit (-1);
}
printf ("Opening output file:%s\n",argv);
df=(f2-f1)/((double)(Nf-1));
fprintf(fp,"# This file produced by program htft\n");
fprintf(fp,"# %20f :f1,starting frequency \n",f1);
fprintf(fp,"# %20f :f2,ending frequency \n",f2);
fprintf(fp,"# %20f :df,frequency step size \n",df);
fprintf(fp,"# %20d :Nf,number of frequency points at which FT was calculated\n",Nf);
fprintf(fp,"# %20d :N,number of samp;es read from the input file,N>1\n",N);
fprintf(fp,"# %20f :sps,samples per second\n",sps);
fprintf(fp,"# %s :input file name\n",argv);
fprintf(fp,"# Left column:real part of FT,Right column:imaginary part\n");
Omega=2.0*M_PI*f1/sps;
dOmega=2.0*M_PI*df/sps;
for (i=0;i<Nf;++i)
{
hrft (N,q,Omega,&Qr,&Qi);
fprintf (fp,"%lf %lf\n",Qr,Qi);
Omega +=dOmega;
}
fclose(fp);
free(q);
return(0);
} 人气呢！！！！
页:
[1]