Re: Fourier Transform with Differential Equation
- To: mathgroup at smc.vnet.net
- Subject: [mg48978] Re: Fourier Transform with Differential Equation
- From: Paul Abbott <paul at physics.uwa.edu.au>
- Date: Fri, 25 Jun 2004 17:52:20 -0400 (EDT)
- Organization: The University of Western Australia
- References: <cbgkes$cls$1@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
In article <cbgkes$cls$1 at smc.vnet.net>, Lee Fisher <lfis at helix.nih.gov>
wrote:
> I'm attempting to find the frequency response of the output of an RC
> circuit (a low-pass filter). For some reason, the output of the Fourier
> transform in mathematica does not seem to correspond to anything I would
> expect from this differential equation. It seems that the placement of
> the initial conditions has an enormous effect on the shape of the
> frequency response, so much so that the phase portion looks nothing like
> i would expect it to. I'm using the following code
>
> w=1
> result=DSolve[{10 V'[t]+V[t]==E^(I w t), V[0]==0},V,t}]
> time=Table[i,{i,0,100}];
> out=Re[Evaluate[V[time]/.result[[1]]]];
> fout=Fourier[out];
> ListPlot[Abs[fout],PlotJoined->True,PlotRange->All]
> ListPlot[Arg[fout],PlotJoined->True,PlotRange->All]
>
> and adjusting w to observe the response at different frequencies. While
> the magnitude graph is somewhat reminiscent of what I would expect, the
> phase graph looks nothing like I'd expect (it should be at Pi highly
> negative frequencies, negative Pi at highly positive frequencies, and at
> 0 at the zero frequency). Could someone please explain to me either
> where my code is wrong, or where my expectation of the magnitude and
> phase output is wrong.
With w=1, what is the period, and what impact does that have on the
sampling rate? Consider the following modification of your code:
For the angular frequency,
w = 2 Pi 100
compute the period:
period = 2 Pi/w
Solve the differential equation as before:
result=DSolve[{10 V'[t]+V[t]==E^(I w t), V[0]==0},V,t}]
Sample the voltage over one period (not repeating the end-point):
time = Range[period, 1, period];
out = Re[V[time] /. First[result]];
Compute the Fourier transform:
fout = Fourier[out];
Plot the magnitude,
ListPlot[Abs[fout], PlotJoined -> True, PlotRange -> All]
and the phase, rotating the list so that zero frequency appears in the
middle:
ListPlot[RotateRight[Arg[fout], Floor[Length[fout]/2]],
PlotJoined -> True, PlotRange -> All]
Hope this helps.
Cheers,
Paul
--
Paul Abbott Phone: +61 8 9380 2734
School of Physics, M013 Fax: +61 8 9380 1014
The University of Western Australia (CRICOS Provider No 00126G)
35 Stirling Highway
Crawley WA 6009 mailto:paul at physics.uwa.edu.au
AUSTRALIA http://physics.uwa.edu.au/~paul