projects
/
ctsim.git
/ blobdiff
commit
grep
author
committer
pickaxe
?
search:
re
summary
|
shortlog
|
log
|
commit
|
commitdiff
|
tree
raw
|
inline
| side by side
r131: *** empty log message ***
[ctsim.git]
/
libctsim
/
filter.cpp
diff --git
a/libctsim/filter.cpp
b/libctsim/filter.cpp
index 27d1829380ae5224feff90fb221a801079884f20..d3a3fbdfaf01ed399ba8e2923a55ee0fd75e3443 100644
(file)
--- a/
libctsim/filter.cpp
+++ b/
libctsim/filter.cpp
@@
-9,7
+9,7
@@
** This is part of the CTSim program
** Copyright (C) 1983-2000 Kevin Rosenberg
**
** This is part of the CTSim program
** Copyright (C) 1983-2000 Kevin Rosenberg
**
-** $Id: filter.cpp,v 1.
8 2000/07/04 18:33:35
kevin Exp $
+** $Id: filter.cpp,v 1.
9 2000/07/04 22:21:01
kevin Exp $
**
** This program is free software; you can redistribute it and/or modify
** it under the terms of the GNU General Public License (version 2) as
**
** This program is free software; you can redistribute it and/or modify
** it under the terms of the GNU General Public License (version 2) as
@@
-115,6
+115,7
@@
SignalFilter::init (const FilterID filterID, const FilterMethodID filterMethodID
m_fail = true;
return;
}
m_fail = true;
return;
}
+ m_traceLevel = TRACE_NONE;
m_nameFilter = convertFilterIDToName (m_idFilter);
m_nameDomain = convertDomainIDToName (m_idDomain);
m_nameFilterMethod = convertFilterMethodIDToName (m_idFilterMethod);
m_nameFilter = convertFilterIDToName (m_idFilter);
m_nameDomain = convertDomainIDToName (m_idDomain);
m_nameFilterMethod = convertFilterMethodIDToName (m_idFilterMethod);
@@
-133,9
+134,9
@@
SignalFilter::init (const FilterID filterID, const FilterMethodID filterMethodID
m_vecFourierSinTable[i] = sin (angleIncrement * i);
}
m_nFilterPoints = m_nSignalPoints;
m_vecFourierSinTable[i] = sin (angleIncrement * i);
}
m_nFilterPoints = m_nSignalPoints;
- m_filterMin =
0
;
- m_filterMax =
m_nSignalPoints * m_signalInc
;
- m_filterInc = (m_filterMax - m_filterMin) /
(m_nFilterPoints - 1)
;
+ m_filterMin =
-1. / (2 * m_signalInc)
;
+ m_filterMax =
1. / (2 * m_signalInc)
;
+ m_filterInc = (m_filterMax - m_filterMin) /
m_nFilterPoints
;
m_vecFilter = new double [m_nFilterPoints];
int halfFilter = m_nFilterPoints / 2;
for (int i = 0; i < halfFilter; i++)
m_vecFilter = new double [m_nFilterPoints];
int halfFilter = m_nFilterPoints / 2;
for (int i = 0; i < halfFilter; i++)
@@
-143,7
+144,7
@@
SignalFilter::init (const FilterID filterID, const FilterMethodID filterMethodID
for (int i = 0; i < halfFilter; i++)
m_vecFilter[m_nFilterPoints - i - 1] = static_cast<double>(i) / (halfFilter - 1) / (2 * m_signalInc);
if (halfFilter % 2) // odd
for (int i = 0; i < halfFilter; i++)
m_vecFilter[m_nFilterPoints - i - 1] = static_cast<double>(i) / (halfFilter - 1) / (2 * m_signalInc);
if (halfFilter % 2) // odd
- m_vecFilter[halfFilter] = 1;
+ m_vecFilter[halfFilter] = 1
/ (2 * m_signalInc)
;
} else if (m_idFilterMethod == FILTER_METHOD_FFT || m_idFilterMethod == FILTER_METHOD_FFT_ZEROPAD_2 || m_idFilterMethod == FILTER_METHOD_FFT_ZEROPAD_4) {
m_nFilterPoints = m_nSignalPoints;
if (m_idFilterMethod == FILTER_METHOD_FFT_ZEROPAD_2 || m_idFilterMethod == FILTER_METHOD_FFT_ZEROPAD_4) {
} else if (m_idFilterMethod == FILTER_METHOD_FFT || m_idFilterMethod == FILTER_METHOD_FFT_ZEROPAD_2 || m_idFilterMethod == FILTER_METHOD_FFT_ZEROPAD_4) {
m_nFilterPoints = m_nSignalPoints;
if (m_idFilterMethod == FILTER_METHOD_FFT_ZEROPAD_2 || m_idFilterMethod == FILTER_METHOD_FFT_ZEROPAD_4) {
@@
-156,9
+157,9
@@
SignalFilter::init (const FilterID filterID, const FilterMethodID filterMethodID
m_nFilterPoints = 1 << nextPowerOf2;
cout << "nFilterPoints = " << m_nFilterPoints << endl;
}
m_nFilterPoints = 1 << nextPowerOf2;
cout << "nFilterPoints = " << m_nFilterPoints << endl;
}
- m_filterMin =
0
;
- m_filterMax =
m_nSignalPoints * m_signalInc
;
- m_filterInc = (m_filterMax - m_filterMin) /
(m_nFilterPoints - 1)
;
+ m_filterMin =
-1. / (2 * m_signalInc)
;
+ m_filterMax =
1. / (2 * m_signalInc)
;
+ m_filterInc = (m_filterMax - m_filterMin) /
m_nFilterPoints
;
m_vecFilter = new double [m_nFilterPoints];
int halfFilter = m_nFilterPoints / 2;
for (int i = 0; i < halfFilter; i++)
m_vecFilter = new double [m_nFilterPoints];
int halfFilter = m_nFilterPoints / 2;
for (int i = 0; i < halfFilter; i++)
@@
-166,7
+167,7
@@
SignalFilter::init (const FilterID filterID, const FilterMethodID filterMethodID
for (int i = 0; i < halfFilter; i++)
m_vecFilter[m_nFilterPoints - i - 1] = static_cast<double>(i) / (halfFilter - 1) / (2 * m_signalInc);
if (halfFilter % 2) // odd
for (int i = 0; i < halfFilter; i++)
m_vecFilter[m_nFilterPoints - i - 1] = static_cast<double>(i) / (halfFilter - 1) / (2 * m_signalInc);
if (halfFilter % 2) // odd
- m_vecFilter[halfFilter] = 1;
+ m_vecFilter[halfFilter] = 1
/ (2 * m_signalInc)
;
#if HAVE_FFTW
m_planForward = fftw_create_plan (m_nFilterPoints, FFTW_FORWARD, FFTW_ESTIMATE);
#if HAVE_FFTW
m_planForward = fftw_create_plan (m_nFilterPoints, FFTW_FORWARD, FFTW_ESTIMATE);
@@
-364,11
+365,24
@@
SignalFilter::filterSignal (const float input[], double output[]) const
complex<double> complexOutput[m_nSignalPoints];
complex<double> filteredSignal[m_nSignalPoints];
finiteFourierTransform (input, fftSignal, m_nSignalPoints, -1);
complex<double> complexOutput[m_nSignalPoints];
complex<double> filteredSignal[m_nSignalPoints];
finiteFourierTransform (input, fftSignal, m_nSignalPoints, -1);
+ if (m_traceLevel >= TRACE_PLOT) {
+ double test[m_nSignalPoints];
+ for (int i = 0; i < m_nSignalPoints; i++)
+ test[i] = abs(fftSignal[i]);
+ ezplot_1d(test, m_nSignalPoints);
+ cio_kb_getc();
+ }
dotProduct (m_vecFilter, fftSignal, filteredSignal, m_nSignalPoints);
dotProduct (m_vecFilter, fftSignal, filteredSignal, m_nSignalPoints);
+ if (m_traceLevel >= TRACE_PLOT) {
+ double test[m_nSignalPoints];
+ for (int i = 0; i < m_nSignalPoints; i++)
+ test[i] = abs(filteredSignal[i]);
+ ezplot_1d(test, m_nSignalPoints);
+ cio_kb_getc();
+ }
finiteFourierTransform (filteredSignal, complexOutput, m_nSignalPoints, 1);
finiteFourierTransform (filteredSignal, complexOutput, m_nSignalPoints, 1);
- for (int i = 0; i < m_nSignalPoints; i++)
{
+ for (int i = 0; i < m_nSignalPoints; i++)
output[i] = abs( complexOutput[i] );
output[i] = abs( complexOutput[i] );
- }
} else if (m_idFilterMethod == FILTER_METHOD_FFT || FILTER_METHOD_FFT_ZEROPAD_2 || FILTER_METHOD_FFT_ZEROPAD_4) {
fftw_complex in[m_nFilterPoints], out[m_nFilterPoints];
for (int i = 0; i < m_nSignalPoints; i++) {
} else if (m_idFilterMethod == FILTER_METHOD_FFT || FILTER_METHOD_FFT_ZEROPAD_2 || FILTER_METHOD_FFT_ZEROPAD_4) {
fftw_complex in[m_nFilterPoints], out[m_nFilterPoints];
for (int i = 0; i < m_nSignalPoints; i++) {
@@
-379,11
+393,32
@@
SignalFilter::filterSignal (const float input[], double output[]) const
in[i].re = in[i].im = 0; // ZeroPad
}
fftw_one(m_planForward, in, out);
in[i].re = in[i].im = 0; // ZeroPad
}
fftw_one(m_planForward, in, out);
+ if (m_traceLevel >= TRACE_PLOT) {
+ double test[m_nFilterPoints];
+ for (int i = 0; i < m_nFilterPoints; i++)
+ test[i] = sqrt(out[i].re * out[i].re + out[i].im * out[i].im);
+ ezplot_1d(test, m_nFilterPoints);
+ cio_kb_getc();
+ }
for (int i = 0; i < m_nFilterPoints; i++) {
out[i].re = m_vecFilter[i] * out[i].re / m_nSignalPoints;
out[i].im = m_vecFilter[i] * out[i].im / m_nSignalPoints;
}
for (int i = 0; i < m_nFilterPoints; i++) {
out[i].re = m_vecFilter[i] * out[i].re / m_nSignalPoints;
out[i].im = m_vecFilter[i] * out[i].im / m_nSignalPoints;
}
+ if (m_traceLevel >= TRACE_PLOT) {
+ double test[m_nFilterPoints];
+ for (int i = 0; i < m_nFilterPoints; i++)
+ test[i] = sqrt(out[i].re * out[i].re + out[i].im * out[i].im);
+ ezplot_1d(test, m_nFilterPoints);
+ cio_kb_getc();
+ }
fftw_one(m_planBackward, out, in);
fftw_one(m_planBackward, out, in);
+ if (m_traceLevel >= TRACE_PLOT) {
+ double test[m_nFilterPoints];
+ for (int i = 0; i < m_nFilterPoints; i++)
+ test[i] = sqrt(in[i].re * in[i].re + in[i].im * in[i].im);
+ ezplot_1d(test, m_nFilterPoints);
+ cio_kb_getc();
+ }
for (int i = 0; i < m_nSignalPoints; i++)
output[i] = sqrt (in[i].re * in[i].re + in[i].im * in[i].im);
}
for (int i = 0; i < m_nSignalPoints; i++)
output[i] = sqrt (in[i].re * in[i].re + in[i].im * in[i].im);
}