vtcbandpassfilter.m

VTC file is the file format of the Brainvoyager. The FMRI matchine scan human s brain and get the vt

function [addTo, multiplyBy] = vtcBandpassFilter(path, vtc, lowest, highest,TR);
% syntax: [addTo, multiplyBy] = vtcBandpassFilter(path, vtc, lowest, highest,TR)
% default values: TR = 3;
% 
% this func z-normalize each TC for each voxel and then run a bandpass
% filter, (which its limits are mentioned in lowest & highest) on it.
% after this only frequencies between lowest and highest will remain.
% note: 'lowest' & 'highest' are normalized to the niquest frequency
%  the niquest frequency is calculated from the 'TR' variable
% the returned TCs are added an addTo and multiplied by multiplyBy
% in order to exploit the 65,535 limit of the uint16 that is used in the
% format of the VTC files.
% 
% This function was written by:
% Hagar Gelbard
% Rafi Malach's Lab
% Weizmann Institute of Science
% Rehovot, Israel
% [e-mail]

% lowest = 0.015;
% highest = 0.03;

if nargin < 5, TR =3; end

minVal = 0; maxVal = 0;

fclose('all');
[B1, A1] = butter(10,[lowest highest]/(1/TR)*2) ;

% open a new vtc file for writing
fid = fopen([path vtc '.vtc'],'r');

% READ HEADER
[NoOfVolumes, cols, rows, slices, protocol,TR] = vtcMoveToAfterHeaderPosition(fid);

voxelsNo = rows * cols * slices;
startData = ftell(fid);

% for each voxel - read the TC and filter it using a band filter

for voxel = 1:voxelsNo

    fseek(fid, 0, 'cof');
	TC = fread(fid, NoOfVolumes, 'uint16');
    TC = zscore(TC);
	filteredTC = filtfilt(B1, A1, TC) ;
    minVal = min([minVal min(filteredTC)]);
    maxVal = max([maxVal max(abs(filteredTC))]);
    
end

addTo = abs(floor(minVal));
multiplyBy = floor(64000/(maxVal+addTo));

% open new vtc file for writing
nfid = fopen([path vtc '_bf' num2str(lowest) '_' num2str(highest) 'Hz_add' num2str(addTo) '_multiply' num2str(multiplyBy) '.vtc'],'w');
%  WRITE THE HEADER OF THE NEW VTC
vtcWriteHeader(nfid, NoOfVolumes,protocol,TR);

fseek(fid, startData, 'bof');

for voxel = 1:voxelsNo

    fseek(fid, 0, 'cof');
	TC = fread(fid, NoOfVolumes, 'uint16');
    TC = zscore(TC);
	filteredTC = (filtfilt(B1, A1, TC) + addTo)*multiplyBy ;
    fwrite(nfid, filteredTC , 'uint16');
    
end

fclose('all');