1 /**
2 @file
3 @ingroup processing
4 */
5 #include <stdio.h>
6 #include <bpm/bpm_units.h>
7 #include <bpm/bpm_messages.h>
8 #include <bpm/bpm_process.h>
9 #include <bpm/bpm_dsp.h>
10
11
process_waveform(doublewf_t * signal,bpmconf_t * bpm,bpmproc_t * proc,bpmproc_t * trig,unsigned int mode)12 int process_waveform( doublewf_t *signal, bpmconf_t *bpm, bpmproc_t *proc, bpmproc_t *trig,
13 unsigned int mode ) {
14
15 int have_fft_pars = FALSE;
16 int have_fit_pars = FALSE;
17 double frequency, tdecay, amp, phase;
18 char msg[128];
19 double t0, tSample;
20
21 if ( ! bpm || ! signal || ! proc ) {
22 bpm_error( "Invalid pointer arguments in process_waveform(...)",
23 __FILE__, __LINE__ );
24 return BPM_FAILURE;
25 }
26
27 // reset the success values for this pulse...
28 proc->fft_success = FALSE;
29 proc->fit_success = FALSE;
30 proc->ddc_success = FALSE;
31
32 // do we have a signal ?
33 if ( ! signal ) {
34 sprintf( msg, "No signal present for BPM %s in process_waveform(...)", bpm->name );
35 bpm_error( msg, __FILE__, __LINE__ );
36 return BPM_FAILURE;
37 }
38
39 /* ------------------------------- check for saturation ------------------------------------ */
40 proc->saturated = check_saturation( signal, bpm->digi_nbits, &(proc->iunsat) );
41
42
43 /* ------------------------------- subtract the pedestal ----------------------------------- */
44 // determing voltage offset and amplitude noise in adc channels (pedestal)
45 if ( get_pedestal( signal, 20, &(proc->voltageoffset) ,&(proc->ampnoise) ) == BPM_FAILURE ) {
46 sprintf( msg, "Error getting pedestal of BPM %s in process_waveform(...)", bpm->name );
47 bpm_error( msg, __FILE__, __LINE__ );
48 return BPM_FAILURE;
49 }
50
51 // subtract the pedestal
52 doublewf_bias( -proc->voltageoffset, signal );
53
54
55 /* ----------------------------------- set the t0 ------------------------------------------ */
56 if ( trig ) {
57 // internal clock, we have a trigger
58 proc->t0 = trig->t0;
59 } else {
60 // external clock, no trigger specified, assume t0 is set in configuration
61 proc->t0 = bpm->t0;
62 }
63
64 /* ------------------------------ check whether to do FFT ? -------------------------------- */
65 if ( mode & PROC_DO_FFT ) {
66
67 // compute the ft
68 if ( fft_waveform( signal, proc->ft ) == BPM_FAILURE ) {
69 sprintf( msg, "Could not perform fft for BPM %s in process_waveform(...)", bpm->name );
70 bpm_warning( msg, __FILE__, __LINE__ );
71
72 } else {
73 proc->fft_success = TRUE;
74
75 if ( mode & PROC_FIT_FFT ) {
76 // fft done, fit for frequency and tdecay
77 if( fit_fft( proc->ft, &(proc->fft_freq), &(proc->fft_tdecay), NULL, NULL )
78 == BPM_FAILURE ) {
79 sprintf( msg, "Could not fit the FFT for BPM %s in process_waveform(...)", bpm->name );
80 bpm_warning( msg, __FILE__, __LINE__ );
81 } else {
82 have_fft_pars = TRUE;
83 }
84
85 }
86
87 }
88
89 } /* if ( mode & PROC_DO_FFT ) */
90
91
92 /* ------------------------------ check whether to do FIT ? -------------------------------- */
93 if ( mode & PROC_DO_FIT ) {
94
95 /*
96 This stuff is a bit shaky... need to clean & test in more detail if needed,
97 rely on ROOT/Minuit based fitting for this !!!
98 */
99 bpm_warning( "Libbpm waveform fitting is highly EXPERIMENTAL and contains bugs!",
100 __FILE__, __LINE__ );
101
102 // initial parameters
103 frequency = tdecay = amp = phase = 0.;
104
105 // Initial parameter for the frequency ::
106 if ( proc->fit_freq <= 0.* _MHz__ || proc->fit_freq > ( bpm->digi_freq/2.) ) {
107 // BM. 11.02.2008
108 // cut this out, don't want it to depend on the rf_LOfreq just yet...
109 /*
110 if ( ( bpm->rf_LOfreq > 0. ) && ( bpm->cav_freq > 0. ) ) {
111 freq = ABS( bpm->cav_freq - bpm->rf_LOfreq );
112 } else
113 */
114 frequency = 20.0 * _MHz__;
115 } else frequency = proc->fit_freq; // use last fitted value
116
117 // Initial parameter for the decay constant
118 if ( proc->fit_tdecay <= 10.* _nsec__ || proc->fit_tdecay > 1.0 * _sec__ ) {
119 if ( bpm->cav_decaytime > 0. ) tdecay = bpm->cav_decaytime;
120 else tdecay = 0.2 * _sec__;
121 } else tdecay = proc->fit_tdecay; // use last fitted value
122
123 // probably can think of something more clever here...
124 if( ABS( proc->fit_amp ) < 100000. ) amp = proc->fit_amp; else amp = 3000.;
125 phase = proc->fit_phase;
126
127 // if we have down the FFT and got the frequency from there, use that as an inital parameter
128 // so overwrite the previous
129 if( have_fft_pars ) {
130 frequency = proc->fft_freq;
131 tdecay = proc->fft_tdecay;
132 }
133
134 // Slight hack to sort out saturation and trigger offsets
135 // Change t0, fit, then extrapolate back in both amplitude and phase
136 sample_to_time( bpm->digi_freq, bpm->digi_nsamples, proc->iunsat, &t0 );
137 t0 = ( t0 > ( proc->t0 + bpm->fit_tOffset ) ? t0 : proc->t0 + bpm->fit_tOffset);
138
139 if ( fit_waveform( signal, t0, frequency, tdecay, amp, phase,
140 &(proc->fit_freq), &(proc->fit_tdecay), &(proc->fit_amp) ,
141 &(proc->fit_phase) ) == BPM_FAILURE ) {
142 sprintf( msg, "Could not fit for BPM %s in process_waveform(...)", bpm->name );
143 bpm_warning( msg, __FILE__, __LINE__ );
144 } else {
145 have_fit_pars = TRUE;
146 proc->fit_success = TRUE;
147
148 // Extrapolate back to the real t0 (proc->t0) !!
149 proc->fit_amp *= exp( ( t0 - proc->t0 ) / proc->fit_tdecay );
150 proc->fit_phase -= (t0 - proc->t0) * proc->fit_freq * 2. * PI;
151 norm_phase( &(proc->fit_phase) );
152 }
153
154 } /* if ( mode & PROC_DO_FIT ) */
155
156
157 /* ------------------------------ check whether to do DDC ? -------------------------------- */
158 if ( mode & PROC_DO_DDC ) {
159
160 // initial parameters
161 frequency = tdecay = 0.;
162
163 // where does the ddc get its freq and tdecay parameters from ?
164 if ( mode & PROC_DDC_FITFREQ ) {
165 if ( have_fit_pars ) frequency = proc->fit_freq;
166 else frequency = bpm->ddc_freq; // asked for fit but it failed :(
167 } else if ( mode & PROC_DDC_FFTFREQ ) {
168 if ( have_fft_pars ) frequency = proc->fft_freq;
169 else frequency = bpm->ddc_freq; // asked for fft but it failed :(
170 } else {
171 // use default
172 frequency = bpm->ddc_freq;
173 }
174
175 // same for tdecay
176 if ( mode & PROC_DDC_FITTDECAY ) {
177 if ( have_fit_pars ) tdecay = proc->fit_tdecay;
178 else tdecay = bpm->ddc_tdecay; // asked for fit but it failed :(
179 } else if ( mode & PROC_DDC_FFTTDECAY ) {
180 if ( have_fft_pars ) tdecay = proc->fft_tdecay;
181 else tdecay = bpm->ddc_tdecay; // asked for fft but it failed :(
182 } else {
183 // use default
184 tdecay = bpm->ddc_tdecay;
185 }
186
187 // determine the sample time, handle saturation here...
188 if ( proc->saturated == 1 ) {
189 // We have saturation, get the index of the last unsaturated sample
190 sample_to_time( bpm->digi_freq, bpm->digi_nsamples, proc->iunsat, &(proc->ddc_tSample) );
191 // normally we should add here something which takes into account the bandwidth
192 // of the filter since we shift the filter window, let's just take
193 } else {
194 // The normall case, the t0 (trigger or fixed + the offset)
195 proc->ddc_tSample = proc->t0 + bpm->ddc_tOffset;
196 }
197
198 // convert to sample index !
199 time_to_sample( bpm->digi_freq, bpm->digi_nsamples, proc->ddc_tSample, &(proc->ddc_iSample) );
200
201 if ( mode & PROC_DDC_FULL ) {
202 // if we must to the full DDC, do it and sample afterwards
203 if ( ddc_waveform( signal, frequency, bpm->ddc_filter, proc->dc,
204 bpm->ddc_buffer_re, bpm->ddc_buffer_im ) == BPM_FAILURE ) {
205 sprintf( msg, "Could not ddc BPM %s waveform in process_waveform(...)", bpm->name );
206 bpm_warning( msg, __FILE__, __LINE__ );
207 proc->ddc_amp = 0.;
208 proc->ddc_phase = 0.;
209
210 } else {
211 proc->ddc_success = TRUE;
212
213 // Now sample the ddc waveform at the iSample requested, extrapolate amplitude
214 proc->ddc_amp = exp( ( proc->ddc_tSample - proc->t0 ) / tdecay ) *
215 c_abs( proc->dc->wf[ proc->ddc_iSample ] );
216 proc->ddc_phase = c_arg( proc->dc->wf[ proc->ddc_iSample ] );
217 norm_phase( &(proc->ddc_phase) );
218 }
219
220 } else {
221 // Execute the "short" ddc algorithm to save cpu cycles
222 if ( ddc_sample_waveform( signal, frequency, bpm->ddc_filter, proc->ddc_iSample,
223 proc->t0, tdecay, &(proc->ddc_amp), &(proc->ddc_phase),
224 bpm->ddc_buffer_re, bpm->ddc_buffer_im ) == BPM_FAILURE ) {
225 sprintf( msg, "Could not ddc_sample BPM %s waveform in process_waveform(...)", bpm->name );
226 bpm_warning( msg, __FILE__, __LINE__ );
227 } else {
228 proc->ddc_success = TRUE;
229 norm_phase( &(proc->ddc_phase) );
230 }
231 }
232
233 } /* if ( mode & PROC_DO_DDC ) */
234
235 return BPM_SUCCESS;
236 }
237