% get_Fourier_transform.m

%

% [ft_power,ft,omega,iflag] = ...

%    get_Fourier_transform( ...

%        dt,t0,f_val,iplot,ipad,omega_plot_max, ...

%        filter_factor);

%

% This MATLAB function calculates and plots

% (optionally) the Fourier transform data

% of an input signal.

%

% dt = time sampling interval

% t0 = initial time of signal data

% f_val = vector of f(t) values at sample times

% iplot = if non-zero, plot results

% ipad = if non-zero, pad data vector as necessary

%           to make number of samples a power of 2

% omega_plot_max = if 0, plot entire frequency range,

%           if non-zero, denotes upper frequency for plot

% filter_factor = if not supplied or zero, no filtering

%           of the data is performed.  Otherwise, the

%           routine applies a digital bandpass filter

%           with an upper cutoff of filter_factor times

%           the Nyquist critical sampling fequency, pi/dt.

%           A value of 0.5 should guard stringently against

%           aliasing.

%

% K. Beers. MIT ChE. 8/13/2002

function [ft_power,ft,omega,iflag] = ...

    get_Fourier_transform( ...

        dt,t0,f_val,iplot,ipad,omega_plot_max, ...

        filter_factor);

iflag = 0;

% Get the length of the input signal data vector.

N = length(f_val);

% make sure that f_val is a row vector

if(size(f_val,1)~=1)

    f_val = f_val';

end

% If requested through the optional argument

% filter_factor, a bandpass filter is applied

% to the data to remove any high-frequency

% signals that may cause aliasing.

% The filter coefficients are taken from

% Press et al., Numerical Recipies in C,

% 2nd ed., 1992, p. 562.

if(exist('filter_factor'))

if(filter_factor~=0)

    P = N*dt/2;

    omega_filter_min = pi/(4*P);

    omega_crit = pi/dt;

    omega_filter_max = filter_factor*omega_crit;

    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_val = -filter(b_filter,a_filter,f_val);

end

end

% We now remove any linear trend in the data.

% This is done because the sample measurements

% might have a long-term drift.

f_val = detrend(f_val);

% if N is not even, add a zero at the end

% to ensure that it is

if(rem(N,2))

    N = N + 1;

    f_val = [f_val 0];

end

% If ipad is non-zero, we round up the number of

% points to the next highest power of 2 by adding

% zeros.

if(ipad)

    two_power = log2(N);

    if(rem(two_power,1))  % two_power is not an integer

        N = pow2(floor(two_power)+ipad);

    end

    num_pad = N - length(f_val);

    f_val = [f_val zeros(1,num_pad)];

end

% calculate the sampled angular frequencies.

P = N*dt/2;  % assummed half-period of input signal

omega_Nyquist = pi/dt;  % max. sampled frequency

d_omega = pi/P;  % difference between successive frequencies

omega_max = (N-1)*d_omega;

omega = [0 : d_omega : omega_max];

if(omega_plot_max==0)

    omega_plot_max = omega_max;

end

% We invoke the MATLAB function fft() to calculate

% to calculate the Fourier transform coefficients.

if(ipad)

    ft = fft(f_val,N);

else

    ft = fft(f_val);

end

% correct for time step

ft = (dt/sqrt(2*pi)).*ft;

% we compute now the power spectrum

ft_power = conj(ft).*ft;

% plot the time signal and Fourier series

if(iplot)

    % if a filter was used, we only report at most the

    % data in the frequency range below the cutoff.

    if(exist('filter_factor'))

        if(omega_plot_max > omega_filter_max)

            omega_plot_max = omega_filter_max;

        end

    end

    j_plot = find(omega <= omega_plot_max);

    figure;

    % plot input signal

    subplot(2,2,1);

    t_end = t0 + (N-1)*dt;

    plot([t0:dt:t_end],f_val);

    xlabel('t');

    ylabel('f(t)');

    title('Input signal');

    axis tight;

    % real part of Fourier transform

    subplot(2,2,3);

    plot(omega(j_plot),real(ft(j_plot)));

    xlabel('\omega (radians/time)');

    ylabel('Real part of FFT');

    axis tight;

    % imaginary part

    subplot(2,2,4);

    plot(omega(j_plot),imag(ft(j_plot)));

    xlabel('\omega (radians/time)');

    ylabel('Imaginary part of FFT');

    axis tight;

    % power spectrum

    subplot(2,2,2);

    plot(omega(j_plot),ft_power(j_plot));

    xlabel('\omega (radians/time)');

    ylabel('|F(\omega)|^2');

    title('Power spectrum');

    axis tight;

end

iflag = 1;

return;