哪位兄弟有c\c++的互相关的程序！！！

xueyongan405

先谢谢了

hyl2323

自己编写，不难的

风花雪月

1. #include <math.h>
   
2. #define M_PI    3.14159265358979323846
   
3. #define FALSE    0
   
4. #define TRUE    1
   
5. #define BIG    1e10
   
6. #define SMALL    1e-10
   
7. 
   
8. typedef struct {
   
9.     float r, i;
   
10. } complex;
    
11. 
    
12. 
    
13. 
    
14. /* FAST CORRELATION OF X(0:L) AND Y(0:L).  FINDS RXY(0) THRU RXY(NMAX). */
    
15. /* L=LAST INDEX IN BOTH X AND Y.  MUST BE (POWER OF 2)+1 AND AT LEAST 5. */
    
16. /* ITYPE=TYPE OF CORRELATION=0 IF X AND Y ARE THE SAME VECTOR (AUTO- */
    
17. /*         CORRELATION), OR NOT 0 IF X AND Y ARE DIFFERENT VECTORS. */
    
18. /* NMAX=MAXIMUM LAG OF INTEREST IN THE CORRELATION FUNCTION. */
    
19. /* FFT LENGTH ,N, USED INTERNALLY, IS L-1. */
    
20. /* LET K=INDEX OF FIRST NONZERO SAMPLE IN Y(0)---Y(N-1).  THEN X(0) */
    
21. /*  到 X(N-1) MUST INCLUDE PADDING OF AT LEAST NMAX-K ZEROS. */
    
22. /* CORRELATION FUNCTION, RXY, REPLACES X(0) THRU X(NMAX). */
    
23. /* Y(0) THRU Y(L) IS REPLACED BY ITS FFT, COMPUTED USING SPFFTR. */
    
24. /* IERROR=0  NO ERROR DETECTED */
    
25. /*        1  L-1 NOT A POWER OF 2 */
    
26. /*        2  NMAX OUT OF RANGE */
    
27. /*        3  INADEQUATE ZERO  */
    
28. 
    
29. 
    
30. void spcorr(float *x, float *y, long *l, long *type, long *nmax, long *error)
    
31. /*
    
32. x:序列X；
    
33. y：序列Y；
    
34. l:序列X与序列Y的长度，不小5，且要为2的幂次方；
    
35. type:相关的类型，0：表示X与Y序列相同，其它值：X与Y序列不相同
    
36. nmax:相关的最大时延；
    
37. error:运行出错提示；0：无错；1：数据长度不是2的幂次方；2：时延超界；3：无足够零填充出错
    
38. */
    
39. 
    
40. {
    
41.      long j, k, m, n;//n:FFT长度；k:序列Y中的首个非零样本的位置序号；在序列Y中必须最少包含有(nmax-k)零填充。
    
42.     complex cx;
    
43.     float test;
    
44. 
    
45.     n = *l - 1;
    
46.     if (*nmax < 0 || *nmax >= n)
    
47.     {
    
48.     *error = 2;
    
49.     return;
    
50.     }
    
51. 
    
52.     test = (float) n;
    
53.     test /= 2.0;
    
54. 
    
55.     while ((test - 2.0) > 0.0)
    
56.     {
    
57.     test /= 2.0;
    
58.     }
    
59. 
    
60.     if ((test - 2.0) == 0)
    
61.     {
    
62.     for (k = 0 ; k < n && y[k] == 0.0 ; ++k) ;
    
63. 
    
64.     for (j = n - 1 ; j >= 0 && x[j] == 0.0 ; --j) ;
    
65. 
    
66.     if ((n - 1 - j) < (*nmax - k))
    
67.     {
    
68.         *error = 3;
    
69.         return;
    
70.     }
    
71. 
    
72.     spfftr(x, &n);//对X序列FFT变换
    
73.     if (*type != 0)
    
74.     {
    
75.         spfftr(y, &n);//如果X、Y是相同序列，则对Y序列也进行FFT
    
76.     }
    
77. 
    
78.     for (m = 0 ; m <= (n / 2) ; ++m)
    
79.     {
    
80.         cx.r = x[m * 2] * y[m * 2] - -x[(m * 2) + 1] * y[(m * 2) + 1];
    
81.         cx.i = x[m * 2] * y[(m * 2) + 1] + -x[(m * 2) + 1] * y[m * 2];
    
82. 
    
83.         x[m * 2] = cx.r / n;
    
84.         x[(m * 2) + 1] = cx.i / n;
    
85.     }
    
86. 
    
87.     spiftr(x, &n);
    
88. 
    
89.     *error = 0;
    
90.     }
    
91.     else if ((test - 2.0) < 0.0)
    
92.     {
    
93.     *error = 1;
    
94.     }
    
95. 
    
96.     return;
    
97. } /* spcorr */
    
98. 
    
99. 
    
100. /* SPFFTR     11/12/85 */
     
101. /* FFT ROUTINE FOR REAL TIME SERIES (X) WITH N=2**K SAMPLES. */
     
102. /* COMPUTATION IS IN PLACE, OUTPUT REPLACES INPUT. */
     
103. /* INPUT:  REAL VECTOR X(0:N+1) WITH REAL DATA SEQUENCE IN FIRST N */
     
104. /*         ELEMENTS; ANYTHING IN LAST 2.  NOTE: X MAY BE DECLARED */
     
105. /*         REAL IN MAIN PROGRAM PROVIDED THIS ROUTINE IS COMPILED  */
     
106. /*         SEPARATELY ... COMPLEX OUTPUT REPLACES REAL INPUT HERE. */
     
107. /* OUTPUT: COMPLEX VECTOR XX(O:N/2), SUCH THAT X(0)=REAL(XX(0)),X(1)= */
     
108. /*         IMAG(XX(0)), X(2)=REAL(XX(1)), ..., X(N+1)=IMAG(XX(N/2). */
     
109. /* IMPORTANT:  N MUST BE AT LEAST 4 AND MUST BE A POWER OF 2. */
     
110. 
     
111. //FFT计算函数
     
112. void spfftr(complex *x, long *n)
     
113. {
     
114.     /* Builtin functions */
     
115.     void r_cnjg();
     
116. 
     
117.     /* Local variables */
     
118.     void spfftc();
     
119. 
     
120.     long m, tmp_int;
     
121.     complex u, tmp, tmp_complex;
     
122.     float tpn, tmp_float;
     
123. 
     
124.     tpn = (float) (2.0 * M_PI / (double) *n);
     
125. 
     
126.     tmp_int = *n / 2;
     
127.     spfftc(x, &tmp_int, &neg_i1);
     
128. 
     
129.     x[*n / 2].r = x[0].r;
     
130.     x[*n / 2].i = x[0].i;
     
131. 
     
132.     for (m = 0 ; m <= (*n / 4) ; ++m)
     
133.     {
     
134.     u.r = (float) sin((double) m * tpn);
     
135.     u.i = (float) cos((double) m * tpn);
     
136. 
     
137.     r_cnjg(&tmp_complex, &x[*n / 2 - m]);
     
138. 
     
139.     tmp.r = (((1.0 + u.r) * x[m].r - u.i * x[m].i)
     
140.         + (1.0 - u.r) * tmp_complex.r - -u.i * tmp_complex.i) / 2.0;
     
141. 
     
142.     tmp.i = (((1.0 + u.r) * x[m].i + u.i * x[m].r)
     
143.         + (1.0 - u.r) * tmp_complex.i + -u.i * tmp_complex.r) / 2.0;
     
144. 
     
145.     tmp_float = ((1.0 - u.r) * x[m].r - -u.i * x[m].i
     
146.             + (1.0 + u.r) * tmp_complex.r - u.i * tmp_complex.i) / 2.0;
     
147.     x[m].i = ((1.0 - u.r) * x[m].i + -u.i * x[m].r
     
148.          + (1.0 + u.r) * tmp_complex.i + u.i * tmp_complex.r) / 2.0;
     
149.     x[m].r = tmp_float;
     
150. 
     
151.     r_cnjg(&x[*n / 2 - m], &tmp);
     
152.     }
     
153. 
     
154.     return;
     
155. } /* spfftr */
     
156. 
     
157. 
     
158. /* SPIFTR     02/20/87 */
     
159. /* INVERSE FFT OF THE COMPLEX SPECTRUM OF A REAL TIME SERIES. */
     
160. /* X AND N ARE THE SAME AS IN SPFFTR.  IMPORTANT: N MUST BE A POWER */
     
161. /* OF 2 AND X MUST BE DIMENSIONED X(0:N+1) (REAL ARRAY, NOT COMPLEX). */
     
162. /* THIS ROUTINE TRANSFORMS THE OUTPUT OF SPFFTR BACK INTO THE INPUT, */
     
163. /* SCALED BY N.  COMPUTATION IS IN PLACE, AS IN SPFFTR. */
     
164. 
     
165. //逆FFT变换函数
     
166. void spiftr(complex *x, long *n)
     
167. {
     
168.     long m, tmp_int;
     
169.     complex u, tmp_complex, tmp;
     
170.     float tpn, tmp_float;
     
171. 
     
172.     tpn = (float) (2.0 * M_PI / (double) *n);
     
173. 
     
174.     for (m = 0 ; m <= (*n / 4) ; ++m)
     
175.     {
     
176.     u.r = (float) sin((double) m * tpn);
     
177.     u.i = (float) -cos((double) m * tpn);
     
178. 
     
179.     r_cnjg(&tmp_complex, &x[*n / 2 - m]);
     
180. 
     
181. 
     
182.     tmp.r = ((1.0 + u.r) * x[m].r - u.i * x[m].i)
     
183.         + ((1.0 - u.r) * tmp_complex.r - -u.i * tmp_complex.i);
     
184.     tmp.i = ((1.0 + u.r) * x[m].i + u.i * x[m].r)
     
185.         + ((1.0 - u.r) * tmp_complex.i + -u.i * tmp_complex.r);
     
186. 
     
187.     r_cnjg(&tmp_complex, &x[*n / 2 - m]);
     
188. 
     
189.     tmp_float = ((1.0 - u.r) * x[m].r - -u.i * x[m].i)
     
190.             + ((1.0 + u.r) * tmp_complex.r - u.i * tmp_complex.i);
     
191.     x[m].i = ((1.0 - u.r) * x[m].i + -u.i * x[m].r)
     
192.         + ((1.0 + u.r) * tmp_complex.i + u.i * tmp_complex.r);
     
193. 
     
194.     x[m].r = tmp_float;
     
195. 
     
196.     r_cnjg(&x[*n / 2 - m], &tmp);
     
197.     }
     
198.     tmp_int = *n / 2;
     
199. 
     
200.     spfftc(x, &tmp_int, &pos_i1);
     
201. 
     
202.     return;
     
203. } /* spiftr *
     
204. 
     
205. 
     
206. 
     
207. void r_cnjg(complex *r, complex *z)
     
208. {
     
209.     r->r = z->r;
     
210.     r->i = -z->i;
     
211. }

复制代码

sweatwbf

很不错，谢了！！！！！！！！！！！！！１１