1 /* -------------------------------------------------------------------------- *
2 * OpenSim: Signal.cpp *
3 * -------------------------------------------------------------------------- *
4 * The OpenSim API is a toolkit for musculoskeletal modeling and simulation. *
5 * See http://opensim.stanford.edu and the NOTICE file for more information. *
6 * OpenSim is developed at Stanford University and supported by the US *
7 * National Institutes of Health (U54 GM072970, R24 HD065690) and by DARPA *
8 * through the Warrior Web program. *
9 * *
10 * Copyright (c) 2005-2017 Stanford University and the Authors *
11 * Author(s): Frank C. Anderson *
12 * *
13 * Licensed under the Apache License, Version 2.0 (the "License"); you may *
14 * not use this file except in compliance with the License. You may obtain a *
15 * copy of the License at http://www.apache.org/licenses/LICENSE-2.0. *
16 * *
17 * Unless required by applicable law or agreed to in writing, software *
18 * distributed under the License is distributed on an "AS IS" BASIS, *
19 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. *
20 * See the License for the specific language governing permissions and *
21 * limitations under the License. *
22 * -------------------------------------------------------------------------- */
23
24 /* Note: This code was originally developed by Realistic Dynamics Inc.
25 * Author: Frank C. Anderson
26 */
27
28
29 // INCLUDES
30 #include <math.h>
31 #include "Signal.h"
32 #include "Array.h"
33 #include "SimTKcommon/Constants.h"
34 #include "SimTKcommon/Orientation.h"
35 #include "SimTKcommon/Scalar.h"
36 #include "SimTKcommon/SmallMatrix.h"
37 #include "simmath/internal/Spline.h"
38 #include "simmath/internal/SplineFitter.h"
39
40 using namespace OpenSim;
41 using namespace std;
42
43 //=============================================================================
44 // FILTERS
45 //=============================================================================
46 //-----------------------------------------------------------------------------
47 // GCVSPL
48 //-----------------------------------------------------------------------------
49 //_____________________________________________________________________________
50 /**
51 * Generalized, cross-validatory spline smoothing (Adv. Engng. Softw. 8:104-113)
52 *
53 * @param degree Degree of the spline.
54 * @param T Sample interval in seconds.
55 * @param fc Cutoff frequency in Hz.
56 * @param N Number of data points in the signal.
57 * @param times The times for a specified signal samples.
58 * @param sig The sampled signal.
59 * @param sigf The filtered signal.
60 *
61 * @return 0 on success, and -1 on failure.
62 */
63 int Signal::
SmoothSpline(int degree,double T,double fc,int N,double * times,double * sig,double * sigf)64 SmoothSpline(int degree,double T,double fc,int N,double *times,double *sig,double *sigf)
65 {
66 /* Background for choice of smoothing parameter (by Ton van den Bogert):
67
68 In the GCVSPL software, VAL is the "smoothing parameter" of the spline,
69 which occurs as a weighting factor in a cost function which is
70 minimized. This produces a compromise between how well the spline
71 fits the measurements and its smoothness (for a more detailed
72 explanation, see Ton van den Bogert's notes in
73 http://isbweb.org/software/sigproc/bogert/filter.pdf).
74
75 Woltring describes the relationship between VAL and frequency domain
76 filter characteristics in his release notes:
77 http://isbweb.org/software/sigproc/gcvspl/gcvspl.memo.
78 This states that a smoothing spline is equivalent to a
79
80 "double, phase-symmetric Butterworth filter, with transfer function
81 H(w) = [1 + (w/wo)^2M]^-1, where w is the frequency,
82 wo = (p*T)^(-0.5/M) the filter's cut-off frequency, p the smoothing
83 parameter, T the sampling interval, and 2M the order of the spline.
84 If T is expressed in seconds, the frequencies are expressed in
85 radians/second."
86
87 The VAL calculations can be derived from this. The equation for the
88 transfer function has the property H(w0)= 0.5. This is because it is
89 a double Butterworth filter (applied twice). Cut-off frequency is
90 usually defined as the frequency at which H=1/sqrt(2), and this is why
91 Tony Reina and Ton van den Bogert have an extra factor in the equation.
92
93 You can easily verify the correctness of your spline smoothing by
94 processing a sine wave signal at the cut-off frequency. If it comes
95 out as a sine wave with its amplitude reduced by a factor 1.41 (=sqrt(2)),
96 you have used the correct VAL.
97
98 */
99
100 SimTK::Vector x(N);
101 SimTK::Vector_<SimTK::Vec<1> > y(N);
102 for (int i = 0; i < N; ++i)
103 x[i] = times[i];
104 for (int i = 0; i < N; ++i)
105 y[i] = SimTK::Vec<1>(sig[i]);
106
107 int M = (degree+1)/2; // Half-order
108 SimTK::Real p; // Smoothing parameter
109 p = (1.0/T) /
110 (pow(
111 (2.0*SimTK_PI*fc) /
112 (pow(
113 (sqrt(2.0)-1),
114 (0.5/M)
115 )),
116 (2*M)));
117
118 // SMOOTH SPLINE FIT TWICE TO SIMULATE DOUBLE, PHASE-SYMMETRIC BUTTERWORTH FILTER
119 SimTK::Spline_<SimTK::Vec<1> > smoothSpline1;
120 smoothSpline1 = SimTK::SplineFitter<SimTK::Vec<1> >::fitForSmoothingParameter(degree,x,y,p).getSpline();
121 SimTK::Spline_<SimTK::Vec<1> > smoothSpline2;
122 smoothSpline2 = SimTK::SplineFitter<SimTK::Vec<1> >::fitForSmoothingParameter(degree,x,smoothSpline1.getControlPointValues(),p).getSpline();
123
124 for (int i = 0; i < N; ++i)
125 sigf[i] = smoothSpline2.getControlPointValues()[i][0];
126
127 // CHECK THAT P IS NOT BEYOND BOUND
128 SimTK::Real pActual;
129 pActual = SimTK::SplineFitter<SimTK::Vec<1> >::fitForSmoothingParameter(degree,x,smoothSpline2.getControlPointValues(),p).getSmoothingParameter();
130 //cout << "Requested smoothing parameter = " << p << endl;
131 //cout << "Actual smoothing parameter = " << pActual << endl;
132 if(p!=pActual) {
133 printf("Signal.SmoothSpline: ERROR- The cutoff frequency (%lf)",fc);
134 printf(" produced a smoothing parameter (%le) beyond its bound (%le).\n",p,pActual);
135 return(-1);
136 }
137
138 return(0);
139 }
140
141 //-----------------------------------------------------------------------------
142 // IIR
143 //-----------------------------------------------------------------------------
144 //_____________________________________________________________________________
145 /**
146 * 3rd ORDER LOWPASS IIR BUTTERWORTH DIGITAL FILTER
147 *
148 * It is assumed that enough memory is allocated at sigf.
149 * Note also that the first and last three data points are not filtered.
150 *
151 * @param T Sample interval in seconds.
152 * @param fc Cutoff frequency in Hz.
153 * @param N Number of data points in the signal.
154 * @param sig The sampled signal.
155 * @param sigf The filtered signal.
156 *
157 * @return 0 on success, and -1 on failure.
158 */
159 int Signal::
LowpassIIR(double T,double fc,int N,double * sig,double * sigf)160 LowpassIIR(double T,double fc,int N,double *sig,double *sigf)
161 {
162 int i,j;
163 double fs/*,ws*/,wc,wa,wa2,wa3;
164 double a[4],b[4],denom;
165 double *sigr;
166
167 // ERROR CHECK
168 if(T==0) return(-1);
169 if(N==0) return(-1);
170 if(sig==NULL) return(-1);
171 if(sigf==NULL) return(-1);
172
173 // CHECK THAT THE CUTOFF FREQUENCY IS LESS THAN HALF THE SAMPLE FREQUENCY
174 fs = 1 / T;
175 if (fc >= 0.5 * fs) {
176 printf("\nCutoff frequency should be less than half sample frequency.");
177 printf("\nchanging the cutoff frequency to 0.49*(Sample Frequency)...");
178 fc = 0.49 * fs;
179 printf("\ncutoff = %lf\n\n",fc);
180 }
181
182 // INITIALIZE SOME VARIABLES
183 //ws = 2*SimTK_PI*fs;
184 wc = 2*SimTK_PI*fc;
185
186 // CALCULATE THE FREQUENCY WARPING
187 wa = tan(wc*T/2.0);
188 wa2 = wa*wa;
189 wa3 = wa*wa*wa;
190
191 // GET COEFFICIENTS FOR THE FILTER
192 denom = (wa+1) * (wa*wa + wa + 1.0);
193 a[0] = wa3 / denom;
194 a[1] = 3*wa3 / denom;
195 a[2] = 3*wa3 / denom;
196 a[3] = wa3 / denom;
197 b[0] = 1;
198 b[1] = (3*wa3 + 2*wa2 - 2*wa - 3) / denom;
199 b[2] = (3*wa3 - 2*wa2 - 2*wa + 3) / denom;
200 b[3] = (wa - 1) * (wa2 - wa + 1) / denom;
201
202 // ALLOCATE MEMORY FOR sigr[]
203 sigr = new double[N];
204 if(sigr==NULL) {
205 printf("\nSignal.lowpassIIR: ERROR- Not enough memory.\n");
206 return(-1);
207 }
208
209 // FILTER THE DATA
210 // FILL THE 1ST THREE TERMS OF sigf
211 for (i=0;i<=3;i++) sigf[i] = sig[i];
212
213 // IMPLEMENT THE FORMULA
214 for (i=3;i<N;i++) {
215 sigf[i] = a[0]*sig[i] + a[1]*sig[i-1] + a[2]*sig[i-2] + a[3]*sig[i-3]
216 - b[1]*sigf[i-1] - b[2]*sigf[i-2] - b[3]*sigf[i-3];
217 }
218
219 // REVERSE THE FILTERED ARRAY
220 for (i=0,j=N-1;i<N;i++,j--) sigr[i] = sigf[j];
221
222 // FILL THE 1ST THREE TERMS OF sigf
223 for (i=0;i<=3;i++) sigf[i] = sigr[i];
224
225 // IMPLEMENT THE FORMULA AGAIN
226 for (i=3;i<N;i++) {
227 sigf[i] = a[0]*sigr[i] + a[1]*sigr[i-1] + a[2]*sigr[i-2] + a[3]*sigr[i-3]
228 - b[1]*sigf[i-1] - b[2]*sigf[i-2] - b[3]*sigf[i-3];
229 }
230
231 // REVERSE THE FILTERED ARRAY AGAIN
232 for (i=0,j=N-1;i<N;i++,j--) sigr[i] = sigf[j];
233
234 // ASSIGN sigf TO sigr
235 for (i=0;i<N;i++) sigf[i] = sigr[i];
236
237 // CLEANUP
238 if(sigr!=NULL) delete[] sigr;
239
240 return(0);
241 }
242
243 //-----------------------------------------------------------------------------
244 // FIR
245 //-----------------------------------------------------------------------------
246 //_____________________________________________________________________________
247 /**
248 * LOWPASS FIR NON-RECURSIVE DIGITAL FILTER
249 *
250 * It is permissible for sig and sigf to be the same array or overlap.
251 *
252 * PARAMETERS
253 * @param M Order of filter (should be 30 or greater).
254 * @param T Sample interval in seconds.
255 * @param f Cutoff frequency in Hz.
256 * @param N Number of data points in the signal.
257 * @param sig The sampled signal.
258 * @param sigf The filtered signal.
259 *
260 * @return 0 on success, and -1 on failure.
261 */
262 int Signal::
LowpassFIR(int M,double T,double f,int N,double * sig,double * sigf)263 LowpassFIR(int M,double T,double f,int N,double *sig,double *sigf)
264 {
265 int n,k;
266 double w,x;
267
268 // CHECK THAT M IS NOT TOO LARGE RELATIVE TO N
269 if((M+M)>N) {
270 printf("rdSingal.lowpassFIR: ERROR- The number of data points (%d)",N);
271 printf(" should be at least twice the order of the filter (%d).\n",M);
272 return(-1);
273 }
274
275 // PAD THE SIGNAL SO FILTERING CAN BEGIN AT THE FIRST DATA POINT
276 double *s = Pad(M,N,sig);
277 if(s==NULL) return(-1);
278
279 // CALCULATE THE ANGULAR CUTOFF FREQUENCY
280 w = 2.0*SimTK_PI*f;
281
282 // FILTER THE DATA
283 double sum_coef,coef;
284 for(n=0;n<N;n++) {
285 sum_coef = 0.0;
286 sigf[n] = 0.0;
287 for(k=-M;k<=M;k++) {
288 x = (double)k*w*T; // k*T = time (seconds) and w scales sinc input argument using filter cutoff
289 coef = (sinc(x)*T*w/SimTK_PI)*hamming(k,M); // scale lowpass sinc amplitude by 2*f*T = T*w/pi
290 sigf[n] = sigf[n] + coef*s[M+n-k];
291 sum_coef = sum_coef + coef;
292 }
293 sigf[n] = sigf[n] / sum_coef; // normalize for unity gain at DC
294 }
295
296 // Filter check derived from http://www.dspguide.com/CH16.PDF
297 //double sum_coef,coef;
298 //for(n=0;n<N;n++) {
299 // sum_coef = 0.0;
300 // sigf[n] = 0.0;
301 // for(k=-M;k<=M;k++) {
302 // if(k==0) {
303 // coef = 2.0*SimTK_PI*f;
304 // } else {
305 // coef = sin((double)k*2.0*SimTK_PI*f) / (double)k;
306 // }
307 // coef = coef*hamming(k,M);
308 // sigf[n] = sigf[n] + coef*s[M+n-k];
309 // sum_coef = sum_coef + coef;
310 // }
311 // sigf[n] = sigf[n] / sum_coef;
312 //}
313
314 // CLEANUP
315 delete[] s;
316
317 return(0);
318 }
319
320
321
322 //_____________________________________________________________________________
323 /**
324 * BANDPASS FIR NON-RECURSIVE DIGITAL FILTER
325 *
326 * Note that sig and sigf must point to distinct locations in memory which do
327 * not overlap.
328 *
329 * PARAMETERS
330 * @param M Order of filter (should be 30 or greater).
331 * @param T Sample interval in seconds.
332 * @paramf1 lowend cutoff frequency in Hz
333 * @paramf2 highend cutoff frequency in Hz
334 * @param N Number of data points in the signal.
335 * @param sig The sampled signal.
336 * @param sigf The filtered signal.
337 *
338 * @return 0 on success, and -1 on failure.
339 */
340 int Signal::
BandpassFIR(int M,double T,double f1,double f2,int N,double * sig,double * sigf)341 BandpassFIR(int M,double T,double f1,double f2,int N,double *sig,
342 double *sigf)
343 {
344 size_t size;
345 int i,j;
346 int n,k;
347 double w1,w2,x1,x2;
348 double *s;
349
350
351 // CHECK THAT M IS NOT TOO LARGE RELATIVE TO N
352 if((M+M)>N) {
353 printf("\n\nThe number of data points (%d) should be at least twice\n",N);
354 printf("the order of the filter (%d).\n\n",M);
355 return(-1);
356 }
357
358 // ALLOCATE MEMORY FOR s
359 size = N + M + M;
360 s = (double *) calloc(size,sizeof(double));
361 if (s == NULL) {
362 printf("\n\nlowpass() -> Not enough memory to process your sampled data.");
363 return(-1);
364 }
365
366 // CALCULATE THE ANGULAR CUTOFF FREQUENCY
367 w1 = 2*SimTK_PI*f1;
368 w2 = 2*SimTK_PI*f2;
369
370 // PAD THE SIGNAL SO FILTERING CAN BEGIN AT THE FIRST DATA POINT
371 for (i=0,j=M;i<M;i++,j--) s[i] = sig[j];
372 for (i=M,j=0;i<M+N;i++,j++) s[i] = sig[j];
373 for (i=M+N,j=N-2;i<M+M+N;i++,j--) s[i] = sig[j];
374
375
376 // FILTER THE DATA
377 double sum_coef,coef;
378 for (n=0;n<N;n++) {
379 sum_coef = 0.0;
380 sigf[n] = 0.0;
381 for (k=-M;k<=M;k++) {
382 x1 = (double)k*w1*T; // k*T = time (seconds) and w scales sinc input argument using filter cutoff
383 x2 = (double)k*w2*T; // k*T = time (seconds) and w scales sinc input argument using filter cutoff
384 coef = (sinc(x2)*T*w2/SimTK_PI - sinc(x1)*T*w1/SimTK_PI)*hamming(k,M); // scale lowpass sinc amplitude by 2*f*T = T*w/pi
385 sigf[n] = sigf[n] + coef*s[M+n-k];
386 sum_coef = sum_coef + coef;
387 }
388 sigf[n] = sigf[n] / sum_coef; // normalize for unity gain at DC
389 }
390
391 // CLEANUP
392 if(s!=NULL) delete s;
393
394 return(0);
395 }
396
397 //_____________________________________________________________________________
398 /**
399 * Pad a signal with a specified number of data points.
400 *
401 * The signal is prepended and appended with a reflected and negated
402 * portion of the signal of the appropriate size so as to preserve the value
403 * and slope of the signal.
404 *
405 * PARAMETERS
406 * @param aPad Size of the pad-- number of points to prepend and append.
407 * @param aN Number of data points in the signal.
408 * @param aSignal Signal to be padded.
409 * @return Padded signal. The size is aN + 2*aPad. NULL is returned
410 * on an error. The caller is responsible for deleting the returned
411 * array.
412 */
413 double* Signal::
Pad(int aPad,int aN,const double aSignal[])414 Pad(int aPad,int aN,const double aSignal[])
415 {
416 if(aPad<=0) return(NULL);
417
418 // COMPUTE FINAL SIZE
419 int size = aN + 2*aPad;
420 if(aPad>aN) {
421 cout<<"\nSignal.Pad(double[]): ERROR- requested pad size ("<<aPad<<") is greater than the number of points ("<<aN<<").\n";
422 return(NULL);
423 }
424
425 // ALLOCATE
426 double *s = new double[size];
427 if (s == NULL) {
428 printf("\n\nSignal.Pad: Failed to allocate memory.\n");
429 return(NULL);
430 }
431
432 // PREPEND
433 int i,j;
434 for(i=0,j=aPad;i<aPad;i++,j--) s[i] = aSignal[j];
435 for(i=0;i<aPad;i++) s[i] = 2.0*aSignal[0] - s[i];
436
437 // SIGNAL
438 for(i=aPad,j=0;i<aPad+aN;i++,j++) s[i] = aSignal[j];
439
440 // APPEND
441 for(i=aPad+aN,j=aN-2;i<aPad+aPad+aN;i++,j--) s[i] = aSignal[j];
442 for(i=aPad+aN;i<aPad+aPad+aN;i++) s[i] = 2.0*aSignal[aN-1] - s[i];
443
444 return(s);
445 }
446 //_____________________________________________________________________________
447 /**
448 * Pad a signal with a specified number of data points.
449 *
450 * The signal is prepended and appended with a reflected and negated
451 * portion of the signal of the appropriate size so as to preserve the value
452 * and slope of the signal.
453 *
454 * PARAMETERS
455 * @param aPad Size of the pad-- number of points to prepend and append.
456 * @param rSignal Signal to be padded.
457 */
458 void Signal::
Pad(int aPad,Array<double> & rSignal)459 Pad(int aPad,Array<double> &rSignal)
460 {
461 if(aPad<=0) return;
462
463 // COMPUTE NEW SIZE
464 int size = rSignal.getSize();
465 int newSize = size + 2*aPad;
466
467 // HANDLE PAD GREATER THAN SIZE
468 if(aPad>=size) {
469 int pad = size - 1;
470 while(size<newSize) {
471 Pad(pad,rSignal);
472 size = rSignal.getSize();
473 pad = (newSize - size) / 2;
474 if(pad>=size) pad = size - 1;
475 }
476 return;
477 }
478
479 // ALLOCATE
480 Array<double> s(0.0,newSize);
481
482 // PREPEND
483 int i,j;
484 for(i=0,j=aPad;i<aPad;i++,j--) s[i] = rSignal[j];
485 for(i=0;i<aPad;i++) s[i] = 2.0*rSignal[0] - s[i];
486
487 // SIGNAL
488 for(i=aPad,j=0;i<aPad+size;i++,j++) s[i] = rSignal[j];
489
490 // APPEND
491 for(i=aPad+size,j=size-2;i<aPad+aPad+size;i++,j--) s[i] = rSignal[j];
492 for(i=aPad+size;i<aPad+aPad+size;i++) s[i] = 2.0*rSignal[size-1] - s[i];
493
494 // ALTER SIGNAL
495 rSignal = s;
496 }
497
498
499 //-----------------------------------------------------------------------------
500 // POINT REDUCTION
501 //-----------------------------------------------------------------------------
502 //_____________________________________________________________________________
503 /**
504 * Remove points in a signal based on the angle between adjacent segments in
505 * the signal. The end points of the signal are always retained.
506 *
507 * @param aAngle If the angle between two adjacent segments is less than
508 * aAngle the point in common between the two segments is removed. This
509 * is evaluate for each point in the signal.
510 * @param rTime Array of time values. This array is altered.
511 * @param rSignal Array of signal values. This array is altered.
512 * @return Number of points removed.
513 */
514 int Signal::
ReduceNumberOfPoints(double aDistance,Array<double> & rTime,Array<double> & rSignal)515 ReduceNumberOfPoints(double aDistance,
516 Array<double> &rTime,Array<double> &rSignal)
517 {
518 // CHECK SIZES
519 int size = rTime.getSize();
520 int sizeSignal = rSignal.getSize();
521 if(size!=sizeSignal) {
522 cout<<"\n\nSignal.ReduceNumberOfPoints:: quitting. The time and signal ";
523 cout<<"arrays have different numbers of points.\n\n";
524 return(0);
525 }
526
527 // CHECK ANGLE
528 if(aDistance<SimTK::Zero) aDistance = SimTK::Zero;
529
530 // APPEND FIRST POINT
531 Array<double> t(0.0,0,size),s(0.0,0,size);
532 t.append(rSignal[0]);
533 s.append(rSignal[0]);
534 int iLast=0;
535
536 // REMOVE POINTS
537 int i,imid;
538 SimTK::Vec3 p1(0.0,0.0,0.0);
539 SimTK::Vec3 p2(0.0,0.0,0.0);
540 SimTK::Vec3 p3(0.0,0.0,0.0);
541 SimTK::Vec3 v1(0.0,0.0,0.0);
542 SimTK::Vec3 v2(0.0,0.0,0.0);
543
544 double tmid,mv1,mv2,cos,dsq;
545 for(i=1;i<(size-1);i++) {
546
547 // FIRST POINT
548 // LAST POINT IN SIMPLIFIED ARRAY
549 p1[0] = t.getLast();
550 p1[1] = s.getLast();
551
552 // THIRD POINT
553 p3[0] = rTime[i+1];
554 p3[1] = rSignal[i+1];
555
556 // SECOND POINT
557 // MIDPOINT
558 tmid = 0.5*(p3[0]-p1[0]) + p1[0];
559 imid = rTime.searchBinary(tmid,false,iLast,i+1);
560 if(imid<=iLast) imid++;
561 p2[0] = rTime[imid];
562 p2[1] = rSignal[imid];
563
564 v1 = p2 - p1; //Mtx::Subtract(1,3,p2,p1,v1);
565 v2 = p3 - p1; //Mtx::Subtract(1,3,p3,p1,v2);
566
567 mv1 = v1.norm(); //Mtx::Magnitude(3,v1);
568 mv2 = v2.norm(); //Mtx::Magnitude(3,v2);
569 cos = (~v1*v2)/(mv1*mv2); //Mtx::DotProduct(3,v1,v2) / (mv1*mv2);
570
571 dsq = mv1 * mv1 * (1.0 - cos*cos);
572
573 if(dsq < (aDistance*aDistance)) continue;
574
575 iLast = i;
576 t.append(rTime[i]);
577 s.append(rSignal[i]);
578 }
579
580 // APPEND ENDPOINT
581 t.append(rTime.getLast());
582 s.append(rSignal.getLast());
583
584 // NUMBER REMOVED
585 int numberRemoved = rTime.getSize() - t.getSize();
586
587 // ALTER ARRAYS
588 rTime = t;
589 rSignal = s;
590
591 return(numberRemoved);
592 }
593
594
595
596 //=============================================================================
597 // CORE MATH
598 //=============================================================================
599 //_____________________________________________________________________________
600 /**
601 * SINC(X) = SIN(X) / X
602 *
603 * @return sin(x)/x.
604 */
605 double Signal::
sinc(double x)606 sinc(double x)
607 {
608 if ((x<1.0e-8) && (x>-1.0e-8)) return(1.0);
609 return(sin(x)/x);
610 }
611 //_____________________________________________________________________________
612 /**
613 * Hamming Window- dampens Gibbs phenomenon in filtering functions.
614 */
615 double Signal::
hamming(int k,int M)616 hamming(int k,int M)
617 {
618 double x = (double)k * SimTK_PI / (double)M;
619 double d = 0.54 + 0.46*cos(x);
620 return(d);
621 }
622
623
624