1 // Copyright (c) <2012> <Leif Asbrink>
2 //
3 // Permission is hereby granted, free of charge, to any person
4 // obtaining a copy of this software and associated documentation
5 // files (the "Software"), to deal in the Software without restriction,
6 // including without limitation the rights to use, copy, modify,
7 // merge, publish, distribute, sublicense, and/or sell copies of
8 // the Software, and to permit persons to whom the Software is
9 // furnished to do so, subject to the following conditions:
10 //
11 // The above copyright notice and this permission notice shall be
12 // included in all copies or substantial portions of the Software.
13 //
14 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
15 // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
16 // OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
17 // NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
18 // HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
19 // WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
20 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE
21 // OR OTHER DEALINGS IN THE SOFTWARE.
22
23
24 #include "globdef.h"
25 #include "uidef.h"
26 #include "fft1def.h"
27 #include "fft2def.h"
28 #include "screendef.h"
29 #include "seldef.h"
30
31 void initial_remove_spur(void);
32 int store_new_spur(int pnt);
33 void swap_spurs(int ia, int ib);
34
35 #define SPUR_PERCENTAGE_LIMIT 0.08
36
37 float spur_search_threshold;
38
spursearch_spectrum_cleanup(void)39 void spursearch_spectrum_cleanup(void)
40 {
41 float amp, pos, maxpow, maxamp, refamp, noise;
42 float *spk;
43 int i, j, k, nn;
44 int ia,ib;
45 float t1, t2, t3, t4;
46 float r1;
47 // Collect the smallest out of each group of 32 points in sp_tmp.
48 k=0;
49 for(i=spur_search_first_point; i<spur_search_last_point; i+=32)
50 {
51 sp_tmp[k]=BIG;
52 for(j=0; j<32; j++)
53 {
54 if(spursearch_spectrum[i+j]<sp_tmp[k])sp_tmp[k]=spursearch_spectrum[i+j];
55 }
56 k++;
57 }
58 if(k == 0)return;
59 // Get the average of the smallest points.
60 t1=0;
61 for(i=0; i<k; i++)
62 {
63 t1+=sp_tmp[i];
64 }
65 t1/=k;
66 // Get the average of those smallest points that are below the average of all
67 // the smallest points.
68 noise=0;
69 j=0;
70 for(i=0; i<k; i++)
71 {
72 if(sp_tmp[i] < t1)
73 {
74 noise+=sp_tmp[i];
75 j++;
76 }
77 }
78 noise/=j;
79 // We now have a reasonable estimate of the minimum value of the
80 // noise floor in 32 points in noise.
81 // The RMS value would be somewhere around 7dB/sqrt(3*spur_speknum) higher.
82 // Subtract the noise floor.
83 noise*=pow(10.,0.7/sqrt((float)(3*spur_speknum)));
84 for(i=spur_search_first_point; i<spur_search_last_point; i++)
85 {
86 spursearch_spectrum[i]-=noise;
87 if(spursearch_spectrum[i] < 0)spursearch_spectrum[i]=0;
88 }
89 // anything that is 15dB/sqrt(3*spur_speknum) above the noise floor
90 // has to be a signal.
91 // Locate the max point.
92 spur_search_threshold=noise*pow(10.,1.5/sqrt((float)(3*spur_speknum)));
93 ia=spur_search_first_point;
94 search_signal:;
95 while(ia < spur_search_last_point &&
96 spursearch_spectrum[ia] <spur_search_threshold)ia++;
97 if(ia == spur_search_last_point )goto search_signal_x;
98 ib=ia+1;
99 while(ib < spur_search_last_point &&
100 spursearch_spectrum[ib] >spur_search_threshold)ib++;
101 if(ib==spur_search_last_point && ib-ia < SPUR_SIZE)goto search_signal_x;
102 maxpow=0;
103 for(i=ia; i<ib; i++)
104 {
105 if(spursearch_spectrum[i]>maxpow)
106 {
107 maxpow=spursearch_spectrum[i];
108 k=i;
109 }
110 }
111 parabolic_fit(&, &pos, spursearch_spectrum[k-1],
112 spursearch_spectrum[k], spursearch_spectrum[k+1]);
113 nn=k-SPUR_SIZE/2+1;
114 if(pos < 0)
115 {
116 pos+=1;
117 nn--;
118 }
119 i=pos*NO_OF_SPUR_SPECTRA;
120 if(i>=NO_OF_SPUR_SPECTRA)i=NO_OF_SPUR_SPECTRA-1;
121 spk=&spur_spectra[i*SPUR_SIZE];
122 refamp=fabs(spk[SPUR_SIZE/2]);
123 if(refamp<fabs(spk[SPUR_SIZE/2-1]))refamp=fabs(spk[SPUR_SIZE/2-1]);
124 maxamp=sqrt(maxpow);
125 t2=maxamp/refamp;
126 // Subtract the reference spectrum from the search data.
127 // Compute the percentage of the power that remains within SPUR_SIZE points.
128 t1=0;
129 t3=0;
130 t4=0;
131 for(i=0; i<SPUR_SIZE; i++)
132 {
133 if(spursearch_spectrum[nn+i]<0)
134 {
135 t3=1;
136 t1=1;
137 goto bad_peak;
138 }
139 t1+=spursearch_spectrum[nn+i];
140 r1=pow(sqrt(spursearch_spectrum[nn+i])-t2*fabs(spk[i]),2.0);
141 t3+=r1;
142 if(t4<r1&&(i<2 ||i>=SPUR_SIZE-2))t4=r1;
143 }
144 // t1 is now the total power of the spur/signal.
145 // t3 is the amount of power that we should expect to remain after subtracting a spur.
146 // t4 is the power of the biggest peak we should expect to remain.
147 // A spur might however drift. That would broaden its spectrum over
148 // the long time (3*spur_speknum) that was used to collect spursearch_spectrum.
149 // This means that we may have to accept a residue that is significantly
150 // above the noise floor.
151 if( (t4/noise > 5 && t1/t4 < 1000) ||
152 (t4/noise > 2 && t1/t4 < 300) || t3/t1 > 0.1)
153 {
154 bad_peak:;
155 // More than SPUR_PERCENTAGE_LIMIT of the signal remains.
156 // This means that the spectrum we see is not only a spur.
157 // It might be one of our desired signals!!
158 // Clear the search spectrum from this peak.
159 while( (spursearch_spectrum[ia] > spursearch_spectrum[ia-1] ||
160 spursearch_spectrum[ia] > spursearch_spectrum[ia-2] ||
161 spursearch_spectrum[ia] > spursearch_spectrum[ia-3]) &&
162 ia>spur_search_first_point)ia--;
163 while( (spursearch_spectrum[ib] > spursearch_spectrum[ib+1] ||
164 spursearch_spectrum[ib] > spursearch_spectrum[ib+2] ||
165 spursearch_spectrum[ib] > spursearch_spectrum[ib+3]) &&
166 ib<spur_search_last_point)ib++;
167 for(i=ia; i<ib; i++)spursearch_spectrum[i]=-0.00000001;
168 }
169 ia=ib;
170 goto search_signal;
171 search_signal_x:;
172 if(autospur_point >= spur_search_last_point-SPUR_WIDTH/2)
173 autospur_point=spur_search_first_point+SPUR_WIDTH/2+1;
174 }
175
init_spur_elimination(void)176 void init_spur_elimination(void)
177 {
178 int i, j, k, pnt, ia, ib, fft1_pnt;
179 int spur_search_pnts;
180 float t1;
181 //mailbox[1]=0;
182 if(no_of_spurs >=genparm[MAX_NO_OF_SPURS])goto skip;
183 pnt=-1;
184 fft1_pnt=-1;
185 if(genparm[AFC_ENABLE] == 1)
186 {
187 // if(autospur_point!=spur_search_first_point+SPUR_WIDTH/2+1)return;
188 if(autospur_point!=spur_search_first_point+SPUR_WIDTH/2+1)lirerr(883461);
189 // Check if we are in the main/waterfall graph area.
190 for(i=0; i<no_of_scro; i++)
191 {
192 if( scro[i].x1 <= mouse_x &&
193 scro[i].x2 >= mouse_x &&
194 scro[i].y1 <= mouse_y &&
195 scro[i].y2 >= mouse_y)
196 {
197 if(scro[i].no == WIDE_GRAPH)
198 {
199 pnt=mouse_x-wg_first_xpixel;
200 if(wg.xpoints_per_pixel != 1 && wg.pixels_per_xpoint != 1)
201 {
202 if(wg.xpoints_per_pixel == 0)
203 {
204 pnt=(pnt+(wg.pixels_per_xpoint>>1))/wg.pixels_per_xpoint;
205 }
206 else
207 {
208 pnt*=wg.xpoints_per_pixel;
209 }
210 }
211 pnt+=wg.first_xpoint;
212 fft1_pnt=1;
213 if(genparm[SECOND_FFT_ENABLE] != 0)
214 {
215 pnt*=fft2_to_fft1_ratio;
216 }
217 }
218 if(scro[i].no == HIRES_GRAPH)
219 {
220 pnt=mouse_x-hg_first_xpixel+hg_first_point;
221 fft1_pnt=-1;
222 }
223 }
224 }
225 if(pnt<0 || pnt >= fftx_size)return;
226 // Check that we can locate a spur and move pnt to the peak.
227 k=SPUR_WIDTH/2;
228 if(genparm[SECOND_FFT_ENABLE] != 0) k*=fft2_size/fft1_size;
229 if(fft1_pnt > 0)
230 {
231 if(wg.xpoints_per_pixel > 1)
232 {
233 k=SPUR_WIDTH*wg.xpoints_per_pixel/2;
234 }
235 else
236 {
237 k=0;
238 }
239 }
240 if(k < 1+SPUR_WIDTH/2)k=1+SPUR_WIDTH/2;
241 spur_search_pnts=2*k;
242 if(spur_search_pnts>=fft1_size/2)
243 {
244 spur_search_pnts=fft1_size/2-1;
245 k=spur_search_pnts/2;
246 }
247 if(ffts_nm < spur_speknum-1)return;
248 if(k<9)k=9;
249 ia=pnt-k;
250 ib=pnt+k;
251 if(ia < spur_search_first_point)ia=spur_search_first_point;
252 if(ib > spur_search_last_point)ib=spur_search_last_point;
253 pnt=-1;
254 t1=0;
255 for(i=ia; i<=ib; i++)
256 {
257 if(t1<spursearch_spectrum[i])
258 {
259 t1=spursearch_spectrum[i];
260 pnt=i;
261 }
262 }
263 if(pnt<0)return;
264 // Always point to the first of SPUR_WIDTH points surrounding the peak.
265 pnt-=SPUR_WIDTH/2;
266 if(pnt < spur_search_first_point)return;
267 }
268 else
269 {
270 while(autospur_point < spur_search_last_point &&
271 spursearch_spectrum[autospur_point] <spur_search_threshold)
272 autospur_point++;
273 while(autospur_point < spur_search_last_point &&
274 (spursearch_spectrum[autospur_point] < spursearch_spectrum[autospur_point+1] ||
275 spursearch_spectrum[autospur_point] < spursearch_spectrum[autospur_point+1]) )autospur_point++;
276 if(autospur_point >= spur_search_last_point-SPUR_WIDTH/2)
277 {
278 skip:;
279 autospur_point=spur_search_last_point;
280 return;
281 }
282 pnt=autospur_point-SPUR_WIDTH/2;
283 for(i=0;i<no_of_spurs;i++)
284 {
285 if(abs(pnt-spur_location[i]) <=SPUR_WIDTH/2)goto skip;
286 }
287 autospur_point += SPUR_WIDTH;
288 }
289 //if(pnt > 9390-SPUR_WIDTH/2 && pnt < 9390+SPUR_WIDTH/2)mailbox[1]=1; else mailbox[1]=0;
290 //mailbox[1]=1;
291 spurno=no_of_spurs;
292 spur_ampl[spurno]=1;
293 spur_noise[spurno]=0.001;
294 spur_avgd2[spurno]=0;
295 if(store_new_spur(pnt) != 0)
296 {
297 // if(mailbox[1]==1)fprintf( dmp,"\nSTORE %d failed",pnt);
298 // mailbox[1]=0;
299 return;
300 }
301 //if(mailbox[1]==1)fprintf( dmp,"\nSTORE %d OK",pnt);
302 // Set up phase locked loop parameters.
303 if(spur_phase_lock(ffts_na) != 0)
304 {
305 //if(mailbox[1]==1)fprintf( dmp,"\nPLL failed");
306 //mailbox[1]=0;
307 return;
308 }
309 //if(mailbox[1] == 1)fprintf( dmp,"\nspurno %d loc %d pnt %d",spurno, spur_location[spurno],pnt);
310 no_of_spurs++;
311 // Remove the spur for the latest spur_speknum fftx transforms.
312 initial_remove_spur();
313 // Rearrange spurs so we keep them in order by increasing frequency.
314 for(i=0; i<spurno; i++)
315 {
316 retry_j:;
317 for(j=i+1; j<no_of_spurs; j++)
318 {
319 if(abs(spur_location[i]-spur_location[j]) < 4)
320 {
321 // These two spurs are too close. remove the weakest one.
322 if(spur_ampl[i] > spur_ampl[j])
323 {
324 k=j;
325 }
326 else
327 {
328 k=i;
329 }
330 no_of_spurs--;
331 if(k != no_of_spurs)
332 {
333 remove_spur(k);
334 }
335 goto retry_j;
336 }
337 if(spur_freq[i] > spur_freq[j])
338 {
339 swap_spurs(i,j);
340 }
341 }
342 }
343 //mailbox[1]=0;
344 }
345
346
initial_remove_spur(void)347 void initial_remove_spur(void)
348 {
349 int i, j, ind;
350 int izz, pnt;
351 int ni;
352 int *uind;
353 float t1,t2,t3,t4,r1,r2;
354 float a1, a2, b1, b2, d1, d2;
355 float average_rot,rot;
356 float freq;
357 float phase_slope,phase_curv;
358 float *z;
359 short int *zmmx;
360 // PLL was sucessful. We now have a "local oscillator"
361 // that is locked to the spur. The LO is characterised
362 // by phase and amplitude parameters.
363 // Subtract the LO from the original data in fftx.
364 uind=&spur_ind[spurno*max_fftxn];
365 average_rot=spur_freq[spurno]*spur_freq_factor;
366 a1=spur_ampl[spurno]*cos(spur_d0pha[spurno]);
367 a2=spur_ampl[spurno]*sin(spur_d0pha[spurno]);
368 b1=cos(spur_d1pha[spurno]);
369 b2=sin(spur_d1pha[spurno]);
370 d1=cos(spur_d2pha[spurno]);
371 d2=sin(spur_d2pha[spurno]);
372 phase_slope=spur_d1pha[spurno];
373 phase_curv=spur_d2pha[spurno];
374 ni=(ffts_na+fftxn_mask)&fftxn_mask;
375 pnt=spur_location[spurno];
376 izz=spur_speknum-1;
377 while(izz >= 0)
378 {
379 rot=-0.5*phase_slope/PI_L;
380 i=average_rot-rot+0.5;
381 rot+=i;
382 freq=rot/spur_freq_factor;
383 j=(int)(freq)+2-pnt-SPUR_SIZE/2;
384 if(j<0)j=0;
385 j=1-j;
386 if(j<0)j=0;
387 ind=uind[ni];
388 if(sw_onechan)
389 {
390 if((j^(spur_location[spurno]&1))==1)
391 {
392 t1=-a1;
393 t2=-a2;
394 }
395 else
396 {
397 t1=a1;
398 t2=a2;
399 }
400 if(swmmx_fft2)
401 {
402 zmmx=&fft2_short_int[twice_rxchan*(ni*fft2_size+pnt)];
403 for(i=0; i<SPUR_WIDTH; i++)
404 {
405 zmmx[2*i ]-=spur_spectra[ind+i]*t1;
406 zmmx[2*i+1]-=spur_spectra[ind+i]*t2;
407 }
408 }
409 else
410 {
411 z=&fftx[twice_rxchan*(ni*fftx_size+pnt)];
412 for(i=0; i<SPUR_WIDTH; i++)
413 {
414 z[2*i ]-=spur_spectra[ind+i]*t1;
415 z[2*i+1]-=spur_spectra[ind+i]*t2;
416 }
417 }
418 }
419 else
420 {
421 if((j^(spur_location[spurno]&1))==1)
422 {
423 t1=-sp_c1*a1;
424 t2=-sp_c1*a2;
425 t3=-sp_c2*a1+sp_c3*a2;
426 t4=-sp_c2*a2-sp_c3*a1;
427 }
428 else
429 {
430 t1=sp_c1*a1;
431 t2=sp_c1*a2;
432 t3=sp_c2*a1-sp_c3*a2;
433 t4=sp_c2*a2+sp_c3*a1;
434 }
435 if(swmmx_fft2)
436 {
437 zmmx=&fft2_short_int[twice_rxchan*(ni*fft2_size+pnt)];
438 for(i=0; i<SPUR_WIDTH; i++)
439 {
440 zmmx[4*i ]-=spur_spectra[ind+i]*t1;
441 zmmx[4*i+1]-=spur_spectra[ind+i]*t2;
442 zmmx[4*i+2]-=spur_spectra[ind+i]*t3;
443 zmmx[4*i+3]-=spur_spectra[ind+i]*t4;
444 }
445 }
446 else
447 {
448 z=&fftx[twice_rxchan*(ni*fftx_size+pnt)];
449 for(i=0; i<SPUR_WIDTH; i++)
450 {
451 z[4*i ]-=spur_spectra[ind+i]*t1;
452 z[4*i+1]-=spur_spectra[ind+i]*t2;
453 z[4*i+2]-=spur_spectra[ind+i]*t3;
454 z[4*i+3]-=spur_spectra[ind+i]*t4;
455 }
456 }
457 }
458 ni=(ni+fftxn_mask)&fftxn_mask;
459 izz--;
460 r1=a1*b1+a2*b2;
461 a2=a2*b1-a1*b2;
462 a1=r1;
463 r2=b1*d1+b2*d2;
464 b2=b2*d1-b1*d2;
465 b1=r2;
466 phase_slope-=phase_curv;
467 }
468 }
469
470
spurspek_norm(void)471 void spurspek_norm(void)
472 {
473 int i;
474 float t1,t2;
475 // We have a power spectrum in spur_power[]
476 // Assume that the level at the spectrum end points correspond
477 // to the white noise floor.
478 // Remove the average of the end points and normalise for the sum
479 // of all points = 1
480 t1=0.5*(spur_power[0]+spur_power[SPUR_WIDTH-1]);
481 t2=0;
482 for(i=0; i<SPUR_WIDTH; i++)
483 {
484 spur_power[i]-=t1;
485 if(spur_power[i] <0)spur_power[i]=0;
486 t2+=spur_power[i];
487 }
488 for(i=0; i<SPUR_WIDTH; i++)
489 {
490 spur_power[i]/=t2;
491 }
492 }
493
make_spur_pol(void)494 void make_spur_pol(void)
495 {
496 int i;
497 float t1,t2;
498 float x2,y2,re_xy,im_xy;
499 float noi2, x2s, y2s, sina;
500 x2=0;
501 y2=0;
502 re_xy=0;
503 im_xy=0;
504 for(i=0; i<SPUR_WIDTH; i++)
505 {
506 x2+=spur_power[i]*spur_pxy[i].x2;
507 y2+=spur_power[i]*spur_pxy[i].y2;
508 re_xy+=spur_power[i]*spur_pxy[i].re_xy;
509 im_xy+=spur_power[i]*spur_pxy[i].im_xy;
510 }
511 t1=x2+y2;
512 x2/=t1;
513 y2/=t1;
514 re_xy/=t1;
515 im_xy/=t1;
516 // Extract polarization from the power weighted average of
517 // powers and cross channel correlation.
518 // Now we have x2,y2 (real values) and xy (complex).
519 // For explanation purposes, assume im_xy == 0, which corresponds to linear
520 // polarization. The signal vill then be polarised in a plane.
521 // a = angle between polarization plane and the horisontal antenna.
522 // Assume that the noise level n is the same in the two antennas, and that
523 // the noise is uncorrelated.
524 // We then find:
525 // x2 = cos(a)**2 + n**2
526 // y2 = sin(a)**2 + n**2
527 // xy = sin(a)*cos(a)
528 // From this we find: x2 * y2 - xy*xy = n**2 + n**4
529 // Neglect n**4:
530 // cos(a)=sqr( x2 - ( x2 * y2 - xy*xy) )
531 // sin(a)=sqr( y2 - ( x2 * y2 - xy*xy) )
532 // The transformation formula to use for rotating the polarization
533 // plane to produce new signals A and B, where A has all the signal and B
534 // only noise, will then be:
535 // A = X * cos(a) + Y * sin(a)
536 // B = Y * cos(a) - X * sin(a)
537 // Extending to im_xy != 0 the transformation becomes
538 // re_A=C1*re_X+C2*re_Y-C3*im_Y
539 // im_A=C1*im_X+C2*im_Y+C3*re_Y
540 // re_B=C1*re_Y-C2*re_X-C3*im_X
541 // im_B=C1*im_Y-C2*im_X+C3*re_X
542 // C1 = cos(a)
543 // C2 = sin(a) * re_xy / sqr( re_xy**2 + im_xy**2)
544 // C3 = sin(a) * im_xy / sqr( re_xy**2 + im_xy**2)
545 t2=re_xy*re_xy+im_xy*im_xy;
546 noi2=x2*y2-t2;
547 x2s=x2-noi2;
548 //if(mailbox[1]==1)fprintf( dmp,"\nXY x2 %f y2 %f ",x2,y2);
549 if(x2s > 0)
550 {
551 y2s=y2-noi2;
552 sp_c1=sqrt(x2s);
553 if(y2s > 0 && t2 > 0)
554 {
555 sina=sqrt(y2s);
556 sp_c2=sina*re_xy/sqrt(t2);
557 sp_c3=-sina*im_xy/sqrt(t2);
558 t1=sqrt(sp_c1*sp_c1+sp_c2*sp_c2+sp_c3*sp_c3);
559 if(sp_c2 < 0)t1=-t1;
560 sp_c1/=t1;
561 sp_c2/=t1;
562 sp_c3/=t1;
563 }
564 else
565 {
566 if(x2 > y2)
567 {
568 sp_c1=1;
569 sp_c2=0;
570 }
571 else
572 {
573 sp_c1=0;
574 sp_c2=1;
575 }
576 sp_c3=0;
577 }
578 }
579 else
580 {
581 sp_c1=1;
582 sp_c2=0;
583 sp_c3=0;
584 }
585 i=3*spurno;
586 spur_pol[i ]=sp_c1;
587 spur_pol[i+1]=sp_c2;
588 spur_pol[i+2]=sp_c3;
589 //if(mailbox[1]==1)fprintf( dmp,"\nPOL: %f %f %f",sp_c1,sp_c2,sp_c3);
590 }
591
make_spur_freq(void)592 int make_spur_freq(void)
593 {
594 int i,k;
595 float amp,pos;
596 float t1,t2,t3;
597 // Get the spur frequency with decimals
598 t2=0;
599 k=0;
600 for(i=0; i<SPUR_WIDTH; i++)
601 {
602 if(t2<spur_power[i])
603 {
604 t2=spur_power[i];
605 k=i;
606 }
607 }
608 if(k == 0 || k == SPUR_WIDTH-1)return -1;
609 t1=sqrt(spur_power[k-1]);
610 t3=sqrt(spur_power[k+1]);
611 t2=sqrt(t2);
612 parabolic_fit( &, &pos, t1, t2, t3);
613 spur_freq[spurno]=spur_location[spurno]+k+pos;
614 return 0;
615 }
616
617
store_new_spur(int pnt)618 int store_new_spur(int pnt)
619 {
620 int i, k, np;
621 float t1;
622 short int *zmmx;
623 float *z, *pwra;
624 TWOCHAN_POWER *pxy;
625 spur_flag[spurno]=0;
626 // Store old transforms in spur table for the new spur and
627 // accumulate power spectrum.
628 spurp0=spurno*spur_block;
629 np=(ffts_na-spur_speknum+max_fftxn)&fftxn_mask;
630 if(swmmx_fft2)
631 {
632 if(sw_onechan)
633 {
634 // With one channel only, just use the power spectrum.
635 for(i=0; i<SPUR_WIDTH; i++)spur_power[i]=0;
636 while(np != ffts_na)
637 {
638 k=spurp0+np*SPUR_WIDTH*twice_rxchan;
639 zmmx=&fft2_short_int[twice_rxchan*(np*fft2_size+pnt)];
640 pwra=&fftx_pwr[np*fftx_size+pnt];
641 for(i=0; i<SPUR_WIDTH; i++)
642 {
643 spur_table_mmx[2*i+k ]=zmmx[2*i ];
644 spur_table_mmx[2*i+k+1]=zmmx[2*i+1];
645 spur_power[i]+=pwra[i];
646 }
647 np=(np+1)&fftxn_mask;
648 }
649 }
650 else //ui.rx_channels = 2
651 {
652 // We have two channels and polarization is unknown.
653 // First make an average of channel powers and correlations.
654 for(i=0; i<SPUR_WIDTH; i++)
655 {
656 spur_pxy[i].x2=0;
657 spur_pxy[i].y2=0;
658 spur_pxy[i].re_xy=0;
659 spur_pxy[i].im_xy=0;
660 }
661 while(np != ffts_na)
662 {
663 pxy=(TWOCHAN_POWER*)(&fftx_xypower[np*fftx_size+pnt].x2);
664 k=spurp0+np*SPUR_WIDTH*twice_rxchan;
665 zmmx=&fft2_short_int[twice_rxchan*(np*fft2_size+pnt)];
666 for(i=0; i<SPUR_WIDTH; i++)
667 {
668 spur_table_mmx[4*i+k ]=zmmx[4*i ];
669 spur_table_mmx[4*i+k+1]=zmmx[4*i+1];
670 spur_table_mmx[4*i+k+2]=zmmx[4*i+2];
671 spur_table_mmx[4*i+k+3]=zmmx[4*i+3];
672 spur_pxy[i].x2+=pxy[i].x2;
673 spur_pxy[i].y2+=pxy[i].y2;
674 spur_pxy[i].re_xy+=pxy[i].re_xy;
675 spur_pxy[i].im_xy+=pxy[i].im_xy;
676 }
677 np=(np+1)&fftxn_mask;
678 }
679 }
680 }
681 else
682 {
683 if(sw_onechan)
684 {
685 // With one channel only, just use the power spectrum.
686 for(i=0; i<SPUR_WIDTH; i++)spur_power[i]=0;
687 while(np != ffts_na)
688 {
689 k=spurp0+np*SPUR_WIDTH*twice_rxchan;
690 z=&fftx[twice_rxchan*(np*fftx_size+pnt)];
691 pwra=&fftx_pwr[np*fftx_size+pnt];
692 for(i=0; i<SPUR_WIDTH; i++)
693 {
694 spur_table[2*i+k ]=z[2*i ];
695 spur_table[2*i+k+1]=z[2*i+1];
696 spur_power[i]+=pwra[i];
697 }
698 np=(np+1)&fftxn_mask;
699 }
700 }
701 else //ui.rx_channels = 2
702 {
703 // We have two channels and polarization is unknown.
704 // First make an average of channel powers and correlations.
705 for(i=0; i<SPUR_WIDTH; i++)
706 {
707 spur_pxy[i].x2=0;
708 spur_pxy[i].y2=0;
709 spur_pxy[i].re_xy=0;
710 spur_pxy[i].im_xy=0;
711 }
712 while(np != ffts_na)
713 {
714 pxy=(TWOCHAN_POWER*)(&fftx_xypower[np*fftx_size+pnt].x2);
715 k=spurp0+np*SPUR_WIDTH*twice_rxchan;
716 z=&fftx[twice_rxchan*(np*fftx_size+pnt)];
717 for(i=0; i<SPUR_WIDTH; i++)
718 {
719 spur_table[4*i+k ]=z[4*i ];
720 spur_table[4*i+k+1]=z[4*i+1];
721 spur_table[4*i+k+2]=z[4*i+2];
722 spur_table[4*i+k+3]=z[4*i+3];
723 spur_pxy[i].x2+=pxy[i].x2;
724 spur_pxy[i].y2+=pxy[i].y2;
725 spur_pxy[i].re_xy+=pxy[i].re_xy;
726 spur_pxy[i].im_xy+=pxy[i].im_xy;
727 }
728 np=(np+1)&fftxn_mask;
729 }
730 }
731 }
732 if(sw_onechan)
733 {
734 spurspek_norm();
735 }
736 else
737 {
738 // Two channels. Extract the power for a signal using the channel
739 // correlation (Slightly better than just adding x2 and y2)
740 for(i=0; i<SPUR_WIDTH; i++)
741 {
742 t1=spur_pxy[i].x2+spur_pxy[i].y2;
743 spur_power[i]=t1+2*(spur_pxy[i].re_xy*spur_pxy[i].re_xy+
744 spur_pxy[i].im_xy*spur_pxy[i].im_xy
745 -spur_pxy[i].x2*spur_pxy[i].y2 ) / t1;
746 }
747 spurspek_norm();
748 make_spur_pol();
749 }
750 spur_location[spurno]=pnt;
751 return make_spur_freq();
752 }
753
swap_spurs(int ia,int ib)754 void swap_spurs(int ia, int ib)
755 {
756 int i, k;
757 short int m;
758 float t1;
759 k=spur_flag[ia];
760 spur_flag[ia]=spur_flag[ib];
761 spur_flag[ib]=k;
762 k=spur_location[ia];
763 spur_location[ia]=spur_location[ib];
764 spur_location[ib]=k;
765 t1=spur_d0pha[ia];
766 spur_d0pha[ia]=spur_d0pha[ib];
767 spur_d0pha[ib]=t1;
768 t1=spur_d1pha[ia];
769 spur_d1pha[ia]=spur_d1pha[ib];
770 spur_d1pha[ib]=t1;
771 t1=spur_d2pha[ia];
772 spur_d2pha[ia]=spur_d2pha[ib];
773 spur_d2pha[ib]=t1;
774 t1=spur_ampl[ia];
775 spur_ampl[ia]=spur_ampl[ib];
776 spur_ampl[ib]=t1;
777 t1=spur_noise[ia];
778 spur_noise[ia]=spur_noise[ib];
779 spur_noise[ib]=t1;
780 t1=spur_avgd2[ia];
781 spur_avgd2[ia]=spur_avgd2[ib];
782 spur_avgd2[ib]=t1;
783 t1=spur_pol[ia];
784 spur_pol[ia]=spur_pol[ib];
785 spur_pol[ib]=t1;
786 t1=spur_freq[ia];
787 spur_freq[ia]=spur_freq[ib];
788 spur_freq[ib]=t1;
789 if(swmmx_fft2)
790 {
791 for(i=0; i<spur_block; i++)
792 {
793 m=spur_table_mmx[ia*spur_block+i];
794 spur_table_mmx[ia*spur_block+i]=
795 spur_table_mmx[ib*spur_block+i];
796 spur_table_mmx[ib*spur_block+i]=m;
797 }
798 }
799 else
800 {
801 for(i=0; i<spur_block; i++)
802 {
803 t1=spur_table[ia*spur_block+i];
804 spur_table[ia*spur_block+i]=spur_table[ib*spur_block+i];
805 spur_table[ib*spur_block+i]=t1;
806 }
807 }
808 for(i=0; i<twice_rxchan*max_fftxn; i++)
809 {
810 t1=spur_signal[twice_rxchan*max_fftxn*ia+i];
811 spur_signal[twice_rxchan*max_fftxn*ia+i]=
812 spur_signal[twice_rxchan*max_fftxn*ib+i];
813 spur_signal[twice_rxchan*max_fftxn*ib+i]=t1;
814 }
815 for(i=0; i<max_fftxn; i++)
816 {
817 k=spur_ind[max_fftxn*ia+i];
818 spur_ind[max_fftxn*ia+i]=spur_ind[max_fftxn*ib+i];
819 spur_ind[max_fftxn*ib+i]=k;
820 }
821 }
822
init_spur_spectra(void)823 void init_spur_spectra(void)
824 {
825 int i, j, ifreq;
826 int iwin;
827 float *tmp_window,*spk, *tmp_buf;
828 float t1, t2, t3;
829 #define SPI_N 8
830 #define SPI_SIZE (1<<SPI_N)
831 unsigned short int *tmp_permute;
832 COSIN_TABLE *tmp_tab;
833 tmp_permute=malloc(SPI_SIZE*sizeof(short int));
834 if(tmp_permute == NULL)
835 {
836 alloc_err:;
837 lirerr(1103);
838 return;
839 }
840 tmp_tab=malloc(SPI_SIZE*sizeof(COSIN_TABLE));
841 if(tmp_tab == NULL)
842 {
843 alloc_err1:;
844 free(tmp_permute);
845 goto alloc_err;
846 }
847 tmp_window=malloc(SPI_SIZE*sizeof(float));
848 if(tmp_window == NULL)
849 {
850 alloc_err2:;
851 free(tmp_tab);
852 goto alloc_err1;
853 }
854 tmp_buf=malloc(2*SPI_SIZE*sizeof(float));
855 if(tmp_buf == NULL)
856 {
857 goto alloc_err2;
858 }
859 if(genparm[SECOND_FFT_ENABLE] != 0)
860 {
861 iwin=genparm[SECOND_FFT_SINPOW];
862 }
863 else
864 {
865 iwin=genparm[FIRST_FFT_SINPOW];
866 }
867 init_fft(0,SPI_N, SPI_SIZE, tmp_tab, tmp_permute);
868 if(iwin !=0)
869 {
870 make_window(4,SPI_SIZE,iwin, tmp_window);
871 }
872 else
873 {
874 for(i=0; i<SPI_SIZE; i++)
875 {
876 tmp_window[i]=1;
877 }
878 }
879 // Generate a sine wave with frequency SPUR_SIZE/4 to SPUR_SIZE/4+1
880 // in NO_OF_SPUR_SPECTRA steps.
881 for(ifreq=0; ifreq<NO_OF_SPUR_SPECTRA; ifreq++)
882 {
883 t1=1+SPUR_SIZE/4+(float)(ifreq)/NO_OF_SPUR_SPECTRA;
884 t1=2*PI_L*t1/SPI_SIZE;
885 t2=0;
886 for(i=0; i<SPI_SIZE; i++)
887 {
888 tmp_buf[2*i ]=tmp_window[i]*sin(t2);
889 tmp_buf[2*i+1]=tmp_window[i]*cos(t2);
890 t2+=t1;
891 }
892 fftforward(SPI_SIZE, SPI_N, tmp_buf, tmp_tab, tmp_permute,FALSE);
893 // For simplicity we determine phases from the weighted sum of
894 // the individual spectral lines.
895 // Calculate the phase of the spectra we have got here with the
896 // same procedure and twist the phase angle so we store these
897 // spectra in phase zero with the simplified phase method.
898 t1=0;
899 t2=0;
900 j=1;
901 for(i=0; i<SPUR_SIZE; i++)
902 {
903 t3=j*(tmp_buf[2*i]*tmp_buf[2*i]+tmp_buf[2*i+1]*tmp_buf[2*i+1]);
904 j*=-1;
905 t1+=t3*tmp_buf[2*i ];
906 t2+=t3*tmp_buf[2*i+1];
907 }
908 t3=sqrt(t1*t1+t2*t2);
909 t1/=t3;
910 t2/=t3;
911 for(i=0; i<SPUR_SIZE; i++)
912 {
913 t3=t1*tmp_buf[2*i]+t2*tmp_buf[2*i+1];
914 tmp_buf[2*i+1]=t1*tmp_buf[2*i+1]-t2*tmp_buf[2*i];
915 tmp_buf[2*i]=t3;
916 }
917 // Store the real part only
918 // The imaginary part we have is just because a too small size
919 // was used for SPI_SIZE to make program start fast
920 // (The real size, fftx_n may go very slow).
921 spk=&spur_spectra[ifreq*SPUR_SIZE];
922 for(i=0; i<SPUR_SIZE; i++)spk[i]=tmp_buf[2*i];
923 // Normalise the spectrum for sum of squares equal to 1
924 t1=0;
925 for(i=0; i<SPUR_SIZE; i++)
926 {
927 t1+=spk[i]*spk[i];
928 }
929 t1=sqrt(t1);
930 for(i=0; i<SPUR_SIZE; i++)
931 {
932 spk[i]/=t1;
933 }
934 }
935 free(tmp_permute);
936 free(tmp_tab);
937 free(tmp_window);
938 for(i=0; i<fftx_size; i++)spursearch_spectrum[i]=0;
939 }
940
average_curvature(float * z,int siz)941 float average_curvature(float *z, int siz)
942 {
943 int i, k, m0, m1, m2;
944 float t1,t2,t3;
945 // Calculate the average at 3 well separated points.
946 // Three values are needed for the second derivative.
947 k=(siz-1)/3;
948 if(k<1)k=1;
949 m1=(siz-k-1)/2;
950 m0=m1-k;
951 if(m0<0)
952 {
953 m0=0;
954 m1=m0+k;
955 }
956 m2=m1+k;
957 if(m2+k < siz-1)
958 {
959 m0++;
960 m1++;
961 m2++;
962 }
963 t1=0;
964 t2=0;
965 t3=0;
966 for(i=0; i<k; i++)
967 {
968 t1+=z[m0+i];
969 t2+=z[m1+i];
970 t3+=z[m2+i];
971 }
972 return (t3-2*t2+t1)/(k*k*k);
973 }
974
average_slope(float * z,int siz)975 float average_slope(float *z, int siz)
976 {
977 int i, k, m1, m2;
978 float t2,t3;
979 // Calculate the average at 2 well separated points.
980 // Two values are needed for the first derivative.
981 k=siz/2;
982 m1=0;
983 m2=k;
984 t2=0;
985 t3=0;
986 for(i=0; i<k; i++)
987 {
988 t2+=z[m1+i];
989 t3+=z[m2+i];
990 }
991 return (t3-t2)/(k*k);
992 }
993
994
995
remove_phasejumps(float * z,int siz)996 void remove_phasejumps(float *z, int siz)
997 {
998 float t1;
999 int i;
1000 // Remove any 2*PI discontinuities in a phase function.
1001 t1=0;
1002 for(i=1; i<siz; i++)
1003 {
1004 z[i]+=t1;
1005 if(z[i]-z[i-1] > PI_L)
1006 {
1007 z[i]-=2*PI_L;
1008 t1-=2*PI_L;
1009 }
1010 if(z[i]-z[i-1] < -PI_L)
1011 {
1012 z[i]+=2*PI_L;
1013 t1+=2*PI_L;
1014 }
1015 }
1016 }
1017
1018
complex_lowpass(float * zin,float * zout,int nn,int siz)1019 void complex_lowpass(float *zin, float *zout, int nn, int siz)
1020 {
1021 int i,j,k,avgnum;
1022 float r1,r2,t1,t2,t3;
1023 // Average a complex function over avgnum points
1024 avgnum=nn|1;
1025 if(avgnum > nn && avgnum > siz/4)avgnum-=2;
1026 if(avgnum < 1 )
1027 {
1028 lirerr(1105);
1029 return;
1030 }
1031 t1=0;
1032 t2=0;
1033 t3=1.0/avgnum;
1034 for(i=0; i<avgnum; i++)
1035 {
1036 t1+=zin[2*i ];
1037 t2+=zin[2*i+1];
1038 }
1039 j=1+avgnum/2;
1040 r1=t1*t3;
1041 r2=t2*t3;
1042 for(i=0; i<j; i++)
1043 {
1044 zout[2*i ]=r1;
1045 zout[2*i+1]=r2;
1046 }
1047 i=0;
1048 k=avgnum;
1049 while(k < siz)
1050 {
1051 t1+=zin[2*k ]-zin[2*i ];
1052 t2+=zin[2*k+1]-zin[2*i+1];
1053 zout[2*j ]=t3*t1;
1054 zout[2*j+1]=t3*t2;
1055 i++;
1056 j++;
1057 k++;
1058 }
1059 t1*=t3;
1060 t2*=t3;
1061 while(j < siz)
1062 {
1063 zout[2*j ]=t1;
1064 zout[2*j+1]=t2;
1065 j++;
1066 }
1067 }
1068
shift_spur_table(int j)1069 void shift_spur_table(int j)
1070 {
1071 int i, knd, ni, nj, shift;
1072 // The spur is no longer centered in the SPUR_WIDTH points
1073 // that we are saving.
1074 nj=(ffts_na+1)&fftxn_mask;
1075 if(j < 0)shift=-1; else shift=1;
1076 spur_location[spurno]+=shift;
1077 if(spur_location[spurno] < SPUR_WIDTH)
1078 {
1079 spur_flag[spurno]=1;
1080 spur_location[spurno]=2*SPUR_WIDTH;
1081 return;
1082 }
1083 if(spur_location[spurno] > fftx_size-SPUR_WIDTH)
1084 {
1085 spur_flag[spurno]=1;
1086 spur_location[spurno]=fftx_size-2*SPUR_WIDTH;
1087 return;
1088 }
1089 ni=(ffts_na-spur_speknum+max_fftxn)&fftxn_mask;
1090 if(sw_onechan)
1091 {
1092 if(swmmx_fft2)
1093 {
1094 if(shift == 1)
1095 {
1096 while(ni != nj)
1097 {
1098 knd=spurp0+ni*SPUR_WIDTH*twice_rxchan;
1099 for(i=1; i<SPUR_WIDTH; i++)
1100 {
1101 spur_table_mmx[2*i+knd-2]=spur_table_mmx[2*i+knd ];
1102 spur_table_mmx[2*i+knd-1]=spur_table_mmx[2*i+knd+1];
1103 }
1104 spur_table_mmx[knd+2*SPUR_WIDTH-2]=0;
1105 spur_table_mmx[knd+2*SPUR_WIDTH-1]=0;
1106 ni=(ni+1)&fftxn_mask;
1107 }
1108 }
1109 else
1110 {
1111 while(ni != nj)
1112 {
1113 knd=spurp0+ni*SPUR_WIDTH*twice_rxchan;
1114 for(i=SPUR_WIDTH-1; i>0; i--)
1115 {
1116 spur_table_mmx[2*i+knd ]=spur_table_mmx[2*i+knd-2];
1117 spur_table_mmx[2*i+knd+1]=spur_table_mmx[2*i+knd-1];
1118 }
1119 spur_table_mmx[knd ]=0;
1120 spur_table_mmx[knd+1]=0;
1121 ni=(ni+1)&fftxn_mask;
1122 }
1123 }
1124 }
1125 else
1126 {
1127 if(shift == 1)
1128 {
1129 while(ni != nj)
1130 {
1131 knd=spurp0+ni*SPUR_WIDTH*twice_rxchan;
1132 for(i=1; i<SPUR_WIDTH; i++)
1133 {
1134 spur_table[2*i+knd-2]=spur_table[2*i+knd ];
1135 spur_table[2*i+knd-1]=spur_table[2*i+knd+1];
1136 }
1137 spur_table[knd+2*SPUR_WIDTH-2]=0;
1138 spur_table[knd+2*SPUR_WIDTH-1]=0;
1139 ni=(ni+1)&fftxn_mask;
1140 }
1141 }
1142 else
1143 {
1144 while(ni != nj)
1145 {
1146 knd=spurp0+ni*SPUR_WIDTH*twice_rxchan;
1147 for(i=SPUR_WIDTH-1; i>0; i--)
1148 {
1149 spur_table[2*i+knd ]=spur_table[2*i+knd-2];
1150 spur_table[2*i+knd+1]=spur_table[2*i+knd-1];
1151 }
1152 spur_table[knd ]=0;
1153 spur_table[knd+1]=0;
1154 ni=(ni+1)&fftxn_mask;
1155 }
1156 }
1157 }
1158 }
1159 else
1160 {
1161 if(swmmx_fft2)
1162 {
1163 if(shift == 1)
1164 {
1165 while(ni != nj)
1166 {
1167 knd=spurp0+ni*SPUR_WIDTH*twice_rxchan;
1168 for(i=1; i<SPUR_WIDTH; i++)
1169 {
1170 spur_table_mmx[4*i+knd-4]=spur_table_mmx[4*i+knd ];
1171 spur_table_mmx[4*i+knd-3]=spur_table_mmx[4*i+knd+1];
1172 spur_table_mmx[4*i+knd-2]=spur_table_mmx[4*i+knd+2];
1173 spur_table_mmx[4*i+knd-1]=spur_table_mmx[4*i+knd+3];
1174 }
1175 spur_table_mmx[knd+4*SPUR_WIDTH-4]=0;
1176 spur_table_mmx[knd+4*SPUR_WIDTH-3]=0;
1177 spur_table_mmx[knd+4*SPUR_WIDTH-2]=0;
1178 spur_table_mmx[knd+4*SPUR_WIDTH-1]=0;
1179 ni=(ni+1)&fftxn_mask;
1180 }
1181 }
1182 else
1183 {
1184 while(ni != nj)
1185 {
1186 knd=spurp0+ni*SPUR_WIDTH*twice_rxchan;
1187 for(i=SPUR_WIDTH-1; i>0; i--)
1188 {
1189 spur_table_mmx[4*i+knd ]=spur_table_mmx[4*i+knd-4];
1190 spur_table_mmx[4*i+knd+1]=spur_table_mmx[4*i+knd-3];
1191 spur_table_mmx[4*i+knd+2]=spur_table_mmx[4*i+knd-2];
1192 spur_table_mmx[4*i+knd+3]=spur_table_mmx[4*i+knd-1];
1193 }
1194 spur_table_mmx[knd ]=0;
1195 spur_table_mmx[knd+1]=0;
1196 spur_table_mmx[knd+2]=0;
1197 spur_table_mmx[knd+3]=0;
1198 ni=(ni+1)&fftxn_mask;
1199 }
1200 }
1201 }
1202 else
1203 {
1204 if(shift == 1)
1205 {
1206 while(ni != nj)
1207 {
1208 knd=spurp0+ni*SPUR_WIDTH*twice_rxchan;
1209 for(i=1; i<SPUR_WIDTH; i++)
1210 {
1211 spur_table[4*i+knd-4]=spur_table[4*i+knd ];
1212 spur_table[4*i+knd-3]=spur_table[4*i+knd+1];
1213 spur_table[4*i+knd-2]=spur_table[4*i+knd+2];
1214 spur_table[4*i+knd-1]=spur_table[4*i+knd+3];
1215 }
1216 spur_table[knd+4*SPUR_WIDTH-4]=0;
1217 spur_table[knd+4*SPUR_WIDTH-3]=0;
1218 spur_table[knd+4*SPUR_WIDTH-2]=0;
1219 spur_table[knd+4*SPUR_WIDTH-1]=0;
1220 ni=(ni+1)&fftxn_mask;
1221 }
1222 }
1223 else
1224 {
1225 while(ni != nj)
1226 {
1227 knd=spurp0+ni*SPUR_WIDTH*twice_rxchan;
1228 for(i=SPUR_WIDTH-1; i>0; i--)
1229 {
1230 spur_table[4*i+knd ]=spur_table[4*i+knd-4];
1231 spur_table[4*i+knd+1]=spur_table[4*i+knd-3];
1232 spur_table[4*i+knd+2]=spur_table[4*i+knd-2];
1233 spur_table[4*i+knd+3]=spur_table[4*i+knd-1];
1234 }
1235 spur_table[knd ]=0;
1236 spur_table[knd+1]=0;
1237 spur_table[knd+2]=0;
1238 spur_table[knd+3]=0;
1239 ni=(ni+1)&fftxn_mask;
1240 }
1241 }
1242 }
1243 }
1244 }
1245
spur_phase_lock(int nx)1246 int spur_phase_lock(int nx)
1247 {
1248 int i,k,izz;
1249 int np, nn;
1250 float t1,t2,t3,t4,r1,r2,r3,r4;
1251 float *zsig;
1252 sp_d0=0;
1253 sp_d1=0;
1254 sp_d2=0;
1255 spur_d0pha[spurno]=0;
1256 spur_d1pha[spurno]=0;
1257 spur_d2pha[spurno]=0;
1258 // There is a spur present and we have it's spectrum.
1259 // In case we have two channels we also have polarization parameters
1260 // Calculate amplitude and phase data from spectrum weighted amplitudes.
1261 zsig=&spur_signal[twice_rxchan*max_fftxn*spurno];
1262 izz=0;
1263 nn=(nx-spur_speknum+max_fftxn)&fftxn_mask;
1264 np=nn;
1265 if(swmmx_fft2)
1266 {
1267 if(sw_onechan)
1268 {
1269 while(np != nx)
1270 {
1271 k=spurp0+np*SPUR_WIDTH*twice_rxchan;
1272 t1=0;
1273 t2=0;
1274 for(i=0; i<SPUR_WIDTH-1; i+=2)
1275 {
1276 t1+=spur_power[i]*spur_table_mmx[2*i+k ]-
1277 spur_power[i+1]*spur_table_mmx[2*i+k+2];
1278 t2+=spur_power[i]*spur_table_mmx[2*i+k+1]-
1279 spur_power[i+1]*spur_table_mmx[2*i+k+3];
1280 }
1281 if((spur_location[spurno]&1)==1)
1282 {
1283 t1=-t1;
1284 t2=-t2;
1285 }
1286 zsig[2*np ]=t1;
1287 zsig[2*np+1]=t2;
1288 sp_sig[2*izz ]=t1;
1289 sp_sig[2*izz+1]=t2;
1290 izz++;
1291 np=(np+1)&fftxn_mask;
1292 }
1293 }
1294 else
1295 {
1296 while(np != nx)
1297 {
1298 k=spurp0+np*SPUR_WIDTH*twice_rxchan;
1299 t1=0;
1300 t2=0;
1301 t3=0;
1302 t4=0;
1303 for(i=0; i<SPUR_WIDTH-1; i+=2)
1304 {
1305 t1+=spur_power[i ]*spur_table_mmx[4*i+k ]-
1306 spur_power[i+1]*spur_table_mmx[4*i+k+4];
1307 t2+=spur_power[i ]*spur_table_mmx[4*i+k+1]-
1308 spur_power[i+1]*spur_table_mmx[4*i+k+5];
1309 t3+=spur_power[i ]*spur_table_mmx[4*i+k+2]-
1310 spur_power[i+1]*spur_table_mmx[4*i+k+6];
1311 t4+=spur_power[i ]*spur_table_mmx[4*i+k+3]-
1312 spur_power[i+1]*spur_table_mmx[4*i+k+7];
1313 }
1314 if((spur_location[spurno]&1)==1)
1315 {
1316 r1=-sp_c1*t1-sp_c2*t3-sp_c3*t4;
1317 r2=-sp_c1*t2-sp_c2*t4+sp_c3*t3;
1318 r3=-sp_c1*t3+sp_c2*t1-sp_c3*t2;
1319 r4=-sp_c1*t4+sp_c2*t2+sp_c3*t1;
1320 }
1321 else
1322 {
1323 r1=sp_c1*t1+sp_c2*t3+sp_c3*t4;
1324 r2=sp_c1*t2+sp_c2*t4-sp_c3*t3;
1325 r3=sp_c1*t3-sp_c2*t1+sp_c3*t2;
1326 r4=sp_c1*t4-sp_c2*t2-sp_c3*t1;
1327 }
1328 zsig[4*np ]=r1;
1329 zsig[4*np+1]=r2;
1330 zsig[4*np+2]=r3;
1331 zsig[4*np+3]=r4;
1332 sp_sig[2*izz ]=r1;
1333 sp_sig[2*izz+1]=r2;
1334 //if(mailbox[1]==1)fprintf( dmp,"\nA sp_sig %f %f %f %f",r1,r2,r3,r4);
1335 izz++;
1336 np=(np+1)&fftxn_mask;
1337 }
1338 }
1339 }
1340 else
1341 {
1342 if(sw_onechan)
1343 {
1344 while(np != nx)
1345 {
1346 k=spurp0+np*SPUR_WIDTH*twice_rxchan;
1347 t1=0;
1348 t2=0;
1349 for(i=0; i<SPUR_WIDTH-1; i+=2)
1350 {
1351 t1+=spur_power[i]*spur_table[2*i+k ]-
1352 spur_power[i+1]*spur_table[2*i+k+2];
1353 t2+=spur_power[i]*spur_table[2*i+k+1]-
1354 spur_power[i+1]*spur_table[2*i+k+3];
1355 }
1356 if((spur_location[spurno]&1)==1)
1357 {
1358 t1=-t1;
1359 t2=-t2;
1360 }
1361 zsig[2*np ]=t1;
1362 zsig[2*np+1]=t2;
1363 sp_sig[2*izz ]=t1;
1364 sp_sig[2*izz+1]=t2;
1365 izz++;
1366 np=(np+1)&fftxn_mask;
1367 }
1368 }
1369 else
1370 {
1371 while(np != nx)
1372 {
1373 k=spurp0+np*SPUR_WIDTH*twice_rxchan;
1374 t1=0;
1375 t2=0;
1376 t3=0;
1377 t4=0;
1378 for(i=0; i<SPUR_WIDTH-1; i+=2)
1379 {
1380 t1+=spur_power[i]*spur_table[4*i+k ]-
1381 spur_power[i+1]*spur_table[4*i+k+4];
1382 t2+=spur_power[i]*spur_table[4*i+k+1]-
1383 spur_power[i+1]*spur_table[4*i+k+5];
1384 t3+=spur_power[i]*spur_table[4*i+k+2]-
1385 spur_power[i+1]*spur_table[4*i+k+6];
1386 t4+=spur_power[i]*spur_table[4*i+k+3]-
1387 spur_power[i+1]*spur_table[4*i+k+7];
1388 }
1389 if((spur_location[spurno]&1)==1)
1390 {
1391 r1=-sp_c1*t1-sp_c2*t3-sp_c3*t4;
1392 r2=-sp_c1*t2-sp_c2*t4+sp_c3*t3;
1393 r3=-sp_c1*t3+sp_c2*t1-sp_c3*t2;
1394 r4=-sp_c1*t4+sp_c2*t2+sp_c3*t1;
1395 }
1396 else
1397 {
1398 r1=sp_c1*t1+sp_c2*t3+sp_c3*t4;
1399 r2=sp_c1*t2+sp_c2*t4-sp_c3*t3;
1400 r3=sp_c1*t3-sp_c2*t1+sp_c3*t2;
1401 r4=sp_c1*t4-sp_c2*t2-sp_c3*t1;
1402 }
1403 zsig[4*np ]=r1;
1404 zsig[4*np+1]=r2;
1405 zsig[4*np+2]=r3;
1406 zsig[4*np+3]=r4;
1407 sp_sig[2*izz ]=r1;
1408 sp_sig[2*izz+1]=r2;
1409 //if(mailbox[1]==1)fprintf( dmp,"\nB sp_sig %f %f %f %f",r1,r2,r3,r4);
1410 izz++;
1411 np=(np+1)&fftxn_mask;
1412 }
1413 }
1414 }
1415
1416 spur_phase_parameters();
1417 if(spur_ampl[spurno] < 3*spur_noise[spurno]/sqrt((float)(spur_speknum)))
1418 {
1419 // if(mailbox[1]==1)fprintf( dmp,"\nnoise test fail: %f < %f ",spur_ampl[spurno],3*spur_noise[spurno]/sqrt((float)(spur_speknum)));
1420 return -1;
1421 }
1422 i=verify_spur_pll();
1423 spur_avgd2[spurno]=spur_d2pha[spurno];
1424 return i;
1425 }
1426
verify_spur_pll(void)1427 int verify_spur_pll(void)
1428 {
1429 int iter, i, j, ind, knd;
1430 int izz;
1431 int ni, pnt;
1432 float sigspek[SPUR_WIDTH+2],errspek[SPUR_WIDTH+2];
1433 float t1,t2,t3,t4,r1,r2,r3,r4;
1434 float a1, a2, b1, b2, d1, d2;
1435 float average_rot,rot;
1436 float freq, d2err, d1err, d0err;
1437 float phase_slope,phase_curv;
1438 float *zsig, *z;
1439 short int *zmmx;
1440 int *uind;
1441 // Use the phases to transform the spur signal.
1442 // Ideally the signal should then be entirely in I.
1443 // Make average_rot the phase shift from transform to transform
1444 // based on the average power spectrum.
1445 // Note that this phase rotation is many turns.
1446 // We use average_rot to determine the number of whole turns
1447 // while the decimals are determined from the transform to transform
1448 // phase shifts we get from the phase first derivative.
1449 // When we know rot with accurate decimals, we can go back to
1450 // frequency (in fft bins) carrying the improved accuracy.
1451 // The decimals of the frequency are used to index spur_spectra
1452 // which contains normalised reference spectra with the same
1453 // frequency.
1454 // Use reference spectra and 90 degree phase shifted reference
1455 // spectra to get I and Q, normalised amplitudes for the
1456 // in phase and out of phase components of the spur signal.
1457 // If our phases were correct I is amplitude while Q is small
1458 // and contains the phase error.
1459 // The procedure is a PLL loop.
1460 // We will use Q to improve the phase function, our local oscillator.
1461 // The phase function is low pass filtered (PLL loop filter) which
1462 // gives the spur removal a narrow bandwidth.
1463 zsig=&spur_signal[twice_rxchan*max_fftxn*spurno];
1464 uind=&spur_ind[spurno*max_fftxn];
1465 freq=spur_freq[spurno];
1466 average_rot=freq*spur_freq_factor;
1467 d0err=0;
1468 d1err=0;
1469 d2err=0;
1470 iter=0;
1471 /*
1472 if(mailbox[1]==1)fprintf( dmp,"\n\n VERIFY ");
1473 if(mailbox[1]==1)
1474 {
1475 fprintf( dmp,"\n************ %d ********************",spur_location[spurno]);
1476 for(i=0; i<spur_speknum; i++)
1477 {
1478 fprintf( dmp,"\n%d %f %f angle %f",i,sp_sig[2*i],sp_sig[2*i+1],atan2(sp_sig[2*i+1],sp_sig[2*i])/PI_L);
1479 }
1480 }
1481 */
1482 restart:;
1483 iter++;
1484 a1=cos(spur_d0pha[spurno]);
1485 a2=sin(spur_d0pha[spurno]);
1486 b1=cos(spur_d1pha[spurno]);
1487 b2=sin(spur_d1pha[spurno]);
1488 d1=cos(spur_d2pha[spurno]);
1489 d2=sin(spur_d2pha[spurno]);
1490 phase_slope=spur_d1pha[spurno];
1491 phase_curv=spur_d2pha[spurno];
1492 ni=(ffts_na+fftxn_mask)&fftxn_mask;
1493 izz=spur_speknum-1;
1494 //if(mailbox[1]==1)fprintf( dmp,"\n-------------------------");
1495 while(izz >= 0)
1496 {
1497 rot=-0.5*phase_slope/PI_L;
1498 i=average_rot-rot+0.5;
1499 rot+=i;
1500 freq=rot/spur_freq_factor;
1501 j=(int)(freq)+2-spur_location[spurno]-SPUR_SIZE/2;
1502 if(j<0)j=0;
1503 j=1-j;
1504 if(j<0)j=0;
1505 ind=NO_OF_SPUR_SPECTRA*(freq-(int)(freq));
1506 if(ind==NO_OF_SPUR_SPECTRA)ind=NO_OF_SPUR_SPECTRA-1;
1507 ind=ind*SPUR_SIZE+j;
1508 uind[ni]=ind;
1509 knd=spurp0+ni*SPUR_WIDTH*twice_rxchan;
1510 if(sw_onechan)
1511 {
1512 r1=0;
1513 r2=0;
1514 if(swmmx_fft2)
1515 {
1516 for(i=0; i<SPUR_WIDTH; i++)
1517 {
1518 r1+=spur_spectra[ind+i]*spur_table_mmx[2*i+knd ];
1519 r2+=spur_spectra[ind+i]*spur_table_mmx[2*i+knd+1];
1520 }
1521 }
1522 else
1523 {
1524 for(i=0; i<SPUR_WIDTH; i++)
1525 {
1526 r1+=spur_spectra[ind+i]*spur_table[2*i+knd ];
1527 r2+=spur_spectra[ind+i]*spur_table[2*i+knd+1];
1528 }
1529 }
1530 }
1531 else
1532 {
1533 t1=0;
1534 t2=0;
1535 t3=0;
1536 t4=0;
1537 if(swmmx_fft2)
1538 {
1539 for(i=0; i<SPUR_WIDTH; i++)
1540 {
1541 t1+=spur_spectra[ind+i]*spur_table_mmx[4*i+knd ];
1542 t2+=spur_spectra[ind+i]*spur_table_mmx[4*i+knd+1];
1543 t3+=spur_spectra[ind+i]*spur_table_mmx[4*i+knd+2];
1544 t4+=spur_spectra[ind+i]*spur_table_mmx[4*i+knd+3];
1545 }
1546 }
1547 else
1548 {
1549 for(i=0; i<SPUR_WIDTH; i++)
1550 {
1551 t1+=spur_spectra[ind+i]*spur_table[4*i+knd ];
1552 t2+=spur_spectra[ind+i]*spur_table[4*i+knd+1];
1553 t3+=spur_spectra[ind+i]*spur_table[4*i+knd+2];
1554 t4+=spur_spectra[ind+i]*spur_table[4*i+knd+3];
1555 }
1556 }
1557 r1=sp_c1*t1+sp_c2*t3+sp_c3*t4;
1558 r2=sp_c1*t2+sp_c2*t4-sp_c3*t3;
1559 r3=sp_c1*t3-sp_c2*t1+sp_c3*t2;
1560 r4=sp_c1*t4-sp_c2*t2-sp_c3*t1;
1561 if((j^(spur_location[spurno]&1))==1)
1562 {
1563 r3=-r3;
1564 r4=-r4;
1565 }
1566 zsig[4*ni+2]=r3;
1567 zsig[4*ni+3]=r4;
1568 }
1569 if((j^(spur_location[spurno]&1))==1)
1570 {
1571 r1=-r1;
1572 r2=-r2;
1573 }
1574 zsig[twice_rxchan*ni ]=r1;
1575 zsig[twice_rxchan*ni+1]=r2;
1576 sp_sig[2*izz ]=r1*a1+r2*a2;
1577 sp_sig[2*izz+1]=r2*a1-r1*a2;
1578 ni=(ni+fftxn_mask)&fftxn_mask;
1579 izz--;
1580 r1=a1*b1+a2*b2;
1581 a2=a2*b1-a1*b2;
1582 a1=r1;
1583 r2=b1*d1+b2*d2;
1584 b2=b2*d1-b1*d2;
1585 b1=r2;
1586 phase_slope-=phase_curv;
1587 }
1588 spur_phase_parameters();
1589 //if(mailbox[1]==1)fprintf( dmp,"\nnoise test: %f < %f ",spur_ampl[spurno],3*spur_noise[spurno]/sqrt((float)(spur_speknum)));
1590
1591 if(spur_ampl[spurno] < 3*spur_noise[spurno]/sqrt((float)(spur_speknum)))
1592 {
1593 // if(mailbox[1]==1)fprintf( dmp,"\nnoise test: %f < %f ",spur_ampl[spurno],3*spur_noise[spurno]/sqrt((float)(spur_speknum)));
1594 return -1;
1595 }
1596 spur_d0pha[spurno]+=sp_d0;
1597 spur_d1pha[spurno]+=sp_d1;
1598 spur_d2pha[spurno]+=sp_d2;
1599 while(spur_d0pha[spurno] > PI_L)spur_d0pha[spurno]-=2*PI_L;
1600 while(spur_d0pha[spurno] <-PI_L)spur_d0pha[spurno]+=2*PI_L;
1601 while(spur_d1pha[spurno] > PI_L)spur_d1pha[spurno]-=2*PI_L;
1602 while(spur_d1pha[spurno] <-PI_L) spur_d1pha[spurno]+=2*PI_L;
1603 while(spur_d2pha[spurno] > PI_L)spur_d2pha[spurno]-=2*PI_L;
1604 while(spur_d2pha[spurno] <-PI_L)spur_d2pha[spurno]+=2*PI_L;
1605 phase_slope=spur_d1pha[spurno];
1606 phase_curv=spur_d2pha[spurno];
1607 rot=-0.5*phase_slope/PI_L;
1608 i=average_rot-rot+0.5;
1609 rot+=i;
1610 spur_freq[spurno]=rot/spur_freq_factor;
1611 /*
1612 if(mailbox[1]==1)
1613 {
1614 fprintf( dmp,"\n************ %d ********************",spur_location[spurno]);
1615 for(i=0; i<spur_speknum; i++)
1616 {
1617 fprintf( dmp,"\n%d %f %f angle %f",i,sp_sig[2*i],sp_sig[2*i+1],atan2(sp_sig[2*i+1],sp_sig[2*i])/PI_L);
1618 }
1619 fprintf( dmp,"\nd0 %f,%f d1 %f,%f d2 %f,%f",sp_d0,d0err,sp_d1,d1err,sp_d2,d2err);
1620 }
1621 */
1622 if(iter > 1 &&
1623 fabs(sp_d0) < 0.1 &&
1624 fabs(sp_d1) < 0.01 &&
1625 fabs(sp_d2) < 0.001 &&
1626 fabs(d0err) < 0.3 &&
1627 fabs(d1err) < 0.03 &&
1628 fabs(d2err) < 0.003 )
1629 {
1630 // Check the spectrum of what would be left when the spur is subtracted.
1631 /*
1632 if(mailbox[1]==1)
1633 {
1634 fprintf( dmp,"\n------------ %d - %f -----------",spur_location[spurno],spur_ampl[spurno]);
1635 for(i=0; i<spur_speknum; i++)
1636 {
1637 fprintf( dmp,"\n%d %f %f",i,sp_sig[2*i],sp_sig[2*i+1]);
1638 }
1639 }
1640 */
1641 for(i=0; i<SPUR_WIDTH+2; i++)
1642 {
1643 sigspek[i]=0;
1644 errspek[i]=0;
1645 }
1646 #define ZX 0.01
1647 average_rot=spur_freq[spurno]*spur_freq_factor;
1648 a1=spur_ampl[spurno]*cos(spur_d0pha[spurno]);
1649 a2=spur_ampl[spurno]*sin(spur_d0pha[spurno]);
1650 b1=cos(spur_d1pha[spurno]);
1651 b2=sin(spur_d1pha[spurno]);
1652 d1=cos(spur_d2pha[spurno]);
1653 d2=sin(spur_d2pha[spurno]);
1654 // if(mailbox[1]==1)fprintf( dmp,"\nphase: d0 %f d1 %f d2 %f",spur_d0pha[spurno],spur_d1pha[spurno],spur_d2pha[spurno]);
1655 phase_slope=spur_d1pha[spurno];
1656 phase_curv=spur_d2pha[spurno];
1657 ni=(ffts_na+fftxn_mask)&fftxn_mask;
1658 pnt=spur_location[spurno];
1659 izz=spur_speknum-1;
1660 while(izz >= 0)
1661 {
1662 rot=-0.5*phase_slope/PI_L;
1663 i=average_rot-rot+0.5;
1664 rot+=i;
1665 freq=rot/spur_freq_factor;
1666 j=(int)(freq)+2-pnt-SPUR_SIZE/2;
1667 if(j<0)j=0;
1668 j=1-j;
1669 if(j<0)j=0;
1670 ind=uind[ni];
1671 if(sw_onechan)
1672 {
1673 if((j^(spur_location[spurno]&1))==1)
1674 {
1675 t1=-a1;
1676 t2=-a2;
1677 }
1678 else
1679 {
1680 t1=a1;
1681 t2=a2;
1682 }
1683 if(swmmx_fft2)
1684 {
1685 zmmx=&fft2_short_int[twice_rxchan*(ni*fft2_size+pnt)];
1686 for(i=0; i<SPUR_WIDTH; i++)
1687 {
1688 sigspek[i+1]+=spur_spectra[ind+i]*spur_spectra[ind+i]*(t1*t1+t2*t2);
1689 errspek[i+1]+=pow(zmmx[2*i ]-spur_spectra[ind+i]*t1,2.0)+
1690 pow(zmmx[2*i+1]-spur_spectra[ind+i]*t2,2.0);
1691 }
1692 errspek[0]+=pow((float)(zmmx[-2]),2.0)+
1693 pow((float)(zmmx[-1]),2.0);
1694 errspek[SPUR_WIDTH+1]+=pow((float)(zmmx[2*(SPUR_WIDTH+1) ]),2.0)+
1695 pow((float)(zmmx[2*(SPUR_WIDTH+1)+1]),2.0);
1696 }
1697 else
1698 {
1699 z=&fftx[twice_rxchan*(ni*fftx_size+pnt)];
1700 for(i=0; i<SPUR_WIDTH; i++)
1701 {
1702 sigspek[i+1]+=spur_spectra[ind+i]*spur_spectra[ind+i]*(t1*t1+t2*t2);
1703 /*
1704 if(mailbox[1]==1)fprintf( dmp,"\nSUB %d %d: (%f,%.3f) (%f,%.3f) (%f,%.3f)",ni,i,
1705 ZX*sqrt(pow(z[2*i ],2.0)+pow(z[2*i+1],2.0)),atan2(z[2*i+1],z[2*i])/PI_L,
1706 ZX*sqrt(pow(spur_spectra[ind+i]*t1,2.0)+pow(spur_spectra[ind+i]*t2,2.0)),
1707 atan2(spur_spectra[ind+i]*t2,spur_spectra[ind+i]*t1)/PI_L,
1708 ZX*sqrt(pow(z[2*i ]-spur_spectra[ind+i]*t1,2.0)+
1709 pow(z[2*i+1]-spur_spectra[ind+i]*t2,2.0)),
1710 atan2(z[2*i+1]-spur_spectra[ind+i]*t2,z[2*i ]-spur_spectra[ind+i]*t1)/PI_L);
1711 */
1712 errspek[i+1]+=pow(z[2*i ]-spur_spectra[ind+i]*t1,2.0)+
1713 pow(z[2*i+1]-spur_spectra[ind+i]*t2,2.0);
1714 }
1715
1716 errspek[0]+=pow(z[-2],2.0)+pow(z[-1],2.0);
1717 errspek[SPUR_WIDTH+1]+=pow(z[2*(SPUR_WIDTH+1) ],2.0)+
1718 pow(z[2*(SPUR_WIDTH+1)+1],2.0);
1719 }
1720 }
1721 else
1722 {
1723 if((j^(spur_location[spurno]&1))==1)
1724 {
1725 t1=-sp_c1*a1;
1726 t2=-sp_c1*a2;
1727 t3=-sp_c2*a1+sp_c3*a2;
1728 t4=-sp_c2*a2-sp_c3*a1;
1729 }
1730 else
1731 {
1732 t1=sp_c1*a1;
1733 t2=sp_c1*a2;
1734 t3=sp_c2*a1-sp_c3*a2;
1735 t4=sp_c2*a2+sp_c3*a1;
1736 }
1737 if(swmmx_fft2)
1738 {
1739 zmmx=&fft2_short_int[twice_rxchan*(ni*fft2_size+pnt)];
1740 for(i=0; i<SPUR_WIDTH; i++)
1741 {
1742 sigspek[i+1]+=spur_spectra[ind+i]*spur_spectra[ind+i]*
1743 (t1*t1+t2*t2+t3*t3+t4*t4);
1744
1745
1746 errspek[i+1]+=pow(zmmx[4*i ]-spur_spectra[ind+i]*t1,2.0)+
1747 pow(zmmx[4*i+1]-spur_spectra[ind+i]*t2,2.0)+
1748 pow(zmmx[4*i+2]-spur_spectra[ind+i]*t3,2.0)+
1749 pow(zmmx[4*i+3]-spur_spectra[ind+i]*t4,2.0);
1750 }
1751 errspek[0]+=pow((float)(zmmx[-4]),2.0)+
1752 pow((float)(zmmx[-3]),2.0)+
1753 pow((float)(zmmx[-2]),2.0)+
1754 pow((float)(zmmx[-1]),2.0);
1755 errspek[SPUR_WIDTH+1]+=pow((float)(zmmx[4*(SPUR_WIDTH+1) ]),2.0)+
1756 pow((float)(zmmx[4*(SPUR_WIDTH+1)+1]),2.0)+
1757 pow((float)(zmmx[4*(SPUR_WIDTH+1)+2]),2.0)+
1758 pow((float)(zmmx[4*(SPUR_WIDTH+1)+3]),2.0);
1759 }
1760 else
1761 {
1762 z=&fftx[twice_rxchan*(ni*fftx_size+pnt)];
1763 for(i=0; i<SPUR_WIDTH; i++)
1764 {
1765 sigspek[i+1]+=spur_spectra[ind+i]*spur_spectra[ind+i]*
1766 (t1*t1+t2*t2+t3*t3+t4*t4);
1767 /*
1768 if(mailbox[1]==1)fprintf( dmp,"\nSUB %d %d:[(%f,%.3f)(%f,%.3f)] %f {%f,%f,%f,%f)",ni,i,
1769 ZX*z[4*i ],ZX*z[4*i+1], ZX*z[4*i+2],ZX*z[4*i+3],
1770
1771 ZX*spur_spectra[ind+i],
1772 ZX*spur_spectra[ind+i]*t1,
1773 ZX*spur_spectra[ind+i]*t2,
1774 ZX*spur_spectra[ind+i]*t3,
1775 ZX*spur_spectra[ind+i]*t4);
1776 */
1777 errspek[i+1]+=pow(z[4*i ]-spur_spectra[ind+i]*t1,2.0)+
1778 pow(z[4*i+1]-spur_spectra[ind+i]*t2,2.0)+
1779 pow(z[4*i+2]-spur_spectra[ind+i]*t3,2.0)+
1780 pow(z[4*i+3]-spur_spectra[ind+i]*t4,2.0);
1781 }
1782 errspek[0]+=pow(z[-4],2.0)+pow(z[-3],2.0)+pow(z[-2],2.0)+pow(z[-1],2.0);
1783 errspek[SPUR_WIDTH+1]+=pow(z[4*(SPUR_WIDTH+1) ],2.0)+
1784 pow(z[4*(SPUR_WIDTH+1)+1],2.0)+
1785 pow(z[4*(SPUR_WIDTH+1)+2],2.0)+
1786 pow(z[4*(SPUR_WIDTH+1)+3],2.0);
1787 }
1788 }
1789 ni=(ni+fftxn_mask)&fftxn_mask;
1790 izz--;
1791 r1=a1*b1+a2*b2;
1792 a2=a2*b1-a1*b2;
1793 a1=r1;
1794 r2=b1*d1+b2*d2;
1795 b2=b2*d1-b1*d2;
1796 b1=r2;
1797 phase_slope-=phase_curv;
1798 }
1799 // Compute the RMS noise amplitude within the fftx bandwidth
1800 // while converting errspek[] to amplitude.
1801 r1=0;
1802 for(i=0; i<SPUR_WIDTH+2; i++)
1803 {
1804 errspek[i]/=spur_speknum;
1805 r1+=errspek[i];
1806 errspek[i]=sqrt(errspek[i]);
1807 //if(mailbox[1]==1) fprintf( dmp,"\n%d N %f",i,errspek[i]);
1808 }
1809 r1=sqrt(r1/(SPUR_WIDTH+2));
1810 // Check whether the noise floor is flat.
1811 r2=0;
1812 for(i=0; i<SPUR_WIDTH+2; i++)
1813 {
1814 r2+=pow(errspek[i]-r1,2.0);
1815 }
1816 r2=sqrt(r2/(SPUR_WIDTH+2));
1817 // r1 is now the amplitude of the noise floor.
1818 // r2 is the standard deviation of the noise floor.
1819 // Spurs may be modulated to some extent, that would lead
1820 // to a non-flat noise floor so we have to accept some
1821 // contribution due to modulation.
1822 // The RMS
1823 //if(mailbox[1]==1)fprintf( dmp,"\nampl %f noi %f",spur_ampl[spurno],spur_noise[spurno]);
1824 //if(mailbox[1]==1)fprintf( dmp,"\nerr: %f spread %f S/N %f",spur_noise[spurno]/r1,r2/r1,spur_ampl[spurno]/r1);
1825 //if(mailbox[1]==1)fprintf( dmp,"\ntest:%f > %f+%f",r2/r1,0.5/sqrt((int)(spur_speknum)),0.02*spur_ampl[spurno]/r1);
1826 if( r2 > 0.5*r1/sqrt((int)(spur_speknum))+0.02*spur_ampl[spurno])
1827 {
1828 // if(mailbox[1]==1)fprintf( dmp,"\nSPREAD failed");
1829 return -1;
1830 }
1831
1832
1833
1834 return 0;
1835 }
1836 d0err=sp_d0;
1837 d1err=sp_d1;
1838 d2err=sp_d2;
1839
1840 if(iter < 5 ) goto restart;
1841 return -1;
1842 }
1843