声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 1264|回复: 3

[FFT] 如何提高频谱分辨率

[复制链接]
发表于 2012-8-30 15:01 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
对密集信号在不增加采用长度如何提高频谱分辨率!
找到一篇文章,仿真了下没效果,大家帮忙看看 specmag.pdf (143.8 KB, 下载次数: 15)
回复
分享到:

使用道具 举报

 楼主| 发表于 2012-8-30 15:01 | 显示全部楼层
function [Qr, Qi] = 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
 楼主| 发表于 2012-8-30 15:04 | 显示全部楼层
C版本
  1. #include <stdlib.h>
  2. #include <stdio.h>
  3. #include <string.h>
  4. #include <math.h>
  5. #define MAXLINESIZE 128
  6. /*Calculates the Fourier transform over Nf points from frequency f1 to f2,of N(N>1) samples of
  7. a signal taken at sps sample per second.The values of cos(n*Omega)&sin(n*Omega) are calculated
  8. using Chebyshev polynomials.*/
  9. void hrft(int N,double *q,double Omega,double *Qr,double *Qi)
  10. {
  11. int n;
  12. double C0,S0;
  13. double T0,T1,T2;
  14. double U0,U1,U2;
  15. double Cn,Sn;
  16. C0=cos(Omega);
  17. S0=sin(Omega);
  18. T0=1.0;
  19. T1=C0;
  20. U0=1.0;
  21. U1=2.0*C0;
  22. *Qr=q[0]+q[1]*T1;
  23. *Qi=q[1]*S0*U0;
  24. C0 *=2.0;
  25. for(n=2;n<N;++n)
  26. {
  27. T2=C0*T1-T0;
  28. T0=T1;
  29. T1=T2;
  30. U2=C0*U1-U0;
  31. U0=U1;
  32. U1=U2;
  33. Cn=T2;
  34. Sn=S0*U0;
  35. *Qr +=q[n]*Cn;
  36. *Qi +=q[n]*Sn;
  37. }
  38. }
  39. int main(int argc,char *argv[])
  40. {
  41. double f1,f2;
  42. int Nf;
  43. double df;
  44. int N;
  45. int i,n;
  46. double sps;
  47. double Omega;
  48. double dOmega;
  49. FILE *fp;
  50. char linebuff[MAXLINESIZE];
  51. double *q;
  52. double Qr,Qi;
  53. if (argc !=8)
  54. {
  55. printf ("\nhrft calculates the Fourier transform at Nf equally\n");
  56. printf ("spaced points in the frequency interval [f1,f2]\n");
  57. printf ("usage: hrft f1 f2 Nf N sps infile outfile\n");
  58. printf (" f1=starting frequency [Hz]\n");
  59. printf (" f2=ending frequency [Hz]\n");
  60. printf (" Nf=number of frequengcy points to calculate FT\n");
  61. printf (" N=number of samples to read from the input file\n");
  62. printf (" sps=samples per second\n");
  63. printf (" infile=input file\n");
  64. printf ("outfile=output file\n");
  65. exit (-1);
  66. }
  67. f1=atof(argv[1]);
  68. f2=atof(argv[2]);
  69. Nf=atoi(argv[3]);
  70. N=atoi(argv[4]);
  71. sps=atof(argv[5]);

  72. if (N<2)
  73. {
  74. printf("N must be greater than 1\n");
  75. exit (-1);
  76. }
  77. if ((fp=fopen(argv[6],"r"))==NULL)
  78. {
  79. perror("Error opening input file");
  80. exit (-1);
  81. }
  82. printf("Reading data file: %s\n\n",argv[6]);

  83. for(n=0;n<N;++n)
  84. {
  85. fgets (linebuff,MAXLINESIZE,fp);
  86. if (linebuff[0] != '#') break;
  87. printf("%s",linebuff);
  88. }
  89. q=(double *)malloc(N*sizeof(double));

  90. sscanf (linebuff,"%lf",&q[0]);

  91. for (n=0;n<N;++n)
  92. {
  93. fscanf(fp,"%lf",&q[n]);
  94. }
  95. fclose(fp);

  96. if ((fp=fopen(argv[7],"w"))==NULL)
  97. {
  98. perror("Error opening output file");
  99. exit (-1);
  100. }
  101. printf ("Opening output file:%s\n",argv[7]);
  102. df=(f2-f1)/((double)(Nf-1));


  103. fprintf(fp,"# This file produced by program htft\n");
  104. fprintf(fp,"# %20f :f1,starting frequency [Hz]\n",f1);
  105. fprintf(fp,"# %20f :f2,ending frequency [Hz]\n",f2);
  106. fprintf(fp,"# %20f :df,frequency step size [Hz]\n",df);
  107. fprintf(fp,"# %20d :Nf,number of frequency points at which FT was calculated\n",Nf);
  108. fprintf(fp,"# %20d :N,number of samp;es read from the input file,N>1\n",N);
  109. fprintf(fp,"# %20f :sps,samples per second\n",sps);
  110. fprintf(fp,"# %s :input file name\n",argv[6]);
  111. fprintf(fp,"# Left column:real part of FT,Right column:imaginary part\n");

  112. Omega=2.0*M_PI*f1/sps;
  113. dOmega=2.0*M_PI*df/sps;
  114. for (i=0;i<Nf;++i)
  115. {
  116. hrft (N,q,Omega,&Qr,&Qi);
  117. fprintf (fp,"%lf %lf\n",Qr,Qi);
  118. Omega +=dOmega;
  119. }

  120. fclose(fp);
  121. free(q);
  122. return(0);
  123. }
复制代码
 楼主| 发表于 2012-8-31 13:18 | 显示全部楼层
人气呢!!!!
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-29 08:48 , Processed in 0.063560 second(s), 21 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表