← Back to branch summary

~valavanisalex/ubuntu/maverick/scidavis/fix-604811

« back to all changes in this revision

Viewing changes to scidavis/src/SmoothFilter.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Ruben Molina
  • Date: 2009-09-06 11:34:04 UTC
  • Revision ID: james.westby@ubuntu.com-20090906113404-4awaey82l3686w4q
Tags: upstream-0.2.3
Import upstream version 0.2.3

Show diffs side-by-side

added added

removed removed

Lines of Context:
scidavis/src/SmoothFilter.cpp
 
1
/***************************************************************************
 
2
    File                 : SmoothFilter.cpp
 
3
    Project              : SciDAVis
 
4
    --------------------------------------------------------------------
 
5
    Copyright            : (C) 2007 by Ion Vasilief
 
6
    Email (use @ for *)  : ion_vasilief*yahoo.fr
 
7
    Description          : Numerical smoothing of data sets
 
8
 
 
9
 ***************************************************************************/
 
10
 
 
11
/***************************************************************************
 
12
 *                                                                         *
 
13
 *  This program is free software; you can redistribute it and/or modify   *
 
14
 *  it under the terms of the GNU General Public License as published by   *
 
15
 *  the Free Software Foundation; either version 2 of the License, or      *
 
16
 *  (at your option) any later version.                                    *
 
17
 *                                                                         *
 
18
 *  This program is distributed in the hope that it will be useful,        *
 
19
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of         *
 
20
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          *
 
21
 *  GNU General Public License for more details.                           *
 
22
 *                                                                         *
 
23
 *   You should have received a copy of the GNU General Public License     *
 
24
 *   along with this program; if not, write to the Free Software           *
 
25
 *   Foundation, Inc., 51 Franklin Street, Fifth Floor,                    *
 
26
 *   Boston, MA  02110-1301  USA                                           *
 
27
 *                                                                         *
 
28
 ***************************************************************************/
 
29
#include "SmoothFilter.h"
 
30
#include "nrutil.h"
 
31
 
 
32
#include <QApplication>
 
33
#include <QMessageBox>
 
34
 
 
35
#include <gsl/gsl_fft_halfcomplex.h>
 
36
 
 
37
SmoothFilter::SmoothFilter(ApplicationWindow *parent, Graph *g, const QString& curveTitle, int m)
 
38
: Filter(parent, g)
 
39
{
 
40
        setDataFromCurve(curveTitle);
 
41
        init(m);
 
42
}
 
43
 
 
44
SmoothFilter::SmoothFilter(ApplicationWindow *parent, Graph *g, const QString& curveTitle,
 
45
                             double start, double end, int m)
 
46
: Filter(parent, g)
 
47
{
 
48
        setDataFromCurve(curveTitle, start, end);
 
49
    init(m);
 
50
}
 
51
 
 
52
void SmoothFilter::init (int m)
 
53
{
 
54
    setName(tr("Smoothed"));
 
55
    setMethod(m);
 
56
    d_points = d_n;
 
57
    d_smooth_points = 2;
 
58
    d_sav_gol_points = 2;
 
59
    d_polynom_order = 2;
 
60
}
 
61
 
 
62
 
 
63
void SmoothFilter::setMethod(int m)
 
64
{
 
65
if (m < 1 || m > 3)
 
66
    {
 
67
        QMessageBox::critical((ApplicationWindow *)parent(), tr("SciDAVis") + " - " + tr("Error"),
 
68
        tr("Unknown smooth filter. Valid values are: 1 - Savitky-Golay, 2 - FFT, 3 - Moving Window Average."));
 
69
        d_init_err = true;
 
70
        return;
 
71
    }
 
72
d_method = (SmoothMethod)m;
 
73
}
 
74
 
 
75
void SmoothFilter::calculateOutputData(double *x, double *y)
 
76
{
 
77
    for (int i = 0; i < d_points; i++)
 
78
        {
 
79
           x[i] = d_x[i];
 
80
           y[i] = d_y[i];//filtering frequencies
 
81
        }
 
82
 
 
83
        switch((int)d_method)
 
84
        {
 
85
                case 1:
 
86
            d_explanation = QString::number(d_smooth_points) + " " + tr("points") + " " + tr("Savitzky-Golay smoothing");
 
87
            smoothSavGol(x, y);
 
88
                        break;
 
89
                case 2:
 
90
            d_explanation = QString::number(d_smooth_points) + " " + tr("points") + " " + tr("FFT smoothing");
 
91
                smoothFFT(x, y);
 
92
                        break;
 
93
                case 3:
 
94
            d_explanation = QString::number(d_smooth_points) + " " + tr("points") + " " + tr("average smoothing");
 
95
                smoothAverage(x, y);
 
96
                        break;
 
97
        }
 
98
}
 
99
 
 
100
void SmoothFilter::smoothFFT(double *x, double *y)
 
101
{
 
102
        gsl_fft_real_workspace *work = gsl_fft_real_workspace_alloc(d_n);
 
103
        gsl_fft_real_wavetable *real = gsl_fft_real_wavetable_alloc(d_n);
 
104
        gsl_fft_real_transform (y, 1, d_n, real, work);//FFT forward
 
105
        gsl_fft_real_wavetable_free (real);
 
106
 
 
107
        double df = 1.0/(double)(x[1] - x[0]);
 
108
        double lf = df/(double)d_smooth_points;//frequency cutoff
 
109
        df = 0.5*df/(double)d_n;
 
110
 
 
111
    for (int i = 0; i < d_n; i++)
 
112
        {
 
113
           x[i] = d_x[i];
 
114
           y[i] = i*df > lf ? 0 : y[i];//filtering frequencies
 
115
        }
 
116
 
 
117
        gsl_fft_halfcomplex_wavetable *hc = gsl_fft_halfcomplex_wavetable_alloc (d_n);
 
118
        gsl_fft_halfcomplex_inverse (y, 1, d_n, hc, work);//FFT inverse
 
119
        gsl_fft_halfcomplex_wavetable_free (hc);
 
120
        gsl_fft_real_workspace_free (work);
 
121
}
 
122
 
 
123
void SmoothFilter::smoothAverage(double *, double *y)
 
124
{
 
125
        int p2 = d_smooth_points/2;
 
126
        double m = double(2*p2+1);
 
127
        double aux = 0.0;
 
128
    double *s = new double[d_n];
 
129
 
 
130
        s[0] = y[0];
 
131
        for (int i=1; i<p2; i++)
 
132
        {
 
133
                aux = 0.0;
 
134
                for (int j=-i; j<=i; j++)
 
135
                        aux += y[i+j];
 
136
 
 
137
                s[i] = aux/(double)(2*i+1);
 
138
        }
 
139
        for (int i=p2; i<d_n-p2; i++)
 
140
        {
 
141
                aux = 0.0;
 
142
                for (int j=-p2; j<=p2; j++)
 
143
                        aux += y[i+j];
 
144
 
 
145
                s[i] = aux/m;
 
146
        }
 
147
        for (int i=d_n-p2; i<d_n-1; i++)
 
148
        {
 
149
                aux = 0.0;
 
150
                for (int j=d_n-i-1; j>=i-d_n+1; j--)
 
151
                        aux += y[i+j];
 
152
 
 
153
                s[i] = aux/(double)(2*(d_n-i-1)+1);
 
154
        }
 
155
        s[d_n-1] = y[d_n-1];
 
156
 
 
157
    for (int i = 0; i<d_n; i++)
 
158
        y[i] = s[i];
 
159
 
 
160
    delete[] s;
 
161
}
 
162
 
 
163
void SmoothFilter::smoothSavGol(double *, double *y)
 
164
{
 
165
        double *s = new double[d_n];
 
166
    int nl = d_smooth_points;
 
167
    int nr = d_sav_gol_points;
 
168
        int np = nl+nr+1;
 
169
        double *c = vector(1, np);
 
170
 
 
171
        //seek shift index for given case nl, nr, m (see savgol).
 
172
        int *index = intvector(1, np);
 
173
        index[1]=0;
 
174
        int i, j=3;
 
175
        for (i=2; i<=nl+1; i++)
 
176
        {// index(2)=-1; index(3)=-2; index(4)=-3; index(5)=-4; index(6)=-5
 
177
                index[i]=i-j;
 
178
                j += 2;
 
179
        }
 
180
        j=2;
 
181
        for (i=nl+2; i<=np; i++)
 
182
        {// index(7)= 5; index(8)= 4; index(9)= 3; index(10)=2; index(11)=1
 
183
                index[i]=i-j;
 
184
                j += 2;
 
185
        }
 
186
 
 
187
        //calculate Savitzky-Golay filter coefficients.
 
188
        savgol(c, np, nl, nr, 0, d_polynom_order);
 
189
 
 
190
        for (i=0; i<d_n; i++)
 
191
        {// Apply filter to input data.
 
192
                s[i]=0.0;
 
193
                for (j=1; j<=np; j++)
 
194
                {
 
195
                        int it = i+index[j];
 
196
                        if (it >=0 && it < d_n)//skip left points that do not exist.
 
197
                                s[i] += c[j]*y[i+index[j]];
 
198
                }
 
199
        }
 
200
 
 
201
    for (i = 0; i<d_n; i++)
 
202
        y[i] = s[i];
 
203
 
 
204
        delete[] s;
 
205
        free_vector(c, 1, np);
 
206
        free_intvector(index, 1, np);
 
207
}
 
208
 
 
209
void SmoothFilter::setSmoothPoints(int points, int left_points)
 
210
{
 
211
    if (points < 0 || left_points < 0)
 
212
    {
 
213
        QMessageBox::critical((ApplicationWindow *)parent(), tr("SciDAVis") + " - " + tr("Error"),
 
214
                                tr("The number of points must be positive!"));
 
215
                d_init_err = true;
 
216
                return;
 
217
    }
 
218
    else if (d_polynom_order > points + left_points)
 
219
    {
 
220
        QMessageBox::critical((ApplicationWindow *)parent(), tr("SciDAVis") + " - " + tr("Error"),
 
221
                                tr("The polynomial order must be lower than the number of left points plus the number of right points!"));
 
222
                d_init_err = true;
 
223
                return;
 
224
    }
 
225
 
 
226
    d_smooth_points = points;
 
227
    d_sav_gol_points = left_points;
 
228
}
 
229
 
 
230
void SmoothFilter::setPolynomOrder(int order)
 
231
{
 
232
        if (d_method != SavitzkyGolay)
 
233
    {
 
234
        QMessageBox::critical((ApplicationWindow *)parent(), tr("SciDAVis") + " - " + tr("Error"),
 
235
                                tr("Setting polynomial order is only available for Savitzky-Golay smooth filters! Ignored option!"));
 
236
                return;
 
237
    }
 
238
        
 
239
    if (order > d_smooth_points + d_sav_gol_points)
 
240
    {
 
241
        QMessageBox::critical((ApplicationWindow *)parent(), tr("SciDAVis") + " - " + tr("Error"),
 
242
                                tr("The polynomial order must be lower than the number of left points plus the number of right points!"));
 
243
                d_init_err = true;
 
244
                return;
 
245
    }
 
246
    d_polynom_order = order;
 
247
}