Recreating Convolution Animation on Wikipedia

Hey All,

I’m trying to recreate the convolution animation for the two square waves on Wikipedia found here:
en.wikipedia.org/wiki/Convolution

I’ve pretty well done everything but I can’t get the convolution function (fg)(t) to be drawn progressively as the waves pass each other. Instead ROOT draws (fg)(t) over the entire range of the plot resulting as a constant line at the value of the integral of the product of the two functions.

I think this is due to my in ability to find a way to make the function convolvedFunc() (see code below) depend on x. I could store the values of the integral in a TGraph and then draw it afterwards, but I want to be able to draw the resulting convoluted function in real time and use the resulting TF1 to do fits.

Thanks for your thoughts!
-Chris

//This can be run in ROOT by doing .x convolution.C

#include <iostream>
using namespace std;

//Instantiate the functions to be used
//Their properties will be set later
TF1 *func1;  //The stationary Function
TF1 *func2;  //The Function to be moved 
TF1 *func3;  //The product of the two functions
TF1 *func4;  //The Resulting function

//Square Wave
Double_t squareWave(Double_t *x, Double_t *par){
  Double_t xx = x[0];
  if (xx < par[0]-par[1]/2.)
    return 0;
  else if (xx >= par[0]-par[1]/2. && xx <= par[0]+par[1]/2.)
    return par[2];
  else if (xx > par[0]+par[1]/2.)
    return 0;

} //End Square Wave

//Does the multiplication of the two functions
Double_t multipliedFuncs(Double_t *x, Double_t *par){
  
  Double_t xx = x[0];
  Double_t f_x = func1->Eval(xx)*func2->Eval(xx);
  return f_x;

}//End multipliedFuncs

//The resulting function of the convolution
Double_t convolvedFunc(Double_t *x, Double_t *par){
  return func3->Integral(-2,2);
}

void convolution(){

  //Create the First Function
  func1 = new TF1("SquareWave",squareWave,-2,2,3);
  func1->SetNpx(10000);
  
  //Set the Parameters
  func1->SetParameter(0,0);  //Center of the Square
  func1->SetParameter(1,1);  //Width of the Square
  func1->SetParameter(2,1);  //Height of the Square

  //Create the Second Function
  func2 = new TF1("SquareWave2",squareWave,-2,2,3);
  func2->SetNpx(10000);

  //Set the Parameters
  func2->SetParameter(0,-2);  //Starting Center of the Square
  func2->SetParameter(1,1);   //Width of the Square
  func2->SetParameter(2,1);   //Height of the Square

  //Amount by which to move the second function every step
  Double_t moveOver = 0 .05;

  //Create the Multiplied Function
  func3 = new TF1("multipliedFuncs",multipliedFuncs,-2,2,0);
  func3->SetNpx(10000);
  func3->SetLineColor(kYellow);
  func3->SetLineWidth(0);
  func3->SetFillColor(kYellow);
  func3->SetFillStyle(1001);

  //Create Convolution Function
  func4 = new TF1("convolvedFunc",convolvedFunc,-2,2,0);
  func4->SetLineColor(kBlack);

  //Loop to move the second function toward the first
  for (Int_t i=0; i<75; i++){

    //Set the Mean of the Guassian - This moves the function
    func2->SetParameter(0,func2->GetParameter(0)+moveOver);

    //Draw the Functions
    func1->Draw();
    func2->Draw("SAME");
    func3->Draw("SAME");
    func4->Draw("SAME");
    gPad->Update();

  }
}
//This can be run in ROOT by doing .x convolution.C

//Instantiate the functions and graphs to be used
//Their properties will be set later
TF1 *func1;  //The stationary Function
TF1 *func2;  //The Function to be moved
TF1 *func3;  //The product of the two functions
TGraph *graph4;  //The Resulting Convolution Graph

//Distorted Gauss Wave
Double_t gaussWave(Double_t *x, Double_t *par) {
  if (par[1] > 0) {
    if (x[0] < par[0])
      return par[2] * TMath::Gaus(x[0], par[0], (2.0 * par[1]), kFALSE);
    else
      return par[2] * TMath::Gaus(x[0], par[0], (0.5 * par[1]), kFALSE);
  } else {
    if (x[0] < par[0])
      return par[2] * TMath::Gaus(x[0], par[0], (-0.5 * par[1]), kFALSE);
    else
      return par[2] * TMath::Gaus(x[0], par[0], (-2.0 * par[1]), kFALSE);
  }
} //End Distorted Gauss Wave

//Does the multiplication of the two functions
Double_t multipliedFuncs(Double_t *x, Double_t *par){
  Double_t xx = x[0];
  Double_t f_x = func1->Eval(xx)*func2->Eval(xx);
  return f_x;
}//End multipliedFuncs

void convolution(){
  //Create the First Function
  func1 = new TF1("GaussWave",gaussWave,-2,2,3);
  func1->SetNpx(1000);
  
  //Set the Parameters
  func1->SetParameter(0, 0.5); //Center of the Gauss
  func1->SetParameter(1, 0.3); //Sigma of the Gauss
  func1->SetParameter(2, 1);   //Height of the Gauss
  
  //Create the Second Function
  func2 = new TF1("GaussWave2",gaussWave,-2,2,3);
  func2->SetNpx(1000);
  
  //Set the Parameters
  func2->SetParameter(0, -2);   //Starting Center of the Gauss
  func2->SetParameter(1, -0.1); //Sigma of the Gauss
  func2->SetParameter(2, 1);    //Height of the Gauss
  
  //Amount by which to move the second function every step
  Double_t moveOver = 0.05;
  
  //Create the Multiplied Function
  func3 = new TF1("multipliedFuncs",multipliedFuncs,-2,2,0);
  func3->SetNpx(1000);
  func3->SetLineColor(kYellow);
  func3->SetLineWidth(0);
  func3->SetFillColor(kYellow);
  func3->SetFillStyle(1001);
  
  //Create Convolution Graph
  graph4 = new TGraph();
  graph4->SetLineColor(kRed);
  graph4->SetLineWidth(3);
  graph4->SetPoint(0, -2, func3->Integral(-2, 2));
  
  //Draw the Functions and Graphs
  func1->Draw("L");
  func2->Draw("L SAME");
  func3->Draw("L SAME");
  graph4->Draw("L SAME");
  gPad->Modified(); gPad->Update();
  
  //Loop to move the second function toward the first
  for (Int_t i = 0; i < 80; i++) {
    //Set the Mean of the Guassian - This moves the function
    Double_t x = func2->GetParameter(0) + moveOver;
    func2->SetParameter(0, x);
    graph4->SetPoint((graph4->GetN()), x, func3->Integral(-2, 2));
    
    //Update the picture
    gPad->Modified(); gPad->Update();
  }
}

P.S. A small warning … the TF1::Integral routine “hates” rapidly changing functions (e.g. “square pulses”).

Brilliant! Thanks!

I wasn’t so interested in the square wave…it was just a simple temp function for use until I could everything working.

Thanks Again!
-Chris

It turns out that there was a small problem with the original code (and as a result Wile’s). The yellow fill color is supposed to represent the overlap region of the two functions being convoluted. It seems to work in the code that Wile provided, but when I changed the functions to the ones I needed the yellow area was incorrect.

The code posted below corrects the issue with the overlap region and represents the final product. My original goal was to show the animation of the convolution of a Breit-Wigner and Gaussian function. For those interested and for reference here it is.

Cheers,
Chris


#include <iostream>
#include <TF1.h>
#include <TCanvas.h>
#include <TAttLine.h>
#include <TLegend.h>
#include <TLegendEntry.h>
#include <TLatex.h>
#include <TGraph.h>


//Instantiate the functions and graphs to be used
//Their properties will be set later
TF1 *func1;       //The function to be translated
TF1 *func2;       //The Function to be reflected
TF1 *func3;       //The product of func1 and func2
TGraph *graph4;   //The Resulting Convolution Graph

//Set Precision of the Integrals (number of steps)
Int_t precision = 100;

//Create an animated gif?
Bool_t createGIF = true;



//Does the multiplication of the two functions
Double_t multipliedFuncs(Double_t *x, Double_t *par){
  Double_t xx = x[0];
  Double_t f_x = func1->Eval(xx)*func2->Eval(xx);
  return f_x;
}//End multipliedFuncs




//Find the larger value between of the functions
//This is used to shade the overlap region
Double_t smallerValue(Double_t x){
  
  //Evaluate the Functions
  Double_t val1 = func1->Eval(x);
  Double_t val2 = func2->Eval(x);
  
  if (val1 < val2)
    return val1;
  else
    return val2;

}//End largerValue




//Makes Graph 5 for each x position
//Graph 5 is responsible for shading the overlap region
TGraph *graph5 = new TGraph();

void makeGraph5(){

  Double_t tempX =-10.0;
  Int_t point = 0;

  //Loop Over all values of x and find the overlap region to shade
  do {

    graph5->SetPoint(point, tempX, smallerValue(tempX));
    
    tempX += 0.05;
    point++;

  } while(tempX < 10.0);

}//End makeGraph5





//Runs the Convolution Animation
void convolutionAnimation(){

  //Breit-Wigner Properties
  Double_t bwMean = 0.0;
  Double_t bwGamma = 2.5;
  
  //Gaussian Properties
  Double_t gWidth = 2.0;

  //Create the Gaussian (func1)
  func1 = new TF1("Gaussian","TMath::Gaus(x,[0],[1],1)",-10,10);
  func1->SetParameters(-10,gWidth);
  func1->SetNpx(precision);
  func1->SetLineColor(kBlack);
  func1->SetLineWidth(3);

  //Define the Breit-Wigner (func2)
  func2 = new TF1("BreitWigner","TMath::BreitWigner(x,[0],[1])",-10,10);
  func2->SetParameters(bwMean,bwGamma);
  func2->SetNpx(precision);
  func2->SetLineColor(kBlue);
  func2->SetLineWidth(3);

  //Define the Product of the two functions
  func3 = new TF1("MultipliedFuncs",multipliedFuncs,-10,10,0);
  func3->SetNpx(precision);

  //Create Convolution Graph
  graph4 = new TGraph();
  graph4->SetName("Convolution");
  graph4->SetLineColor(kRed);
  graph4->SetLineWidth(3);
  graph4->SetPoint(0, -10, func3->Integral(-10, 10));

  //Set Attibutes for Graph5
  graph5->SetFillColor(kYellow);
  graph5->SetFillStyle(1001);

  //Set the Title to be Drawn on the Canvas
  func2->SetTitle("Convolution of a Gaussian and Breit-Wigner");
  
  //Create a Legend
  TLegend *legend = new TLegend(0.55,0.72,0.88,0.88,NULL,"brNDC");
  TLegendEntry *entry = new TLegendEntry();
  entry = legend->AddEntry(func2,"Breit-Wigner (#mu = 0, #Gamma=2.5)","l");
  entry = legend->AddEntry(func1,"Gaussian (#mu = 0, #sigma = 2)","l");
  entry = legend->AddEntry(graph4,"Convolution of BW and Gaussian","l");
  entry = legend->AddEntry(graph5,"Overlap Area of Gauss and BW","f");
  
  //Set Legend Properties
  legend->SetFillColor(kWhite);
  entry->SetFillColor(kWhite);
  
  //Create a Canvas to Draw to
  TCanvas *canvas = new TCanvas("canvas","canvas");

  //Draw the Functions and Graphs
  func2->Draw("L");
  graph5->Draw("SAME F");
  func1->Draw("L SAME");
  graph4->Draw("L SAME");
  legend->Draw();

  gPad->Modified(); 
  gPad->Update();

  //Amount by which to move the second function every step
  Double_t moveOver = 0.05;

  //Loop to move the second function toward the first
  for (Int_t i = 0; i < 400; i++) {
    //Set the Mean of the Guassian - This moves the function
    Double_t x = func1->GetParameter(0) + moveOver;
    func1->SetParameter(0, x);
    graph4->SetPoint((graph4->GetN()), x, func3->Integral(-10, 10));
    makeGraph5();

    //Update the picture
    gPad->Modified(); 
    gPad->Update();
 
    //Create an animated gif
    if (createGIF){
      canvas->SaveAs("Gaus_BW_Convolution.gif+50");
    }

  }
}
//This can be run in ROOT by doing .x convolution.C[++]

#include <TMath.h>
#include <TStyle.h>
#include <TF1.h>
#include <TGraph.h>
#include <TLegend.h>
#include <TCanvas.h>

//Breit-Wigner Properties
Double_t bwMean = 0.0;
Double_t bwGamma = 2.5;

//Gaussian Properties
Double_t gWidth = 2.0;

//Set Precision of the Integrals (number of steps)
Int_t precision = 200;

//Create an animated gif?
Bool_t createGIF = kFALSE; // kFALSE or kTRUE

//Instantiate the functions and graphs to be used
//Their properties will be set later
TF1 *func1;     //The Function to be reflected ("stationary")
TF1 *func2;     //The function to be translated ("moving")
TF1 *func3;     //The product of the two functions (func1 and func2)
TF1 *func4;     //The overlap of the two functions (func1 and func2)
TGraph *graph4; //The Resulting Convolution Graph

//Does the product of the two functions (func1 and func2)
Double_t multipliedFuncs(Double_t *x, Double_t *par) {
#if !defined(__CINT__)
  par = par; // get rid of compiler's warning: unused parameter ‘par’
#endif /* !defined(__CINT__) */
  Double_t xx = x[0];
  return func1->Eval(xx) * func2->Eval(xx);
}//End multipliedFuncs

//Find the smaller value of the functions (func1 and func2)
//This is used to shade the overlap region
Double_t smallerValue(Double_t *x, Double_t *par) {
#if !defined(__CINT__)
  par = par; // get rid of compiler's warning: unused parameter ‘par’
#endif /* !defined(__CINT__) */
  Double_t xx = x[0];
  return TMath::Min(func1->Eval(xx), func2->Eval(xx));
}//End smallerValue

void convolution() {
  //Define the Breit-Wigner (func1)
  func1 = new TF1("BreitWigner","TMath::BreitWigner(x,[0],[1])",-10,10);
  func1->SetTitle("Convolution of a Breit-Wigner and a Gaussian");
  func1->SetParameters(bwMean,bwGamma);
  func1->SetNpx(precision);
  func1->SetLineColor(kBlue);
  func1->SetLineWidth(3);
  
  //Create the Gaussian (func2)
  func2 = new TF1("Gaussian","TMath::Gaus(x,[0],[1],1)",-10,10);
  func2->SetParameters(-10,gWidth);
  func2->SetNpx(precision);
  func2->SetLineColor(kBlack);
  func2->SetLineWidth(3);
  
  //Define the Product of the two functions (func1 and func2)
  func3 = new TF1("MultipliedFuncs",multipliedFuncs,-10,10,0);
  func3->SetNpx(precision);
  func3->SetLineColor(kOrange);
  func3->SetLineWidth(3);
  func3->SetFillColor(kOrange);
  func3->SetFillStyle(1001);
  
  //Define the Overlap of the two functions (func1 and func2)
  func4 = new TF1("OverlapFuncs",smallerValue,-10,10,0);
  func4->SetNpx(precision);
  func4->SetLineColor(kYellow);
  func4->SetLineWidth(3);
  func4->SetFillColor(kYellow);
  func4->SetFillStyle(1001);
  
  //Create Convolution Graph
  graph4 = new TGraph();
  graph4->SetName("Convolution");
  graph4->SetLineColor(kRed);
  graph4->SetLineWidth(3);
  //Just a trick ... we will draw a line, so we need at least two points ...
  // graph4->SetPoint((graph4->GetN()), -10, (0.9 * func3->Integral(-10, 10)));
  graph4->SetPoint((graph4->GetN()), -10, func3->Integral(-10, 10));
  
  //Create a Legend
  TLegend *legend = new TLegend(0.55, 0.72, 0.88, 0.88, "", "brNDC");
  legend->SetFillColor(gStyle->GetCanvasColor());
  legend->AddEntry(func1,
                   TString::Format("Breit-Wigner ( #mu = %.1f, #Gamma = %.1f )", bwMean, bwGamma),
                   "l");
  legend->AddEntry(func2,
                   TString::Format("Gaussian ( #sigma = %.1f )", gWidth),
                   "l");
  legend->AddEntry(graph4, "Convolution of BW and Gaussian", "l");
  legend->AddEntry(func4, "Overlap Area of BW and Gaussian", "f");
  legend->AddEntry(func3, "Product of BW and Gaussian", "f");
  
  //Create a Canvas to Draw to
  TCanvas *canvas = new TCanvas("canvas","canvas");
  
  //Draw the Functions and Graphs, Legend
  func1->Draw("L");
  func4->Draw("L SAME");
  func3->Draw("L SAME");
  func1->Draw("L SAME");
  func2->Draw("L SAME");
  graph4->Draw("L");
  legend->Draw();
  gPad->RedrawAxis(); // Redraw axes
  // gPad->RedrawAxis("g"); // Redraw gridx and/or gridy
  gPad->Modified(); gPad->Update();
  
  //Create an animated gif
  if (createGIF){
    canvas->SaveAs("BW_Gaus_Convolution.gif");
  }
  
  //Amount by which to move the second function every step
  Double_t moveOver = 0.05;
  
  //Loop to move the second function toward the first
  for (Int_t i = 0; i < 400; i++) {
    //Set the Mean of the Guassian - This moves the function
    Double_t xx = func2->GetParameter(0) + moveOver;
    func2->SetParameter(0, xx);
    graph4->SetPoint((graph4->GetN()), xx, func3->Integral(-10, 10));
    
    //Update the picture
    gPad->Modified(); gPad->Update();
    
    //Update an animated gif
    if (createGIF){
      canvas->SaveAs("BW_Gaus_Convolution.gif+50");
    }
  }
}

BTW. In your original post, you were explicitly asking to draw the “product” of both functions as the “yellow area”, not the “overlap” (and so, the “yellow area” was correct, but it was simply the “product”, not the “overlap”).

I think we are having notation problems.

In my original post I was writing in the notation of the wiki article. I was complaining that I couldn’t get the resulting function of the convolution to be drawn in real time. In the wiki article the resulting function is written in the notation (f*g)(t). Where as the product of f and g would be written as f(t)*g(t) or in this case since one of the functions is offset f(t)*g(t-t). The convolution is therefore defined as (f*g)(t) = Integral[f(t)*g(t-t), {t`, -Infinity, Infinity}].

Nevertheless, in my original code I was shading the area under the product of f(t)*g(t-t`) which gives the overlap area for two square waves of unit height and width…but not much else :wink: !