Commit 7a0eea82 authored by Sam Harrison's avatar Sam Harrison
Browse files

Replace `spm_hrf` with FLOBS

Removes another hidden dependency.
parent ea792f87
......@@ -52,9 +52,32 @@ end
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ HRF, t_HRF ] = loadFLOBS(dt_HRF)
% Loads the FLOBS HRF basis, and resamples to `dt_HRF`
% Optionally returns the sample times, `t_HRF`
% Load the basis
fsldir = getenv('FSLDIR');
assert(~isempty(fsldir), 'Could not find `$FSLDIR`!');
hrf = load(fullfile(fsldir, 'etc', 'default_flobs.flobs', 'hrfbasisfns.txt'));
dt_hrf = 0.05; %FLOBS sampling rate
% Establish time samplings
n_hrf = length(hrf);
t_hrf = linspace(0.0, dt_hrf * (n_hrf - 1), n_hrf);
n_HRF = floor(t_hrf(end) / dt_HRF);
t_HRF = linspace(0.0, dt_HRF * (n_HRF - 1), n_HRF);
% Interpolate HRF to requested rate
HRF = interp1(t_hrf, hrf, t_HRF, 'pchip', 0.0);
end
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ PA, A ] = generateBoldSignal_LinearHRF(P, An, params, options, plotFigures)
%Convolves timecourses with the SPM linear HRF and mixes with spatial maps
%to get BOLD signal
%Convolves timecourses with the mean FLOBS linear HRF and mixes with spatial
%maps to get BOLD signal
%
% options.BS.SNR_P - SNR (in terms of power) of the maps
% options.BS.SNR_A - SNR (in terms of power) of the timecourses
......@@ -63,7 +86,7 @@ function [ PA, A ] = generateBoldSignal_LinearHRF(P, An, params, options, plotFi
%voxelwise timecourses
%Load the HRF
hrf = spm_hrf(params.dt);
hrf = loadFLOBS(params.dt); hrf = hrf(:,1);
PA = cell(params.S,1);
A = cell(params.S,1);
......@@ -118,17 +141,11 @@ function [ PA, A ] = generateBoldSignal_FlobsHRF(P, An, params, options, plotFig
%voxelwise timecourses
%Load the HRF
hrf = load([getenv('FSLDIR') '/etc/default_flobs.flobs/hrfbasisfns.txt']);
dt = 0.05; %FLOBS sampling rate
%Interpolate HRF to requested rate
tHRF = dt*(1:length(hrf)) - dt;
tNew = params.dt * (0:floor(tHRF(end)/params.dt));
hrf = interp1(tHRF, hrf, tNew, 'pchip', 0);
[hrf, t_hrf] = loadFLOBS(params.dt);
if plotFigures
figure; hold all; nPlots = 0;
xlim([tNew(1) tNew(end)]); xlabel('Time'); ylabel('HRF')
xlim([t_hrf(1) t_hrf(end)]); xlabel('Time'); ylabel('HRF')
title('FLOBS HRF: example HRFs')
end
......@@ -147,7 +164,7 @@ for s = 1:params.S
end
end
if plotFigures && (nPlots<100)
plot(tNew, randHrf); nPlots = nPlots + params.N;
plot(t_hrf, randHrf); nPlots = nPlots + params.N;
end
PA{s} = cell(params.R(s),1);
......@@ -286,7 +303,7 @@ NFFT = params.T;
Fs = 1/params.TR; %sampling frequency
f = Fs/2*linspace(0,1,NFFT/2+1);
%Example HRF
hrf = spm_hrf(params.TR);
hrf = loadFLOBS(params.TR); hrf = hrf(:,1);
fHRF = abs(fft(hrf, NFFT))';
% mean frequency content for each scan
figure; hold on
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment