I am wondering what is the correct function in Mathematica to plot the true impulse function, better known as the DiracDelta[] function. When using this inside of a function or just the function itself when plotting, it renders output = zero. Quick example:

Plot[DiracDelta[x], {x,-1,1}]

I am wondering, is this the correct delta function which is infinite in height at zero and zero everywhere else. I have seen other functions such as KroneckerDelta[], but this seems to do the same exact thing. The code is below:

eqn1wb := (2 ((I Pi f) Exp[(-I) 2 Pi f]) ((1/2) DiracDelta[f-2] + 
          DiracDelta[f+2]))/(1 + I 2 Pi f); 
wb1 = Plot[Re[eqn1wb], {f, -5, 5}, 
PlotRange :> All, PlotStyle :> {Thick, Red},
AspectRatio :> 1,
GridLines :> {{-4, -2, 2, 4}, {1.0, 0.5, -0.5, -1.0}}, 
PlotLabel :> "Re(w) - Frequency Response", 
AxesLabel :> {"f", "W(f)"}

5 Answers 5


To create the plot you could replace any occurrence of DiracDelta[a] with something like 10000 UnitStep[1/10000 - a^2]], so for example to plot

f[x_] := DiracDelta[x - 2] + DiracDelta[x + 2]

you could do something like

Plot[Evaluate[f[x] /. DiracDelta[a_] :> 10000 UnitStep[1/10000 - a^2]], 
  {x, -4, 4}, Exclusions -> None, PlotPoints -> 800]

Mathematica graphics

Note that for Mathematica to see the discontinuities you need to increase the number of plot points. The number of points needed will depend on the plot range, so you might have to tweak that.

  • 5
    $\begingroup$ A large number of PlotPoints isn't necessary. There is a nice possibility to specify additional points with PlotPoints->{Automatic, {-2, 2}} (ref). Then recursion algorithm will plot it fine. $\endgroup$
    – ybeltukov
    Commented Oct 28, 2015 at 19:13

As others already have written, the Dirac delta is not a real function and it can't be plotted. Other programs that claim to plot it just fake it.

Having said that, you can roll a diracDelta of your own, that more or less mimics the Dirac Delta's behavior but is still continuous. Advantage with respect to Heike's solution is that we don't need to increase the number of PlotPoints. I'll use the NormalDistrution with a very small standard deviation:

diracDelta[x_] = PDF[NormalDistribution[0, 1/100], x];

eqn1wb := (2 ((I Pi f) Exp[(-I) 2 Pi f]) ((1/2) DiracDelta[f - 2] + 
       DiracDelta[f + 2]))/(1 + I 2 Pi f);

wb1 = Plot[
  Re[eqn1wb] /. DiracDelta -> diracDelta // Evaluate, {f, -5, 5}, 
  PlotStyle :> {Thick, Red}, AspectRatio :> 1, 
  GridLines :> {{-4, -2, 2, 4}, {1.0, 0.5, -0.5, -1.0}}, 
  PlotLabel :> "Re(w) - Frequency Response", 
  AxesLabel :> {"f", "W(f)"}, PlotRange -> {Automatic, 1.5}]

Mathematica graphics

DiracDelta's are sometimes plotted by using arrows. Of course we can do that as well. In this case we can set the function at zero everywhere and add arraows in the Epilog part of the plot.

 Plot[0, {f, -5, 5}, PlotStyle :> {Thick, Red}, AspectRatio :> 1,
   GridLines :> {{-4, -2, 2, 4}, {1.0, 0.5, -0.5, -1.0}}, 
  PlotLabel :> "Re(w) - Frequency Response", 
  AxesLabel :> {"f", "W(f)"}, PlotRange -> {Automatic, 1.}, 
  Epilog -> {Red, Thick,Arrow[{{-2, 0}, {-2, 1}}], Arrow[{{2, 0}, {2, 1}}]}]

Mathematica graphics

I believe the convention is to scale the lengths of the arrows in proportion to the factors in front of the DiracDelta. The above can be easily extended to do that.

Factors can be found using Coefficient:

Coefficient[eqn1wb, DiracDelta[-2 + f]]

==> (I E^(-2 I f \[Pi]) f \[Pi])/(1 + 2 I f \[Pi])

Ratio's of coefficients:

Coefficient[eqn1wb, DiracDelta[2 + f]]/Coefficient[eqn1wb, DiracDelta[2 - f]]

==> 2

Locations of the Dirac delta's:

Cases[eqn1wb, DiracDelta[a__] :> Solve[a == 0, f], Infinity]

==> {{{f -> 2}}, {{f -> -2}}}

I'm not going to automate it all, but the idea is clear. In the end you'll get something like:

Mathematica graphics


The Dirac delta, $\delta(x)$ is zero everywhere except at zero, and has an integral of 1 over $\mathbb{R}$. It is not really a function in the true sense and equating $\delta(0)=\infty$ is a rather loose definition; it should technically be considered as a distribution or a delta measure. Mathematica's implementation of DiracDelta remains unevaluated at 0 and has all the other properties.

DiracDelta /@ {-1, 0, 1}
Out[1]= {0, DiracDelta[0], 0}

Integrate[DiracDelta[x] , {x, -Infinity, Infinity}]
Out[2]= 1

This explains why you weren't getting any result with Plot.

In the off-chance that you were trying to plot the impulse response of a discrete-time system, then you might be interested in the DiscreteDelta function:

ListPlot[DiscreteDelta /@ Range[-5, 5], Filling -> Axis, PlotMarkers -> {Automatic, 10}]

enter image description here

You can also reproduce your MATLAB "plot" of the Dirac delta function in Mathematica. Let me first note that you're not using their dirac() function to produce the output you've shown below. In fact, they define it "loosely":

function Y = dirac(X)
%DIRAC  Delta function.
%    DIRAC(X) is zero for all X, except X == 0 where it is infinite.
%    more comments

Y = zeros(size(X));
Y(X == 0) = Inf;

and this would not have given you the plot you showed, because of the Inf. You have probably replaced Inf with 10000 or written a similar function. Heike, image_doctor and Sjoerd have shown you ways of plotting it. You can also roll your own, a la MATLAB style as:

dirac[x_] := Piecewise[{{10000, x == 0}, {0, True}}]
ListLinePlot[{#, dirac@#} & /@ Range[-1, 1, 1/1000], Frame -> True, Axes -> False]

enter image description here

Now use this definition to redefine DiracDelta inside a Block and you can use your original code with some slight modifications:

wb1 = Block[{DiracDelta = Piecewise[{{10000, # == 0}, {0, True}}] &}, 
  ListLinePlot[Re[eqn1wb] /. f -> # & /@ Range[-5, 5, 1/100], 
   DataRange -> {-5, 5}, PlotRange :> {All, {0, 1}}, 
   PlotStyle :> {Thick, Red}, AspectRatio :> 1, 
   GridLines :> {{-4, -2, 2, 4}, {1.0, 0.5, -0.5, -1.0}}, 
   PlotLabel :> "Re(w) - Frequency Response", 
   AxesLabel :> {"f", "W(f)"}]

enter image description here

If I were to use a schematic to show the delta function, I would do it similar to Sjoerd's arrow plots.

  • $\begingroup$ Okay, so what your saying is, is that it cannot be plotted over $\mathbb{R}$? I can accomplish it in another program, but wanted to do it in Mathematica. $\endgroup$
    – night owl
    Commented Mar 26, 2012 at 8:26
  • $\begingroup$ Can you show your output plot that you generated in your other program? $\endgroup$
    – rm -rf
    Commented Mar 26, 2012 at 8:30
  • $\begingroup$ You can define your own Dirac Delta function which has the behaviour that you would like. DiracD[x_] := DiracDelta[x] /. DiracDelta[0] -> Infinity . You can swap infinity for any large value that suits your application if you wish. $\endgroup$ Commented Mar 26, 2012 at 8:36
  • $\begingroup$ @R.M: Here is the output graph of the Delta function. i.sstatic.net/VSEr9.png $\endgroup$
    – night owl
    Commented Mar 26, 2012 at 9:00
  • $\begingroup$ @image_doctor: This seems to not work. I used Plot[DiracD[x], {x,-1,1}, PlotStyle:> Thick] to get this: i.sstatic.net/Yh6r8.png $\endgroup$
    – night owl
    Commented Mar 26, 2012 at 9:03

Building on Sjoerd C. de Vries's solution, I wrote this little function to make the task easy for anyone interested in plotting Dirac's Delta as arrows in a Plot command in Mathematica. The function receives an equation and a variable, and returns a list of {Arrow[]}'s in the correct location of the Diract Deltas, that can be included in the Epilog of a Plot command.

ArrowsDeltaFunction[eqn_, x_] := Module[{xsubs, listDeltas, coefDeltas, locationDeltas},
  xsubs = (x /. Cases[eqn, DiracDelta[a__] :> Solve[a == 0, x], Infinity]) /. x -> {};
  listDeltas = DiracDelta[x - x0] /. x0 -> xsubs;
  coefDeltas = Flatten[Coefficient[eqn, listDeltas]];
  locationDeltas = Flatten[xsubs];

  Arrow[Table[{{locationDeltas[[i]], 0}, {locationDeltas[[i]], 
      coefDeltas[[i]]}}, {i, 1, Length[locationDeltas]}]]

So you can use this function throughout your code by including it in the Plot Epilog, like this:

eqn1 := x + (3/4) DiracDelta[x - 0.2] + (\[Pi]/3) DiracDelta[x + 0.5]

Plot[eqn1, {x, -5, 5},
   PlotStyle -> {Thick, Red},
   Epilog -> {Thick, Red, ArrowsDeltaFunction[eqn1, x]}

Output of Plot Command

You can see the arrows appear in the "gaps" of the function where the Delta's exist, and their height correspond to the coefficient of the Delta. The code works even if the function does not have Deltas (no arrows are plotted).

Edit: This can have nice applications for teachers to explain Fourier Transforms, its properties, etc.

(* Time Domain Signal (in this case complex) *)
Ft :=  (1 - I/2) Exp[-I t] +  Cos[2 t]  + 2 Sin[4 t]

(* Fourier Transform *)
Fw := FourierTransform[Ft, t, w] 
ReFw := ComplexExpand[Re[Fw]]
ImFw := ComplexExpand[Im[Fw]]

(* Plot *)
Plot[{ReFw, ImFw}, {w, -5, 5},
 PlotStyle -> {Thick, Red, Blue},
 PlotRange -> 3,
 AxesLabel -> {"w", "F(w)"},
 PlotLegends -> LineLegend[{Red, Blue}, {"Real part", "Imag part"}],
 Epilog -> {Thick, 
    Red, ArrowsDeltaFunction[ReFw, w], 
    Blue, ArrowsDeltaFunction[ImFw, w] 

Fourier Transform Examples

  • $\begingroup$ In the case of functions like eqn=DiracDelta[2x-2], the ArrowsDeltaFunction function attempts to extract the coefficient of DiracDelta[x-1], instead, even though no such term exists. This results in some arrows incorrectly having 0 height. I implemented a fix, which I'll post below. $\endgroup$ Commented May 18, 2016 at 10:11
  • 1
    $\begingroup$ ArrowsDeltaFunction[eqn_, x_] := Module[{xsubs, listDeltas, coefDeltas, locationDeltas}, xsubs = (x /. Cases[eqn, DiracDelta[a__] :> Solve[a == 0, x], Infinity]) /. x -> {}; locationDeltas = Flatten[xsubs]; Arrow[Table[{{locationDeltas[[i]], 0}, {locationDeltas[[i]], (eqn/ DiracDelta[x - locationDeltas[[i]]]) /. x -> locationDeltas[[i]]}}, {i, 1, Length[locationDeltas]}]]] $\endgroup$ Commented May 18, 2016 at 10:12

At your own risk, you can redefine the DiracDelta function to have the behaviour you want:

DiracDelta[0] := your large number;

The integral given by

Integrate[DiracDelta[x], {x, -Infinity, Infinity}]

still seems to equate to one with the new definition in place.

  • $\begingroup$ Thanks, but I am not really interested in a numerical value for it. I am looking to plot the function. Am I not seeing how to do this?, because I am not getting the correct output using this. $\endgroup$
    – night owl
    Commented Mar 26, 2012 at 9:21
  • $\begingroup$ h[x_] := With[{delta = 0.01, maxVal = 5}, Piecewise[{{maxVal, x <= delta && x >= -delta}, {0, x < -delta}, {0, x > delta}}, x]] and Plot[h[x], {x, -10, 10}, PlotRange -> {{-10, 10}, {0, 10}}, Mesh -> All, PlotPoints -> 5000] Shows you which points plot is trying to display. $\endgroup$ Commented Mar 26, 2012 at 9:59
  • $\begingroup$ Thanks, But seems too difficult for mathematica to deal with, which has surprised me. If I have a function that I have written in matlab for the delta function, can I use that to call it in mathematica, or save it as a .m function file? $\endgroup$
    – night owl
    Commented Mar 26, 2012 at 10:08
  • $\begingroup$ Do you mean a function in Matlab to create the plot you want? If so, Mathematica can execute any operating system command line using Run[]. So using Matlab in batch mode to make the graph, you can import the graph into Mathematica from a saved matlab image file. If you mean using Matlab to compute the Dirac funtion .. I think Mathematica should be perfectly able to duplicate the matlab function. Hope that helps :) $\endgroup$ Commented Mar 26, 2012 at 10:30
  • $\begingroup$ Was there something in the examples discussed above that doesn't replicate the mathematical behaviour of the function you want? $\endgroup$ Commented Mar 26, 2012 at 10:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.