Frequency response function
- To: mathgroup at smc.vnet.net
- Subject: [mg82755] Frequency response function
- From: Will Robertson <wspr81 at gmail.com>
- Date: Tue, 30 Oct 2007 03:36:29 -0500 (EST)
Hello,
People have asked in the past for how to plot a power spectral density
plot, or a frequency response plot, or whatever you like to call it.
Since there isn't an archive of a good solution on this newgroup,
here's a quick and basic solution that I just put together:
FRF[signal_, tt : {t_, tmin_, tmax_, tstep_}] :=
Module[{s, list, fft, N, L, freq},
(* Generate the signal and its DFT: *)
list = Table[signal, tt];
fft = Abs[Chop@Fourier@list]^2;
(* Get the number of samples and throw away the wraparound: *)
N = Length[list];
L = If[OddQ[N], (N - 1)/2, N/2];
fft = Take[fft, L];
(* Scale the frequencies and plot the FRF: *)
freq = Range[L]/(N tstep);
ListLogPlot[Transpose@{freq, fft}, PlotRange -> All,
Joined -> True]
]
With[{\[Omega] = 2},
FRF[Sin[2 \[Pi] \[Omega] t] + Sin[2 \[Pi] 10 \[Omega] t], {t, 0, 10,
1/50}]]
(Modifiable & distributable under the Apache License v2.)
Obviously this is just a "first pass" at the many many features that
should be added to such a function (such as padding, windowing,
support for complex data even!). I don't even know if it works
correctly for tmin != 0.
I'd be great for people to post their own improvements in this thread.
Eventually I might make a package out of the whole thing. (A port of
the spectral density routines here would be ideally the best: <http://
www.mecheng.adelaide.edu.au/~pvl/octave/>)
Cheers,
Will