如何提高频谱分辨率

a.gain

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

2012-8-30 15:00 上传

点击文件名下载附件

a.gain

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

a.gain

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. 
    
73. if (N<2)
    
74. {
    
75. printf("N must be greater than 1\n");
    
76. exit (-1);
    
77. }
    
78. if ((fp=fopen(argv[6],"r"))==NULL)
    
79. {
    
80. perror("Error opening input file");
    
81. exit (-1);
    
82. }
    
83. printf("Reading data file: %s\n\n",argv[6]);
    
84. 
    
85. for(n=0;n<N;++n)
    
86. {
    
87. fgets (linebuff,MAXLINESIZE,fp);
    
88. if (linebuff[0] != '#') break;
    
89. printf("%s",linebuff);
    
90. }
    
91. q=(double *)malloc(N*sizeof(double));
    
92. 
    
93. sscanf (linebuff,"%lf",&q[0]);
    
94. 
    
95. for (n=0;n<N;++n)
    
96. {
    
97. fscanf(fp,"%lf",&q[n]);
    
98. }
    
99. fclose(fp);
    
100. 
     
101. if ((fp=fopen(argv[7],"w"))==NULL)
     
102. {
     
103. perror("Error opening output file");
     
104. exit (-1);
     
105. }
     
106. printf ("Opening output file:%s\n",argv[7]);
     
107. df=(f2-f1)/((double)(Nf-1));
     
108. 
     
109. 
     
110. fprintf(fp,"# This file produced by program htft\n");
     
111. fprintf(fp,"# %20f :f1,starting frequency [Hz]\n",f1);
     
112. fprintf(fp,"# %20f :f2,ending frequency [Hz]\n",f2);
     
113. fprintf(fp,"# %20f :df,frequency step size [Hz]\n",df);
     
114. fprintf(fp,"# %20d :Nf,number of frequency points at which FT was calculated\n",Nf);
     
115. fprintf(fp,"# %20d :N,number of samp;es read from the input file,N>1\n",N);
     
116. fprintf(fp,"# %20f :sps,samples per second\n",sps);
     
117. fprintf(fp,"# %s :input file name\n",argv[6]);
     
118. fprintf(fp,"# Left column:real part of FT,Right column:imaginary part\n");
     
119. 
     
120. Omega=2.0*M_PI*f1/sps;
     
121. dOmega=2.0*M_PI*df/sps;
     
122. for (i=0;i<Nf;++i)
     
123. {
     
124. hrft (N,q,Omega,&Qr,&Qi);
     
125. fprintf (fp,"%lf %lf\n",Qr,Qi);
     
126. Omega +=dOmega;
     
127. }
     
128. 
     
129. fclose(fp);
     
130. free(q);
     
131. return(0);
     
132. }

复制代码

a.gain

人气呢！！！！