fftdemo_hrf.m

%% FFTDEMO_HRF (adapted from FFTDEMO in Matlab)
% FFT for Spectral Analysis
% This example shows the use of the FFT function for spectral  analysis.  A
% common use of FFT's is to find the frequency components of a signal buried in
% a noisy time domain signal.
%
% Copyright 1984-2005 The MathWorks, Inc. 
% $Revision: 5.8.4.2 $  $Date: 2009/04/03 21:23:24 $

function fftdemo_hrf()

%% Data
% First create some data.  Consider data sampled at 1/1.5 Hz (TR = 1.5s).  Start by forming a
% time axis for our data, running from t=0 until t=250 s in steps of 1.5s
% .  Then form a signal, x, containing sine waves at 0.125 and 0.375
% Hz.

t = 0:1.5:250;
f = 0.02; % in Hz
x = sin(2*pi*f*t) + sin(2*pi*3*f*t);

noiseLevel = 1; % low = 0.1; high = 1.0

%% Noise
% Add some random noise with a standard deviation of 2 to  produce a noisy
% signal y.  Take a look at this noisy  signal y by plotting it.

f_ = figure(1)
%  = figure(1);
set(f_, 'position', [1873         133         560         420])

clf
subplot(3,1,1)
y = x + noiseLevel*randn(size(t));
% plot(t(1:50), y(1:50),'.-')
plot(t, y,'.-'), hold on
plot(t, x,'k-')

title(sprintf('Noisy time domain signal of %.2f and 3*%.2f Hz + some noise', f, f))
ylabel('Signal')
xlabel('Time')

%%
% Clearly, it is difficult to identify the frequency components from looking at
% this signal; that's why spectral analysis is so popular.
%
% Finding the discrete Fourier transform of the noisy signal y is easy; just
% take the fast-Fourier transform (FFT).

Y = fft(y,251);
X = fft(x,251); % clean data


%%
% Compute the power spectral density, a measurement of the energy at various
% frequencies, using the complex conjugate (CONJ).  Form a  frequency axis for
% the first 127 points and use it to plot the result.  (The remainder of the
% points are symmetric.)

subplot(3,1,2)
Pyy = Y.*conj(Y)/251;
Pxx = X.*conj(X)/251;
f = (1/1.5)/251*(0:127); % sampling rate / # points for 1/2 the sequence
plot(f,Pxx(1:128)); hold on
plot(f,Pyy(1:128),'r','linewidth',2)
title('Power spectral density')
xlabel('Frequency (Hz)')
ylabel('Power')

%%
% Zoom in and plot only 1/5 the original range of 1/1.5 Hz.  Notice the peaks at 50 Hz and  120 Hz.
% These are the frequencies of the original signal.

subplot(3,1,3)
plot(f(1:50),Pxx(1:50),'k-'); hold on
plot(f(1:50),Pyy(1:50),'r')
title('Power spectral density (Zommed X-axis)')
xlabel('Frequency (Hz)')
ylabel('Power')

% suptitle('sin(2\pi f t) + sin(2\pi 3 f t) + noise')