% results_get_Fourier_transform.m
%
% This text file contains the results of
% calculations using get_Fourier_transform.
Example 1. Simple linear combination of
functions
-------------------------------------------------
% Set input signal, parameters
>> t_val = linspace(0,6*pi,2^6);
>> dt = t_val(2) - t_val(1);
>> t0 = t_val(1);
>> f_val = sin(t_val) + 2*cos(2*t_val);
>> iplot = 1;
>> ipad = 0;
>> omega_plot_max = 0;
% call routine
>> [ft_power,ft,omega,iflag] = ...
get_Fourier_transform(
...
dt,t0,f_val,iplot,ipad,omega_plot_max);
% results are in plot : FFT_ex1_1.fig
% Then, we plot again only up to a frequency
of 4
>> omega_plot_max = 4;
>> [ft_power,ft,omega,iflag] = ...
get_Fourier_transform(
...
dt,t0,f_val,iplot,ipad,omega_plot_max);
% results are in plot : FFT_ex1_2.fig
% We now use the same data, but sample more
frequently.
>> t_val = linspace(0,6*pi,2^10);
>> dt = t_val(2) - t_val(1);
>> t0 = t_val(1);
>> f_val = sin(t_val) + 2*cos(2*t_val);
>> iplot = 1;
>> ipad = 0;
>> omega_plot_max = 0;
% call routine
>> [ft_power,ft,omega,iflag] = ...
get_Fourier_transform(
...
dt,t0,f_val,iplot,ipad,omega_plot_max);
% results are in plot : FFT_ex1_3.fig
% We now sample for a longer period of time
to decrease the
% spacing between sampled frequencies.
We also pad with zeros
% at the end of the data sample to bring the
number of points
% up to a power of 2. This is done,
because the FFT algorithm
% is more effective when the number of
points is a power of 2.
>> t_val = linspace(0,48*pi,10000);
>> dt = t_val(2) - t_val(1);
>> t0 = t_val(1);
>> f_val = sin(t_val) + 2*cos(2*t_val);
>> iplot = 1;
>> ipad = 1;
>> omega_plot_max = 4;
% call routine
>> [ft_power,ft,omega,iflag] = ...
get_Fourier_transform(
...
dt,t0,f_val,iplot,ipad,omega_plot_max);
% results are in plot : FFT_ex1_4.fig
Example 2. Signal with random noise
-----------------------------------
% Generate signal with random Gaussian noise
>> num_points = 2^17;
>> t_val = linspace(0,4*pi,num_points);
>> dt = t_val(2) - t_val(1);
>> t0 = t_val(1);
>> f_val = zeros(size(t_val));
>> for k=1:10
f_val = f_val + (1/sqrt(k))*sin(k*t_val);
end
>> noise_mag = 0.2;
>> f_val = f_val + noise_mag*randn(size(f_val));
>> iplot = 1;
>> ipad = 1;
>> omega_plot_max = 15;
% call routine
>> [ft_power,ft,omega,iflag] = ...
get_Fourier_transform(
...
dt,t0,f_val,iplot,ipad,omega_plot_max);
% results are in plot : FFT_ex2_1.fig
% We now attempt to remove the noise by
applying
% a digital filter. The following
filter is
% designed to remove all contributions to
the
% signal with frequencies outside of the
% range [omega_filter_max,omega_filter_min].
% This is known as a bandpass filter.
% The filter coefficients are taken from
% Press et al., Numerical Recipies in C,
% 2nd ed., 1992, p. 562.
>> omega_filter_min = 1e-5;
>> omega_filter_max = 50;
>> w_min = tan(pi*omega_filter_min/(2*pi)*dt);
>> w_max = tan(pi*omega_filter_max/(2*pi)*dt);
>> a_filter = zeros(1,3);
>> a_filter(1) = 1;
>> a_filter(2) = ...
-( (1+w_min)*(1-w_max) +
(1-w_min)*(1+w_max) ) ...
/ (1+w_min) / (1 + w_max);
>> a_filter(3) = (1-w_min)*(1-w_max)
...
/ (1+w_min)/(1+w_max);
>> b_filter = zeros(1,3);
>> b_filter(1) = -w_max/(1+w_min)/(1+w_max);
>> b_filter(3) = w_max/(1+w_min)/(1+w_max);
>> f_filtered = -filter(b_filter,a_filter,f_val);
% make a plot of the original and filtered
data
>> figure;
>> subplot(2,1,1);
>> plot(t_val,f_val);
>> xlabel('t');
>> ylabel('f(t)');
>> title(['Bandpass filter applied to
f(t), ', ...
'[\omega_{min},\omega_{max}]
= ', ...
'[
', num2str(omega_filter_min), ' , ', ...
num2str(omega_filter_max),
']']);
>> axis tight;
>> subplot(2,1,2);
>> plot(t_val,f_filtered);
>> xlabel('t');
>> ylabel('f(t) (filtered)');
>> axis tight;
% call FFT on filtered data
>> [ft_power,ft,omega,iflag] = ...
get_Fourier_transform(
...
dt,t0,f_filtered,iplot,ipad,omega_plot_max);
% results are in plots : FFT_ex2_2.fig,
FFT_ex2_3.fig
Example 3. FFT with aliasing
----------------------------
% We now consider a problem with aliasing,
the corruption of
% Fourier transform data due to non-sampled
higher-frequency
% contributions to f(t).
%
% We choose a problem with a half-period P =
3*pi with
% sampling of N = 2^M points.
>> P = 3*pi;
>> M = 6;
>> N = pow2(6);
>> dt = 2*P/N;
% The critical frequency, the maximum value
sampled by
% this data, is
>> omega_crit = pi/dt;
% We now generate a signal with a
contribution that is
% outside of the range of the sampled
frequencies.
>> t_val = linspace(0,2*P,N);
>> alpha = 1.1;
>> omega_high = alpha*omega_crit;
>> f_val = sin(t_val) + 2*cos(2*t_val)
+ ...
sin(omega_high*t_val);
% We now generate the FFT for this problem.
>> t0 = t_val(1);
>> iplot = 1;
>> ipad = 0;
>> omega_plot_max = 0;
>> [ft_power,ft,omega,iflag] = ...
get_Fourier_transform(
...
dt,t0,f_val,iplot,ipad,omega_plot_max);
% result is in plot : FFT_ex3_1.fig
% We see that the Fourier transform has
peaks at
% omega = 1 and omega = 2 from the sin and
cos terms
% respectively. But, it also has peaks
near the
% middle that arise because of the unsampled
high-frequency
% contribution that is placed, by mistake,
in the sampled
% frequency domain. This is the
phenomenon of aliasing,
% and is to be avoided when computing
Fourier transforms.
% In this case, the contribution to f(t) at
a frequency
% of omega_high = 11.7333 is "wrapped"
back into the
% low-frequency sampled domain and appears
as a spurious
% peak at slightly less than a frequency of
10. Simply
% by looking at the form of f(t) in this
example, we
% know that the real Fourier transform
should have no
% peak here.
% We show how to remove aliasing by two
approaches. The
% first is to adjust the sampling rate so
that any
% remaining high-frequency contributions are
sampled.
% We show that this approach removes
aliasing, by repeating
% the Fourier transform calculation with a
smaller dt.
>> N = pow2(8);
>> dt = 2*P/N;
>> t_val = linspace(0,2*P,N);
>> f_val = sin(t_val) + 2*cos(2*t_val)
+ ...
sin(omega_high*t_val);
>> [ft_power,ft,omega,iflag] = ...
get_Fourier_transform(
...
dt,t0,f_val,iplot,ipad,omega_plot_max);
% result is in plot : FFT_ex3_2.fig
% replot with only interesting region
>> omega_plot_max = 20;
>> [ft_power,ft,omega,iflag] = ...
get_Fourier_transform(
...
dt,t0,f_val,iplot,ipad,omega_plot_max);
% results are in plot : FFT_ex3_3.fig
% Since we may not know apriori how to
select the sampling
% rate to remove the bias, an alternative is
to apply
% a bandpass filter to remove any high-frequency
contributions
% before calculating the Fourier transform.
The routine
% get_Fourier_transform is now called with
an argument that
% calls for the digital filter to be applied
to remove aliasing.
% We now invoke the routine with the digital
filter on the
% original data.
>> P = 3*pi;
>> M = 6;
>> N = pow2(6);
>> dt = 2*P/N;
>> omega_crit = pi/dt;
% We now generate a signal with a
contribution that is
% outside of the range of the sampled
frequencies.
>> t_val = linspace(0,2*P,N);
>> alpha = 1.1;
>> omega_high = alpha*omega_crit;
>> f_val = sin(t_val) + 2*cos(2*t_val)
+ ...
sin(omega_high*t_val);
% We now generate the FFT for this problem.
>> t0 = t_val(1);
>> iplot = 1;
>> ipad = 0;
>> omega_plot_max = 0;
% The following number sets the ratio
between the upper cutoff
% of the bandpass filter and the critical
frequency. If filter_factor
% is 0.5, the bandpass filters out all
contributions above 1/2 of the
% critical frequency, which should remove
all trace of aliasing (while
% maintaining accuracy in left-side (low-frequency)
half ot the spectrum.
>> filter_factor = 0.5;
>> [ft_power,ft,omega,iflag] = ...
get_Fourier_transform(
...
dt,t0,f_val,iplot,ipad,omega_plot_max,
...
filter_factor);
% results are in plot : FFT_ex3_4.fig