Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
/*
* Copyright (C) 2019 Aaron Viets <aaron.viets@ligo.org>
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License along
* with this program; if not, write to the Free Software Foundation, Inc.,
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
*/
/*
* ============================================================================
*
* Preamble
*
* ============================================================================
*/
/*
* stuff from C
*/
#include <string.h>
#include <math.h>
#include <complex.h>
/*
* stuff from gobject/gstreamer
*/
#include <glib.h>
Aaron Viets
committed
#include <glib/gprintf.h>
#include <gst/gst.h>
#include <gst/audio/audio.h>
#include <gst/base/gstbasetransform.h>
/*
* our own stuff
*/
#include <gstlal/gstlal_audio_info.h>
#include <gstlal_sensingtdcfs.h>
/*
* ============================================================================
*
* GStreamer Boiler Plate
*
* ============================================================================
*/
#define INCAPS \
"audio/x-raw, " \
"format = (string) {"GST_AUDIO_NE(Z64)", "GST_AUDIO_NE(Z128)"}, " \
"rate = " GST_AUDIO_RATE_RANGE ", " \
"channels = (int) [1, MAX], " \
"layout = (string) interleaved, " \
"channel-mask = (bitmask) 0"
#define OUTCAPS \
"audio/x-raw, " \
"format = (string) {"GST_AUDIO_NE(F32)", "GST_AUDIO_NE(F64)"}, " \
"rate = " GST_AUDIO_RATE_RANGE ", " \
"channels = (int) [2, MAX], " \
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
"layout = (string) interleaved, " \
"channel-mask = (bitmask) 0"
static GstStaticPadTemplate sink_factory = GST_STATIC_PAD_TEMPLATE(
GST_BASE_TRANSFORM_SINK_NAME,
GST_PAD_SINK,
GST_PAD_ALWAYS,
GST_STATIC_CAPS(INCAPS)
);
static GstStaticPadTemplate src_factory = GST_STATIC_PAD_TEMPLATE(
GST_BASE_TRANSFORM_SRC_NAME,
GST_PAD_SRC,
GST_PAD_ALWAYS,
GST_STATIC_CAPS(OUTCAPS)
);
G_DEFINE_TYPE(
GSTLALSensingTDCFs,
gstlal_sensingtdcfs,
GST_TYPE_BASE_TRANSFORM
);
/*
* ============================================================================
*
* Utilities
*
* ============================================================================
*/
/*
* set the metadata on an output buffer
*/
static void set_metadata(GSTLALSensingTDCFs *element, GstBuffer *buf, guint64 outsamples, gboolean gap) {
GST_BUFFER_OFFSET(buf) = element->next_out_offset;
element->next_out_offset += outsamples;
GST_BUFFER_OFFSET_END(buf) = element->next_out_offset;
GST_BUFFER_TIMESTAMP(buf) = element->t0 + gst_util_uint64_scale_int_round(GST_BUFFER_OFFSET(buf) - element->offset0, GST_SECOND, element->rate);
GST_BUFFER_DURATION(buf) = element->t0 + gst_util_uint64_scale_int_round(GST_BUFFER_OFFSET_END(buf) - element->offset0, GST_SECOND, element->rate) - GST_BUFFER_TIMESTAMP(buf);
GST_BUFFER_FLAG_UNSET(buf, GST_BUFFER_FLAG_GAP);
if(G_UNLIKELY(element->need_discont)) {
GST_BUFFER_FLAG_SET(buf, GST_BUFFER_FLAG_DISCONT);
element->need_discont = FALSE;
}
if(gap)
GST_BUFFER_FLAG_SET(buf, GST_BUFFER_FLAG_GAP);
else
GST_BUFFER_FLAG_UNSET(buf, GST_BUFFER_FLAG_GAP);
}
/*
* Macro functions to solve for kappa_C using three different sensing function models: 0, 1, and 2
*/
#define DEFINE_KAPPA_C_0(DTYPE, F_OR_NOT) \
Aaron Viets
committed
void kappa_C_0_ ## DTYPE(DTYPE Cinv_tdep1_real, DTYPE Cinv_tdep1_imag, DTYPE Cinv_tdep2_real, DTYPE Cinv_tdep2_imag, DTYPE f1, DTYPE f2, DTYPE *kc1, DTYPE *kc2, DTYPE *kc3, DTYPE *kc4, char *debug_file) { \
DTYPE H1, H2, I1, I2, Xi, Zeta; \
/*
* Due to numerical instabilities, we typecast to long double precision
* here to avoid potentially large numerical errors.
*/ \
complex long double Q0, S0, complex_kc1, complex_kc2, complex_kc3, complex_kc4; \
long double a, b, c, d, e, Delta0, Delta1, p, q; \
H1 = -f1 * f1 * Cinv_tdep1_real; \
H2 = -f2 * f2 * Cinv_tdep2_real; \
I1 = -f1 * f1 * Cinv_tdep1_imag; \
I2 = -f2 * f2 * Cinv_tdep2_imag; \
Xi = (f2 * f2 * H1 - f1 * f1 * H2) / (f1 * f1 - f2 * f2); \
Zeta = I1 / f1 - I2 / f2; \
a = (long double) f2 * f2 * Xi * (H2 + Xi) * Zeta * Zeta; \
b = (long double) f2 * (f2 * f2 - f1 * f1) * (H2 + Xi) * Zeta * I2 + pow ## F_OR_NOT(f2, 4) * (H2 + 2 * Xi) * Zeta * Zeta; \
c = (long double) pow ## F_OR_NOT(f2, 3) * (f2 * f2 - f1 * f1) * Zeta * I2 + pow ## F_OR_NOT(f2 * f2 - f1 * f1, 2) * pow ## F_OR_NOT(H2 + Xi, 2) + pow ## F_OR_NOT(f2, 6) * Zeta * Zeta; \
d = (long double) 2 * f2 * f2 * pow ## F_OR_NOT(f2 * f2 - f1 * f1, 2) * (H2 + Xi); \
e = (long double) pow ## F_OR_NOT(f2, 4) * pow ## F_OR_NOT(f2 * f2 - f1 * f1, 2); \
Delta0 = c * c - 3 * b * d + 12 * a * e; \
Delta1 = 2 * powl(c, 3) - 9 * b * c * d + 27 * b * b * e + 27 * a * d * d - 72 * a * c * e; \
p = (8 * a * c - 3 * b * b) / (8 * a * a); \
q = (powl(b, 3) - 4 * a * b * c + 8 * a * a * d) / (8 * powl(a, 3)); \
Q0 = cpowl(0.5 * ((complex long double) Delta1 + cpowl((complex long double) (Delta1 * Delta1 - 4 * powl(Delta0, 3)), 0.5)), 1.0 / 3.0); \
S0 = 0.5 * cpowl((-2 * p) / 3.0 + (Q0 + Delta0 / Q0) / (3 * a), 0.5); \
\
complex_kc1 = -b / (4 * a) - S0 + 0.5 * cpowl(-4 * S0 * S0 - 2 * p + q / S0, 0.5); \
complex_kc2 = -b / (4 * a) - S0 - 0.5 * cpowl(-4 * S0 * S0 - 2 * p + q / S0, 0.5); \
complex_kc3 = -b / (4 * a) + S0 + 0.5 * cpowl(-4 * S0 * S0 - 2 * p - q / S0, 0.5); \
complex_kc4 = -b / (4 * a) + S0 - 0.5 * cpowl(-4 * S0 * S0 - 2 * p - q / S0, 0.5); \
\
/* DEBUGGING */ \
Aaron Viets
committed
if(debug_file) { \
FILE *fp; \
fp = fopen(debug_file, "a"); \
complex long double fcc1 = (f2 * f2 - f1 * f1) / (complex_kc1 * (f2 * Cinv_tdep2_imag - f1 * Cinv_tdep1_imag)); \
complex long double fs_squared1 = complex_kc1 * (Cinv_tdep2_real - Cinv_tdep1_real) / (1.0 / (f2 * f2) - 1.0 / (f1 * f1)); \
complex long double fs_over_Q1 = fcc1 * (complex_kc1 * Cinv_tdep1_real - 1.0 - fs_squared1 / (f1 * f1)); \
complex long double fcc2 = (f2 * f2 - f1 * f1) / (complex_kc2 * (f2 * Cinv_tdep2_imag - f1 * Cinv_tdep1_imag)); \
complex long double fs_squared2 = complex_kc2 * (Cinv_tdep2_real - Cinv_tdep1_real) / (1.0 / (f2 * f2) - 1.0 / (f1 * f1)); \
complex long double fs_over_Q2 = fcc2 * (complex_kc2 * Cinv_tdep1_real - 1.0 - fs_squared2 / (f1 * f1)); \
complex long double fcc3 = (f2 * f2 - f1 * f1) / (complex_kc3 * (f2 * Cinv_tdep2_imag - f1 * Cinv_tdep1_imag)); \
complex long double fs_squared3 = complex_kc3 * (Cinv_tdep2_real - Cinv_tdep1_real) / (1.0 / (f2 * f2) - 1.0 / (f1 * f1)); \
complex long double fs_over_Q3 = fcc3 * (complex_kc3 * Cinv_tdep1_real - 1.0 - fs_squared3 / (f1 * f1)); \
complex long double fcc4 = (f2 * f2 - f1 * f1) / (complex_kc4 * (f2 * Cinv_tdep2_imag - f1 * Cinv_tdep1_imag)); \
complex long double fs_squared4 = complex_kc4 * (Cinv_tdep2_real - Cinv_tdep1_real) / (1.0 / (f2 * f2) - 1.0 / (f1 * f1)); \
complex long double fs_over_Q4 = fcc4 * (complex_kc4 * Cinv_tdep1_real - 1.0 - fs_squared4 / (f1 * f1)); \
Aaron Viets
committed
long double Delta = 256 * a * a * a * e * e * e - 192 * a * a * b * d * e * e - 128 * a * a * c * c * e * e + 144 * a * a * c * d * d * e - 27 * a * a * d * d * d * d + 144 * a * b * b * c * e * e - 6 * a * b * b * d * d * e - 80 * a * b * c * c * d * e + 18 * a * b * c * d * d * d + 16 * a * c * c * c * c * e - 4 * a * c * c * c * d * d - 27 * b * b * b * b * e * e + 18 * b * b * b * c * d * e - 4 * b * b * b * d * d * d - 4 * b * b * c * c * c * e + b * b * c * c * d * d; \
long double P = 8 * a * c - 3 * b * b; \
long double D = 64 * a * a * a * e - 16 * a * a * c * c + 16 * a * b * b * c - 16 * a * a * b * d - 3 * b * b * b * b; \
\
Aaron Viets
committed
complex DTYPE condition1 = Cinv_tdep1_real + I * Cinv_tdep1_imag; \
complex DTYPE condition2 = Cinv_tdep2_real + I * Cinv_tdep2_imag; \
complex long double solution11 = (1 + I * f1 / fcc1) / complex_kc1 * (f1 * f1 + fs_squared1 - I * f1 * fs_over_Q1) / (f1 * f1); \
complex long double solution12 = (1 + I * f2 / fcc1) / complex_kc1 * (f2 * f2 + fs_squared1 - I * f2 * fs_over_Q1) / (f2 * f2); \
complex long double solution21 = (1 + I * f1 / fcc2) / complex_kc2 * (f1 * f1 + fs_squared2 - I * f1 * fs_over_Q2) / (f1 * f1); \
complex long double solution22 = (1 + I * f2 / fcc2) / complex_kc2 * (f2 * f2 + fs_squared2 - I * f2 * fs_over_Q2) / (f2 * f2); \
complex long double solution31 = (1 + I * f1 / fcc3) / complex_kc3 * (f1 * f1 + fs_squared3 - I * f1 * fs_over_Q3) / (f1 * f1); \
complex long double solution32 = (1 + I * f2 / fcc3) / complex_kc3 * (f2 * f2 + fs_squared3 - I * f2 * fs_over_Q3) / (f2 * f2); \
complex long double solution41 = (1 + I * f1 / fcc4) / complex_kc4 * (f1 * f1 + fs_squared4 - I * f1 * fs_over_Q4) / (f1 * f1); \
complex long double solution42 = (1 + I * f2 / fcc4) / complex_kc4 * (f2 * f2 + fs_squared4 - I * f2 * fs_over_Q4) / (f2 * f2); \
long double diff1 = cabsl(1.0 - condition1 / solution11) + cabsl(1.0 - condition2 / solution12); \
long double diff2 = cabsl(1.0 - condition1 / solution21) + cabsl(1.0 - condition2 / solution22); \
long double diff3 = cabsl(1.0 - condition1 / solution31) + cabsl(1.0 - condition2 / solution32); \
long double diff4 = cabsl(1.0 - condition1 / solution41) + cabsl(1.0 - condition2 / solution42); \
complex long double zero1 = a * complex_kc1 * complex_kc1 * complex_kc1 * complex_kc1 + b * complex_kc1 * complex_kc1 * complex_kc1 + c * complex_kc1 * complex_kc1 + d * complex_kc1 + e; \
complex long double zero2 = a * complex_kc2 * complex_kc2 * complex_kc2 * complex_kc2 + b * complex_kc2 * complex_kc2 * complex_kc2 + c * complex_kc2 * complex_kc2 + d * complex_kc2 + e; \
complex long double zero3 = a * complex_kc3 * complex_kc3 * complex_kc3 * complex_kc3 + b * complex_kc3 * complex_kc3 * complex_kc3 + c * complex_kc3 * complex_kc3 + d * complex_kc3 + e; \
complex long double zero4 = a * complex_kc4 * complex_kc4 * complex_kc4 * complex_kc4 + b * complex_kc4 * complex_kc4 * complex_kc4 + c * complex_kc4 * complex_kc4 + d * complex_kc4 + e; \
\
Aaron Viets
committed
g_fprintf(fp, "Calculation Details\nDelta=%Le\nP=%Le\nD=%Le\n", Delta, P, D); \
g_fprintf(fp, "kc1=%Le+1j*%Le\nfcc1=%Le+1j*%Le\nfs_squared1=%Le+1j*%Le\nfs_over_Q1=%Le+1j*%Le\n", creall(complex_kc1), cimagl(complex_kc1), creall(fcc1), cimagl(fcc1), creall(fs_squared1), cimagl(fs_squared1), creall(fs_over_Q1), cimagl(fs_over_Q1)); \
g_fprintf(fp, "kc2=%Le+1j*%Le\nfcc2=%Le+1j*%Le\nfs_squared2=%Le+1j*%Le\nfs_over_Q2=%Le+1j*%Le\n", creall(complex_kc2), cimagl(complex_kc2), creall(fcc2), cimagl(fcc2), creall(fs_squared2), cimagl(fs_squared2), creall(fs_over_Q2), cimagl(fs_over_Q2)); \
g_fprintf(fp, "kc3=%Le+1j*%Le\nfcc3=%Le+1j*%Le\nfs_squared3=%Le+1j*%Le\nfs_over_Q3=%Le+1j*%Le\n", creall(complex_kc3), cimagl(complex_kc3), creall(fcc3), cimagl(fcc3), creall(fs_squared3), cimagl(fs_squared3), creall(fs_over_Q3), cimagl(fs_over_Q3)); \
g_fprintf(fp, "kc4=%Le+1j*%Le\nfcc4=%Le+1j*%Le\nfs_squared4=%Le+1j*%Le\nfs_over_Q4=%Le+1j*%Le\n", creall(complex_kc4), cimagl(complex_kc4), creall(fcc4), cimagl(fcc4), creall(fs_squared4), cimagl(fs_squared4), creall(fs_over_Q4), cimagl(fs_over_Q4)); \
if(diff1 > 1e-12) \
Aaron Viets
committed
g_fprintf(fp, "Solution 1 failed: diff=%Le\n", diff1); \
else \
g_fprintf(fp, "diff1=%Le\n", diff1); \
if(diff2 > 1e-12) \
Aaron Viets
committed
g_fprintf(fp, "Solution 2 failed: diff=%Le\n", diff2); \
else \
g_fprintf(fp, "diff2=%Le\n", diff2); \
if(diff3 > 1e-12) \
Aaron Viets
committed
g_fprintf(fp, "Solution 3 failed: diff=%Le\n", diff3); \
else \
g_fprintf(fp, "diff3=%Le\n", diff3); \
if(diff4 > 1e-12) \
Aaron Viets
committed
g_fprintf(fp, "Solution 4 failed: diff=%Le\n", diff4); \
else \
g_fprintf(fp, "diff4=%Le\n", diff4); \
g_fprintf(fp, "a=%Le b=%Le c=%Le d=%Le e=%Le\n", a, b, c, d, e); \
g_fprintf(fp, "Solution 1: 0 = %Le+I%Le\n", creall(zero1), cimagl(zero1)); \
g_fprintf(fp, "Solution 2: 0 = %Le+I%Le\n", creall(zero2), cimagl(zero2)); \
g_fprintf(fp, "Solution 3: 0 = %Le+I%Le\n", creall(zero3), cimagl(zero3)); \
g_fprintf(fp, "Solution 4: 0 = %Le+I%Le\n\n", creall(zero4), cimagl(zero4)); \
fclose(fp); \
} \
/* END DEBUGGGING */ \
\
/* Any solution that is clearly complex or negative is not the solution we want */ \
if(creall(complex_kc1) > 0.0 && fabsl(cimagl(complex_kc1) / creall(complex_kc1)) < 1e-3) \
*kc1 = (DTYPE) creall(complex_kc1); \
if(creall(complex_kc2) > 0.0 && fabsl(cimagl(complex_kc2) / creall(complex_kc2)) < 1e-3) \
*kc2 = (DTYPE) creall(complex_kc2); \
if(creall(complex_kc3) > 0.0 && fabsl(cimagl(complex_kc3) / creall(complex_kc3)) < 1e-3) \
*kc3 = (DTYPE) creall(complex_kc3); \
if(creall(complex_kc4) > 0.0 && fabsl(cimagl(complex_kc4) / creall(complex_kc4)) < 1e-3) \
*kc4 = (DTYPE) creall(complex_kc4); \
}
DEFINE_KAPPA_C_0(float, f);
DEFINE_KAPPA_C_0(double, );
/*
* Macro functions to solve for f_cc using three different sensing function models: 0, 1, and 2
*/
#define DEFINE_F_CC_0(DTYPE) \
DTYPE f_cc_0_ ## DTYPE(DTYPE Cinv_tdep1_imag, DTYPE Cinv_tdep2_imag, DTYPE f1, DTYPE f2, DTYPE kappa_C) { \
return (f2 * f2 - f1 * f1) / (kappa_C * (f2 * Cinv_tdep2_imag - f1 * Cinv_tdep1_imag)); \
}
DEFINE_F_CC_0(float);
DEFINE_F_CC_0(double);
/*
* Macro functions to solve for f_s^2 using three different sensing function models: 0, 1, and 2
*/
#define DEFINE_F_S_SQUARED_0(DTYPE) \
DTYPE f_s_squared_0_ ## DTYPE(DTYPE Cinv_tdep1_real, DTYPE Cinv_tdep2_real, DTYPE f1, DTYPE f2, DTYPE kappa_C) { \
return kappa_C * (Cinv_tdep2_real - Cinv_tdep1_real) / (1.0 / (f2 * f2) - 1.0 / (f1 * f1)); \
}
DEFINE_F_S_SQUARED_0(float);
DEFINE_F_S_SQUARED_0(double);
/*
* Macro functions to solve for f_s / Q using three different sensing function models: 0, 1, and 2
*/
#define DEFINE_FS_OVER_Q_0(DTYPE) \
DTYPE f_s_over_Q_0_ ## DTYPE(DTYPE Cinv_tdep1_real, DTYPE f1, DTYPE kappa_C, DTYPE f_cc, DTYPE f_s_squared) { \
return f_cc * (kappa_C * Cinv_tdep1_real - 1.0 - f_s_squared / (f1 * f1)); \
}
DEFINE_FS_OVER_Q_0(float);
DEFINE_FS_OVER_Q_0(double);
#define DEFINE_FIND_BEST_SOLUTION(DTYPE, F_OR_NOT) \
Aaron Viets
committed
guint find_best_solution_ ## DTYPE(DTYPE kc1, DTYPE kc2, DTYPE kc3, DTYPE kc4, DTYPE fcc1, DTYPE fcc2, DTYPE fcc3, DTYPE fcc4, DTYPE fs_squared1, DTYPE fs_squared2, DTYPE fs_squared3, DTYPE fs_squared4, DTYPE fs_over_Q1, DTYPE fs_over_Q2, DTYPE fs_over_Q3, DTYPE fs_over_Q4, DTYPE f1, DTYPE f2, DTYPE complex Cinv_tdep1, complex DTYPE tdep_sensing_at_f2, DTYPE lowest_error, char *debug_file) { \
* Any of the four solutions can be "correct". The variable gain of the sensing
* function is kappa_C / f_s^2, and there are 3 time-dependent poles, which depend
* on f_cc, f_s, and Q. The algorithm does not care which variable name is
* assigned to which pole or how the gain is distributed, but we need to care
* because we associate a physical meaning with each variable. We require that:
* 1. kappa_C > 0
* 2. f_cc is the highest-frequency pole of the sensing function.
Aaron Viets
committed
* Of the remaining solutions, we choose the one with the smallest sum of square
* errors at the two Pcal line frequencies. Errors can be caused by numerical
* instabilities.
Aaron Viets
committed
DTYPE error1, error2, error3, error4; \
error1 = 2.0; \
error2 = 2.0; \
error3 = 2.0; \
error4 = 2.0; \
complex DTYPE fs1, fs2, fs3, fs4, Qinv1, Qinv2, Qinv3, Qinv4; \
guint best_solution = 0; \
/* Condition 1 */ \
if(kc1 > 0) { \
fs1 = cpow ## F_OR_NOT((complex DTYPE) fs_squared1, 0.5); \
Qinv1 = fs_over_Q1 / fs1; \
/* Condition 2 */ \
if(fcc1 > cabs ## F_OR_NOT(fs1 / 2 * (Qinv1 + cpow ## F_OR_NOT(cpow ## F_OR_NOT(Qinv1, 2) + 4, 0.5))) && fcc1 > cabs ## F_OR_NOT(fs1 / 2 * (Qinv1 - cpow ## F_OR_NOT(cpow ## F_OR_NOT(Qinv1, 2) + 4, 0.5)))) { \
error1 = cpow ## F_OR_NOT(cabs ## F_OR_NOT(((1 + I * f1 / fcc1) / kc1) * ((f1 * f1 + fs_squared1 - I * f1 * fs_over_Q1) / (f1 * f1)) / Cinv_tdep1 - 1), 2); \
error1 += cpow ## F_OR_NOT(cabs ## F_OR_NOT(((1 + I * f2 / fcc1) / kc1) * ((f2 * f2 + fs_squared1 - I * f2 * fs_over_Q1) / (f2 * f2)) / tdep_sensing_at_f2 - 1), 2); \
if(error1 < lowest_error) { \
lowest_error = error1; \
best_solution = 1; \
} \
} \
} \
if(kc2 > 0) { \
fs2 = cpow ## F_OR_NOT((complex DTYPE) fs_squared2, 0.5); \
Qinv2 = fs_over_Q2 / fs2; \
if(fcc2 > cabs ## F_OR_NOT(fs2 / 2 * (Qinv2 + cpow ## F_OR_NOT(cpow ## F_OR_NOT(Qinv2, 2) + 4, 0.5))) && fcc2 > cabs ## F_OR_NOT(fs2 / 2 * (Qinv2 - cpow ## F_OR_NOT(cpow ## F_OR_NOT(Qinv2, 2) + 4, 0.5)))) { \
error2 = cpow ## F_OR_NOT(cabs ## F_OR_NOT(((1 + I * f1 / fcc2) / kc2) * ((f1 * f1 + fs_squared2 - I * f1 * fs_over_Q2) / (f1 * f1)) / Cinv_tdep1 - 1), 2); \
error2 += cpow ## F_OR_NOT(cabs ## F_OR_NOT(((1 + I * f2 / fcc2) / kc2) * ((f2 * f2 + fs_squared2 - I * f2 * fs_over_Q2) / (f2 * f2)) / tdep_sensing_at_f2 - 1), 2); \
if(error2 < lowest_error) { \
lowest_error = error2; \
best_solution = 2; \
} \
} \
} \
if(kc3 > 0) { \
fs3 = cpow ## F_OR_NOT((complex DTYPE) fs_squared3, 0.5); \
Qinv3 = fs_over_Q3 / fs3; \
if(fcc3 > cabs ## F_OR_NOT(fs3 / 2 * (Qinv3 + cpow ## F_OR_NOT(cpow ## F_OR_NOT(Qinv3, 2) + 4, 0.5))) && fcc3 > cabs ## F_OR_NOT(fs3 / 2 * (Qinv3 - cpow ## F_OR_NOT(cpow ## F_OR_NOT(Qinv3, 2) + 4, 0.5)))) { \
error3 = cpow ## F_OR_NOT(cabs ## F_OR_NOT(((1 + I * f1 / fcc3) / kc3) * ((f1 * f1 + fs_squared3 - I * f1 * fs_over_Q3) / (f1 * f1)) / Cinv_tdep1 - 1), 2); \
error3 += cpow ## F_OR_NOT(cabs ## F_OR_NOT(((1 + I * f2 / fcc3) / kc3) * ((f2 * f2 + fs_squared3 - I * f2 * fs_over_Q3) / (f2 * f2)) / tdep_sensing_at_f2 - 1), 2); \
if(error3 < lowest_error) { \
lowest_error = error3; \
best_solution = 3; \
} \
} \
} \
if(kc4 > 0) { \
fs4 = cpow ## F_OR_NOT((complex DTYPE) fs_squared4, 0.5); \
Qinv4 = fs_over_Q4 / fs4; \
if(fcc4 > cabs ## F_OR_NOT(fs4 / 2 * (Qinv4 + cpow ## F_OR_NOT(cpow ## F_OR_NOT(Qinv4, 2) + 4, 0.5))) && fcc4 > cabs ## F_OR_NOT(fs4 / 2 * (Qinv4 - cpow ## F_OR_NOT(cpow ## F_OR_NOT(Qinv4, 2) + 4, 0.5)))) { \
error4 = cpow ## F_OR_NOT(cabs ## F_OR_NOT(((1 + I * f1 / fcc4) / kc4) * ((f1 * f1 + fs_squared4 - I * f1 * fs_over_Q4) / (f1 * f1)) / Cinv_tdep1 - 1), 2); \
error4 += cpow ## F_OR_NOT(cabs ## F_OR_NOT(((1 + I * f2 / fcc4) / kc4) * ((f2 * f2 + fs_squared4 - I * f2 * fs_over_Q4) / (f2 * f2)) / tdep_sensing_at_f2 - 1), 2); \
if(error4 < lowest_error) { \
lowest_error = error4; \
best_solution = 4; \
} \
} \
} \
Aaron Viets
committed
if(debug_file) { \
FILE *fp; \
fp = fopen(debug_file, "a"); \
g_fprintf(fp, "\nFinalized Solutions\n"); \
g_fprintf(fp, "kc1=%e fcc1=%e fs_squared1=%e fs_over_Q1=%e error1=%e\n", kc1, fcc1, fs_squared1, fs_over_Q1, error1); \
g_fprintf(fp, "kc2=%e fcc2=%e fs_squared2=%e fs_over_Q2=%e error2=%e\n", kc2, fcc2, fs_squared2, fs_over_Q2, error2); \
g_fprintf(fp, "kc3=%e fcc3=%e fs_squared3=%e fs_over_Q3=%e error3=%e\n", kc3, fcc3, fs_squared3, fs_over_Q3, error3); \
g_fprintf(fp, "kc4=%e fcc4=%e fs_squared4=%e fs_over_Q4=%e error4=%e\n", kc4, fcc4, fs_squared4, fs_over_Q4, error4); \
g_fprintf(fp, "best_solution=%d lowest_error=%e\n\n\n", best_solution, lowest_error); \
fclose(fp); \
} \
return best_solution; \
}
DEFINE_FIND_BEST_SOLUTION(float, f);
DEFINE_FIND_BEST_SOLUTION(double, );
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
/*
* ============================================================================
*
* GstBaseTransform Method Overrides
*
* ============================================================================
*/
/*
* get_unit_size()
*/
static gboolean get_unit_size(GstBaseTransform *trans, GstCaps *caps, gsize *size)
{
GstAudioInfo info;
gboolean success = TRUE;
success &= gstlal_audio_info_from_caps(&info, caps);
if(success) {
*size = GST_AUDIO_INFO_BPF(&info);
} else
GST_WARNING_OBJECT(trans, "unable to parse caps %" GST_PTR_FORMAT, caps);
return success;
}
/*
* transform_caps()
*/
static GstCaps *transform_caps(GstBaseTransform *trans, GstPadDirection direction, GstCaps *caps, GstCaps *filter) {
guint n;
caps = gst_caps_normalize(gst_caps_copy(caps));
switch(direction) {
case GST_PAD_SRC:
/*
* We know the caps on the source pad, and we want to put constraints on
* the sink pad caps. The sink pad caps are complex, while the source pad
* caps are real with twice as many channels.
*/
for(n = 0; n < gst_caps_get_size(caps); n++) {
GstStructure *str = gst_caps_get_structure(caps, n);
const GValue *v = gst_structure_get_value(str, "channels");
if(GST_VALUE_HOLDS_INT_RANGE(v)) {
gint channels_out_min, channels_out_max;
channels_out_min = gst_value_get_int_range_min(v);
channels_out_max = gst_value_get_int_range_max(v);
g_assert_cmpint(channels_out_max, >=, 2);
gst_structure_set(str, "channels", GST_TYPE_INT_RANGE, channels_out_min / 2, channels_out_max / 2, NULL);
} else if(G_VALUE_HOLDS_INT(v)) {
gint channels_out = g_value_get_int(v);
/* The number of output channels must be even. */
g_assert(channels_out % 2 == 0);
gst_structure_set(str, "channels", G_TYPE_INT, channels_out / 2, NULL);
} else
GST_ELEMENT_ERROR(trans, CORE, NEGOTIATION, (NULL), ("invalid type for channels in caps"));
const gchar *format = gst_structure_get_string(str, "format");
if(!format) {
GST_DEBUG_OBJECT(trans, "unrecognized caps %" GST_PTR_FORMAT, caps);
goto error;
} else if(!strcmp(format, GST_AUDIO_NE(F32)))
gst_structure_set(str, "format", G_TYPE_STRING, GST_AUDIO_NE(Z64), NULL);
else if(!strcmp(format, GST_AUDIO_NE(F64)))
gst_structure_set(str, "format", G_TYPE_STRING, GST_AUDIO_NE(Z128), NULL);
else {
GST_DEBUG_OBJECT(trans, "unrecognized format %s in %" GST_PTR_FORMAT, format, caps);
goto error;
}
}
break;
case GST_PAD_SINK:
/*
* We know the caps on the sink pad, and we want to put constraints on
* the source pad caps. The sink pad caps are complex, while the source pad
* caps are real with twice as many channels.
*/
for(n = 0; n < gst_caps_get_size(caps); n++) {
GstStructure *str = gst_caps_get_structure(caps, n);
const GValue *v = gst_structure_get_value(str, "channels");
if(GST_VALUE_HOLDS_INT_RANGE(v)) {
gint channels_in_min, channels_in_max;
channels_in_min = gst_value_get_int_range_min(v);
channels_in_max = gst_value_get_int_range_max(v);
gst_structure_set(str, "channels", GST_TYPE_INT_RANGE, 2 * channels_in_min, 2 * channels_in_max, NULL);
} else if(G_VALUE_HOLDS_INT(v)) {
gint channels_in;
channels_in = g_value_get_int(v);
gst_structure_set(str, "channels", G_TYPE_INT, 2 * channels_in, NULL);
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
} else
GST_ELEMENT_ERROR(trans, CORE, NEGOTIATION, (NULL), ("invalid type for channels in caps"));
const gchar *format = gst_structure_get_string(str, "format");
if(!format) {
GST_DEBUG_OBJECT(trans, "unrecognized caps %" GST_PTR_FORMAT, caps);
goto error;
} else if(!strcmp(format, GST_AUDIO_NE(Z64)))
gst_structure_set(str, "format", G_TYPE_STRING, GST_AUDIO_NE(F32), NULL);
else if(!strcmp(format, GST_AUDIO_NE(Z128)))
gst_structure_set(str, "format", G_TYPE_STRING, GST_AUDIO_NE(F64), NULL);
else {
GST_DEBUG_OBJECT(trans, "unrecognized format %s in %" GST_PTR_FORMAT, format, caps);
goto error;
}
}
break;
case GST_PAD_UNKNOWN:
GST_ELEMENT_ERROR(trans, CORE, NEGOTIATION, (NULL), ("invalid direction GST_PAD_UNKNOWN"));
goto error;
default:
g_assert_not_reached();
}
if(filter) {
GstCaps *intersection = gst_caps_intersect(caps, filter);
gst_caps_unref(caps);
caps = intersection;
}
return gst_caps_simplify(caps);
error:
gst_caps_unref(caps);
return GST_CAPS_NONE;
}
/*
* set_caps()
*/
static gboolean set_caps(GstBaseTransform *trans, GstCaps *incaps, GstCaps *outcaps)
{
GSTLALSensingTDCFs *element = GSTLAL_SENSINGTDCFS(trans);
gint rate, channels_in, channels_out;
gsize unit_size_in, unit_size_out;
/*
* parse the caps
*/
GstStructure *outstr = gst_caps_get_structure(outcaps, 0);
GstStructure *instr = gst_caps_get_structure(incaps, 0);
const gchar *name = gst_structure_get_string(outstr, "format");
if(!name) {
GST_DEBUG_OBJECT(element, "unable to parse format from %" GST_PTR_FORMAT, outcaps);
return FALSE;
}
if(!get_unit_size(trans, outcaps, &unit_size_out)) {
GST_DEBUG_OBJECT(element, "function 'get_unit_size' failed");
return FALSE;
}
if(!get_unit_size(trans, incaps, &unit_size_in)) {
GST_DEBUG_OBJECT(element, "function 'get_unit_size' failed");
return FALSE;
}
if(!gst_structure_get_int(outstr, "rate", &rate)) {
GST_DEBUG_OBJECT(element, "unable to parse rate from %" GST_PTR_FORMAT, outcaps);
return FALSE;
}
if(!gst_structure_get_int(outstr, "channels", &channels_out)) {
GST_DEBUG_OBJECT(element, "unable to parse channels from %" GST_PTR_FORMAT, outcaps);
return FALSE;
}
if(!gst_structure_get_int(instr, "channels", &channels_in)) {
GST_DEBUG_OBJECT(element, "unable to parse channels from %" GST_PTR_FORMAT, incaps);
return FALSE;
}
/* Requirements for channels */
g_assert(channels_out == 2 * channels_in);
/*
* record stream parameters
*/
if(!strcmp(name, GST_AUDIO_NE(F32))) {
element->data_type = GSTLAL_SENSINGTDCFS_FLOAT;
g_assert_cmpuint(unit_size_out, ==, 4 * (guint) channels_out);
Aaron Viets
committed
g_assert_cmpuint(unit_size_in, ==, 8 * (guint) channels_in);
} else if(!strcmp(name, GST_AUDIO_NE(F64))) {
element->data_type = GSTLAL_SENSINGTDCFS_DOUBLE;
g_assert_cmpuint(unit_size_out, ==, 8 * (guint) channels_out);
Aaron Viets
committed
g_assert_cmpuint(unit_size_in, ==,16 * (guint) channels_in);
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
} else
g_assert_not_reached();
element->rate = rate;
element->channels_out = channels_out;
element->channels_in = channels_in;
element->unit_size_out = unit_size_out;
element->unit_size_in = unit_size_in;
return TRUE;
}
/*
* start()
*/
static gboolean start(GstBaseTransform *trans) {
GSTLALSensingTDCFs *element = GSTLAL_SENSINGTDCFS(trans);
element->t0 = GST_CLOCK_TIME_NONE;
element->offset0 = GST_BUFFER_OFFSET_NONE;
element->next_in_offset = GST_BUFFER_OFFSET_NONE;
element->next_out_offset = GST_BUFFER_OFFSET_NONE;
element->need_discont = TRUE;
/* Sanity checks */
if(element->freq1 == G_MAXDOUBLE)
GST_WARNING_OBJECT(element, "freq1 was not set. It must be set in order to produce sensible output.");
if(element->freq2 == G_MAXDOUBLE)
GST_WARNING_OBJECT(element, "freq2 was not set. It must be set in order to produce sensible output.");
if(element->sensing_model == 2 && element->freq2 == G_MAXDOUBLE)
GST_WARNING_OBJECT(element, "freq4 was not set. When using sensing-model is 2, freq4 must be set in order to produce sensible output.");
Aaron Viets
committed
/* If we are writing debugging output to file, and a file already exists with the same name, remove it */
if(element->debug_file)
remove(element->debug_file);
element->consecutive_failures = 0;
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
return TRUE;
}
/*
* transform()
*/
static GstFlowReturn transform(GstBaseTransform *trans, GstBuffer *inbuf, GstBuffer *outbuf) {
GSTLALSensingTDCFs *element = GSTLAL_SENSINGTDCFS(trans);
GstMapInfo inmap, outmap;
GstFlowReturn result = GST_FLOW_OK;
/*
* check for discontinuity
*/
if(G_UNLIKELY(GST_BUFFER_IS_DISCONT(inbuf) || GST_BUFFER_OFFSET(inbuf) != element->next_in_offset || !GST_CLOCK_TIME_IS_VALID(element->t0))) {
GST_DEBUG_OBJECT(element, "pushing discontinuous buffer at input timestamp %lu", (long unsigned) GST_TIME_AS_SECONDS(GST_BUFFER_PTS(inbuf)));
element->t0 = GST_BUFFER_PTS(inbuf);
element->offset0 = element->next_out_offset = GST_BUFFER_OFFSET(inbuf);
element->need_discont = TRUE;
}
element->next_in_offset = GST_BUFFER_OFFSET_END(inbuf);
/*
* process buffer
*/
if(!GST_BUFFER_FLAG_IS_SET(inbuf, GST_BUFFER_FLAG_GAP)) {
/*
* input is not gap.
*/
gst_buffer_map(inbuf, &inmap, GST_MAP_READ);
gst_buffer_map(outbuf, &outmap, GST_MAP_WRITE);
if(element->data_type == GSTLAL_SENSINGTDCFS_FLOAT) {
complex float *indata = (complex float *) inmap.data;
float *outdata = (float *) outmap.data;
guint64 i, samples = outmap.size / element->unit_size_out;
float f1 = (float) element->freq1;
float f2 = (float) element->freq2;
float kappa_C1, kappa_C2, kappa_C3, kappa_C4, f_cc1, f_cc2, f_cc3, f_cc4, f_s_squared1, f_s_squared2, f_s_squared3, f_s_squared4, f_s_over_Q1, f_s_over_Q2, f_s_over_Q3, f_s_over_Q4;
guint best_solution;
/*
* kappa_C is the solution of a quartic equation, so there are 4 solutions
* for kappa_C. For each solution, a different value for fcc, fs, and Q
* will be computed, and the best solution will be determined at the end.
*/
kappa_C1 = kappa_C2 = kappa_C3 = kappa_C4 = f_cc1 = f_cc2 = f_cc3 = f_cc4 = 0.0;
switch(element->sensing_model) {
case 0: ;
f_s_squared1 = f_s_squared2 = f_s_squared3 = f_s_squared4 = f_s_over_Q1 = f_s_over_Q2 = f_s_over_Q3 = f_s_over_Q4 = 0.0;
for(i = 0; i < samples; i++) {
Aaron Viets
committed
kappa_C_0_float(crealf(indata[2 * i]), cimagf(indata[2 * i]), crealf(indata[2 * i + 1]), cimagf(indata[2 * i + 1]), f1, f2, &kappa_C1, &kappa_C2, &kappa_C3, &kappa_C4, element->debug_file);
/* Only compute f_cc, f_s, and Q if we have to. */
if(kappa_C1 != 0.0) {
f_cc1 = f_cc_0_float(cimagf(indata[2 * i]), cimagf(indata[2 * i + 1]), f1, f2, kappa_C1);
f_s_squared1 = f_s_squared_0_float(crealf(indata[2 * i]), crealf(indata[2 * i + 1]), f1, f2, kappa_C1);
f_s_over_Q1 = f_s_over_Q_0_float(crealf(indata[2 * i]), f1, kappa_C1, f_cc1, f_s_squared1);
}
if(kappa_C2 != 0.0) {
f_cc2 = f_cc_0_float(cimagf(indata[2 * i]), cimagf(indata[2 * i + 1]), f1, f2, kappa_C2);
f_s_squared2 = f_s_squared_0_float(crealf(indata[2 * i]), crealf(indata[2 * i + 1]), f1, f2, kappa_C2);
f_s_over_Q2 = f_s_over_Q_0_float(crealf(indata[2 * i]), f1, kappa_C2, f_cc2, f_s_squared2);
}
if(kappa_C3 != 0.0) {
f_cc3 = f_cc_0_float(cimagf(indata[2 * i]), cimagf(indata[2 * i + 1]), f1, f2, kappa_C3);
f_s_squared3 = f_s_squared_0_float(crealf(indata[2 * i]), crealf(indata[2 * i + 1]), f1, f2, kappa_C3);
f_s_over_Q3 = f_s_over_Q_0_float(crealf(indata[2 * i]), f1, kappa_C3, f_cc3, f_s_squared3);
}
if(kappa_C4 != 0.0) {
f_cc4 = f_cc_0_float(cimagf(indata[2 * i]), cimagf(indata[2 * i + 1]), f1, f2, kappa_C4);
f_s_squared4 = f_s_squared_0_float(crealf(indata[2 * i]), crealf(indata[2 * i + 1]), f1, f2, kappa_C4);
f_s_over_Q4 = f_s_over_Q_0_float(crealf(indata[2 * i]), f1, kappa_C4, f_cc4, f_s_squared4);
/* Determine which solution is our favorite. */
Aaron Viets
committed
best_solution = find_best_solution_float(kappa_C1, kappa_C2, kappa_C3, kappa_C4, f_cc1, f_cc2, f_cc3, f_cc4, f_s_squared1, f_s_squared2, f_s_squared3, f_s_squared4, f_s_over_Q1, f_s_over_Q2, f_s_over_Q3, f_s_over_Q4, f1, f2, indata[2 * i], indata[2 * i + 1], (float) element->max_error, element->debug_file);
switch(best_solution) {
case 1:
outdata[4 * i] = element->current_kc = kappa_C1;
outdata[4 * i + 1] = element->current_fcc = f_cc1;
outdata[4 * i + 2] = element->current_fs_squared = f_s_squared1;
outdata[4 * i + 3] = element->current_fs_over_Q = f_s_over_Q1;
Aaron Viets
committed
element->consecutive_failures = 0;
outdata[4 * i] = element->current_kc = kappa_C2;
outdata[4 * i + 1] = element->current_fcc = f_cc2;
outdata[4 * i + 2] = element->current_fs_squared = f_s_squared2;
outdata[4 * i + 3] = element->current_fs_over_Q = f_s_over_Q2;
Aaron Viets
committed
element->consecutive_failures = 0;
outdata[4 * i] = element->current_kc = kappa_C3;
outdata[4 * i + 1] = element->current_fcc = f_cc3;
outdata[4 * i + 2] = element->current_fs_squared = f_s_squared3;
outdata[4 * i + 3] = element->current_fs_over_Q = f_s_over_Q3;
Aaron Viets
committed
element->consecutive_failures = 0;
outdata[4 * i] = element->current_kc = kappa_C4;
outdata[4 * i + 1] = element->current_fcc = f_cc4;
outdata[4 * i + 2] = element->current_fs_squared = f_s_squared4;
outdata[4 * i + 3] = element->current_fs_over_Q = f_s_over_Q4;
Aaron Viets
committed
element->consecutive_failures = 0;
Aaron Viets
committed
element->consecutive_failures++;
if(!(element->consecutive_failures % 1000))
GST_WARNING_OBJECT(element, "%lu: Unable to find a stable solution for the sensing function TDCFs for %d samples", GST_BUFFER_PTS(inbuf) / GST_SECOND, element->consecutive_failures);
Aaron Viets
committed
else
GST_DEBUG_OBJECT(element, "%lu: Unable to find a stable solution for the sensing function TDCFs", GST_BUFFER_PTS(inbuf) / GST_SECOND);
outdata[4 * i] = (float) element->current_kc;
outdata[4 * i + 1] = (float) element->current_fcc;
outdata[4 * i + 2] = (float) element->current_fs_squared;
outdata[4 * i + 3] = (float) element->current_fs_over_Q;
}
break;
case 1:
/* Solution not developed */
break;
case 2:
/* Solution not developed */
break;
default:
g_assert_not_reached();
}
} else if (element->data_type == GSTLAL_SENSINGTDCFS_DOUBLE) {
complex double *indata = (complex double *) inmap.data;
double *outdata = (double *) outmap.data;
guint64 i, samples = outmap.size / element->unit_size_out;
double f1 = element->freq1;
double f2 = element->freq2;
double kappa_C1, kappa_C2, kappa_C3, kappa_C4, f_cc1, f_cc2, f_cc3, f_cc4, f_s_squared1, f_s_squared2, f_s_squared3, f_s_squared4, f_s_over_Q1, f_s_over_Q2, f_s_over_Q3, f_s_over_Q4;
guint best_solution;
/*
* kappa_C is the solution of a quartic equation, so there are 4 solutions
* for kappa_C. For each solution, a different value for fcc, fs, and Q
* will be computed, and the best solution will be determined at the end.
*/
kappa_C1 = kappa_C2 = kappa_C3 = kappa_C4 = f_cc1 = f_cc2 = f_cc3 = f_cc4 = 0.0;
switch(element->sensing_model) {
case 0: ;
f_s_squared1 = f_s_squared2 = f_s_squared3 = f_s_squared4 = f_s_over_Q1 = f_s_over_Q2 = f_s_over_Q3 = f_s_over_Q4 = 0.0;
for(i = 0; i < samples; i++) {
Aaron Viets
committed
kappa_C_0_double(crealf(indata[2 * i]), cimagf(indata[2 * i]), crealf(indata[2 * i + 1]), cimagf(indata[2 * i + 1]), f1, f2, &kappa_C1, &kappa_C2, &kappa_C3, &kappa_C4, element->debug_file);
/* Only compute f_cc, f_s, and Q if we have to. */
if(kappa_C1 != 0.0) {
f_cc1 = f_cc_0_double(cimagf(indata[2 * i]), cimagf(indata[2 * i + 1]), f1, f2, kappa_C1);
f_s_squared1 = f_s_squared_0_double(crealf(indata[2 * i]), crealf(indata[2 * i + 1]), f1, f2, kappa_C1);
f_s_over_Q1 = f_s_over_Q_0_double(crealf(indata[2 * i]), f1, kappa_C1, f_cc1, f_s_squared1);
}
if(kappa_C2 != 0.0) {
f_cc2 = f_cc_0_double(cimagf(indata[2 * i]), cimagf(indata[2 * i + 1]), f1, f2, kappa_C2);
f_s_squared2 = f_s_squared_0_double(crealf(indata[2 * i]), crealf(indata[2 * i + 1]), f1, f2, kappa_C2);
f_s_over_Q2 = f_s_over_Q_0_double(crealf(indata[2 * i]), f1, kappa_C2, f_cc2, f_s_squared2);
}
if(kappa_C3 != 0.0) {
f_cc3 = f_cc_0_double(cimagf(indata[2 * i]), cimagf(indata[2 * i + 1]), f1, f2, kappa_C3);
f_s_squared3 = f_s_squared_0_double(crealf(indata[2 * i]), crealf(indata[2 * i + 1]), f1, f2, kappa_C3);
f_s_over_Q3 = f_s_over_Q_0_double(crealf(indata[2 * i]), f1, kappa_C3, f_cc3, f_s_squared3);
}
if(kappa_C4 != 0.0) {
f_cc4 = f_cc_0_double(cimagf(indata[2 * i]), cimagf(indata[2 * i + 1]), f1, f2, kappa_C4);
f_s_squared4 = f_s_squared_0_double(crealf(indata[2 * i]), crealf(indata[2 * i + 1]), f1, f2, kappa_C4);
f_s_over_Q4 = f_s_over_Q_0_double(crealf(indata[2 * i]), f1, kappa_C4, f_cc4, f_s_squared4);
/* Determine which solution is our favorite. */
Aaron Viets
committed
best_solution = find_best_solution_double(kappa_C1, kappa_C2, kappa_C3, kappa_C4, f_cc1, f_cc2, f_cc3, f_cc4, f_s_squared1, f_s_squared2, f_s_squared3, f_s_squared4, f_s_over_Q1, f_s_over_Q2, f_s_over_Q3, f_s_over_Q4, f1, f2, indata[2 * i], indata[2 * i + 1], element->max_error, element->debug_file);
switch(best_solution) {
case 1:
outdata[4 * i] = element->current_kc = kappa_C1;
outdata[4 * i + 1] = element->current_fcc = f_cc1;
outdata[4 * i + 2] = element->current_fs_squared = f_s_squared1;
outdata[4 * i + 3] = element->current_fs_over_Q = f_s_over_Q1;
Aaron Viets
committed
element->consecutive_failures = 0;
outdata[4 * i] = element->current_kc = kappa_C2;
outdata[4 * i + 1] = element->current_fcc = f_cc2;
outdata[4 * i + 2] = element->current_fs_squared = f_s_squared2;
outdata[4 * i + 3] = element->current_fs_over_Q = f_s_over_Q2;
Aaron Viets
committed
element->consecutive_failures = 0;
outdata[4 * i] = element->current_kc = kappa_C3;
outdata[4 * i + 1] = element->current_fcc = f_cc3;
outdata[4 * i + 2] = element->current_fs_squared = f_s_squared3;
outdata[4 * i + 3] = element->current_fs_over_Q = f_s_over_Q3;
Aaron Viets
committed
element->consecutive_failures = 0;
outdata[4 * i] = element->current_kc = kappa_C4;
outdata[4 * i + 1] = element->current_fcc = f_cc4;
outdata[4 * i + 2] = element->current_fs_squared = f_s_squared4;
outdata[4 * i + 3] = element->current_fs_over_Q = f_s_over_Q4;
Aaron Viets
committed
element->consecutive_failures = 0;
Aaron Viets
committed
element->consecutive_failures++;
if(!(element->consecutive_failures % 1000))
GST_WARNING_OBJECT(element, "%lu: Unable to find a stable solution for the sensing function TDCFs for %d samples", GST_BUFFER_PTS(inbuf) / GST_SECOND, element->consecutive_failures);
Aaron Viets
committed
else
GST_DEBUG_OBJECT(element, "%lu: Unable to find a stable solution for the sensing function TDCFs", GST_BUFFER_PTS(inbuf) / GST_SECOND);
outdata[4 * i] = element->current_kc;
outdata[4 * i + 1] = element->current_fcc;
outdata[4 * i + 2] = element->current_fs_squared;
outdata[4 * i + 3] = element->current_fs_over_Q;
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
}
break;
case 1:
/* Solution not developed */
break;
case 2:
/* Solution not developed */
break;
default:
g_assert_not_reached();
}
} else {
g_assert_not_reached();
}
set_metadata(element, outbuf, outmap.size / element->unit_size_out, FALSE);
gst_buffer_unmap(outbuf, &outmap);
gst_buffer_unmap(inbuf, &inmap);
} else {
/*
* input is gap.
*/
gst_buffer_map(outbuf, &outmap, GST_MAP_WRITE);
memset(outmap.data, 0, outmap.size);
set_metadata(element, outbuf, outmap.size / element->unit_size_out, TRUE);
gst_buffer_unmap(outbuf, &outmap);
}
/*
* done
*/
return result;
}
/*
* ============================================================================
*
* GObject Method Overrides
*
* ============================================================================
*/
/*
* properties
*/
enum property {
ARG_SENSING_MODEL = 1,
ARG_FREQ1,
ARG_FREQ2,
ARG_CURRENT_KC,
ARG_CURRENT_FCC,
ARG_CURRENT_FS_SQUARED,
Aaron Viets
committed
ARG_CURRENT_FS_OVER_Q,
ARG_MAX_ERROR,
ARG_DEBUG_FILE
};
static void set_property(GObject *object, enum property prop_id, const GValue *value, GParamSpec *pspec) {
GSTLALSensingTDCFs *element = GSTLAL_SENSINGTDCFS(object);
GST_OBJECT_LOCK(element);
switch (prop_id) {
case ARG_SENSING_MODEL:
element->sensing_model = g_value_get_int(value);
break;
case ARG_FREQ1:
element->freq1 = g_value_get_double(value);
break;
case ARG_FREQ2:
element->freq2 = g_value_get_double(value);
break;
case ARG_FREQ4:
element->freq4 = g_value_get_double(value);
break;
case ARG_CURRENT_KC:
element->current_kc = g_value_get_double(value);
case ARG_CURRENT_FCC:
element->current_fcc = g_value_get_double(value);
case ARG_CURRENT_FS_SQUARED:
element->current_fs_squared = g_value_get_double(value);
break;
case ARG_CURRENT_FS_OVER_Q:
element->current_fs_over_Q = g_value_get_double(value);
Aaron Viets
committed
case ARG_MAX_ERROR:
element->max_error = g_value_get_double(value);
break;
case ARG_DEBUG_FILE:
element->debug_file = g_value_dup_string(value);
break;
default:
G_OBJECT_WARN_INVALID_PROPERTY_ID(object, prop_id, pspec);
break;
}
GST_OBJECT_UNLOCK(element);
}
static void get_property(GObject *object, enum property prop_id, GValue *value, GParamSpec *pspec) {
GSTLALSensingTDCFs *element = GSTLAL_SENSINGTDCFS(object);
GST_OBJECT_LOCK(element);
switch (prop_id) {
case ARG_SENSING_MODEL:
g_value_set_int(value, element->sensing_model);
break;
case ARG_FREQ1:
g_value_set_double(value, element->freq1);