RealFFT::groupDelay now computes the group delay of a signal.

.gitignore

@@ -47,4 +47,5 @@ node_modules
 js/libgtkIOStreamNode.js
 webcomponents/gtkiostream/libgtkIOStream.js
 webcomponents/gtkiostream/libgtkIOStreamNode.js
-webcomponents/gtkiostream/test
+webcomponents/gtkiostream/test
+RealFFTExampleGD

include/fft/RealFFT.H

@@ -57,6 +57,13 @@ public:
 
     /// Inverse transform the data (out to in)
     void invTransform();
+
+    /** Find the group delay of a real fft object
+    Uses the Smith algorithm for computing the group delay.
+    \param rfd The DFT data to find the group delay of
+    \returns The group delay in the form of a RealFFTData object where the RealFFTData::in variable is the group delay
+    */
+    RealFFTData groupDelay(RealFFTData &rfd);
 };
 /** \example RealFFTExample.C
  * This is an example of how to use the class.

src/RealFFT.C

@@ -70,6 +70,17 @@ void RealFFT::invTransform() {
         fftw_execute(invPlan);
 }
 
+RealFFTData RealFFT::groupDelay(RealFFTData &rfd){
+  RealFFTData gd(rfd.getSize()); // correctly size the group delay object
+  for (int i=0; i<rfd.getSize(); i++)
+    gd.in[i]=rfd.in[i]*i;
+  switchData(gd); // derive the DFT coefficients
+  fwdTransform();
+  for (int i=0; i<rfd.getSize(); i++)
+    gd.in[i]=(gd.getComplexCoeff(i)/rfd.getComplexCoeff(i)).real();
+  return gd;
+}
+
 #include "gtkiostream_config.h"
 #ifdef HAVE_EMSCRIPTEN
 #include <emscripten/bind.h>

test/RealFFTExampleGD.C

@@ -50,6 +50,8 @@ int main (void){
     // forward transform :
     rfft.fwdTransform();
 
+    RealFFTData gd=rfft.groupDelay(fftData);
+
     // cout<<"y=["<<endl;
     // fftData.dumpOut();
     // cout<<"];"<<endl;
@@ -60,9 +62,19 @@ int main (void){
     }
     cout<<"];"<<endl;
 
+    cout<<"gd=[";
+    for (unsigned int i=0; i<N; i++){
+      cout<<gd.in[i]<<' ';
+    }
+    cout<<"];"<<endl;
+
     cout<<"XHAT=fft(x);"<<endl;
     cout<<"[X; XHAT]"<<endl;
-    cout<<"sum(X-XHAT)/length(X)"<<endl;
+    cout<<"complexError=sum(X-XHAT)/length(X)"<<endl;
+
+    cout<<"gdHat=grpdelay(x,1,"<<N<<",'whole')';"<<endl;
+    cout<<"[gd; gdHat]"<<endl;
+    cout<<"gdError=sum(gd-gdHat)/length(X)"<<endl;
 
     // cout<<"\""<<endl;
   return 0;