forked from Chaogan-Yan/REST
-
Notifications
You must be signed in to change notification settings - Fork 0
/
rest_IdealFilter.m
executable file
·56 lines (44 loc) · 2.33 KB
/
rest_IdealFilter.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
function [Data_Filtered] = rest_IdealFilter(Data, SamplePeriod, Band)
% FORMAT [Data_Filtered] = y_IdealFilter(Data, SamplePeriod, Band)
% Input:
% Data - 2D data matrix (nDimTimePoints * nTimeSeries)
% SamplePeriod - Sample period, i.e., 1/sample frequency. E.g., TR
% Band - The frequency for filtering, 1*2 Array. Could be:
% [LowCutoff_HighPass HighCutoff_LowPass]: band pass filtering
% [0 HighCutoff_LowPass]: low pass filtering
% [LowCutoff_HighPass 0]: high pass filtering
% Output:
% Data_Filtered - The data after filtering
%-----------------------------------------------------------
% Written by YAN Chao-Gan 120504 based on REST.
% The Nathan Kline Institute for Psychiatric Research, 140 Old Orangeburg Road, Orangeburg, NY 10962, USA
% Child Mind Institute, 445 Park Avenue, New York, NY 10022, USA
% The Phyllis Green and Randolph Cowen Institute for Pediatric Neuroscience, New York University Child Study Center, New York, NY 10016, USA
sampleFreq = 1/SamplePeriod;
sampleLength = size(Data,1);
paddedLength = rest_nextpow2_one35(sampleLength); %2^nextpow2(sampleLength);
LowCutoff_HighPass = Band(1);
HighCutoff_LowPass = Band(2);
% Get the frequency index
if (LowCutoff_HighPass >= sampleFreq/2) % All high stop
idxLowCutoff_HighPass = paddedLength/2 + 1;
else % high pass, such as freq > 0.01 Hz
idxLowCutoff_HighPass = ceil(LowCutoff_HighPass * paddedLength * SamplePeriod + 1);
end
if (HighCutoff_LowPass>=sampleFreq/2)||(HighCutoff_LowPass==0) % All low pass
idxHighCutoff_LowPass = paddedLength/2 + 1;
else % Low pass, such as freq < 0.08 Hz
idxHighCutoff_LowPass = fix(HighCutoff_LowPass * paddedLength * SamplePeriod + 1);
end
FrequencyMask = zeros(paddedLength,1);
FrequencyMask(idxLowCutoff_HighPass:idxHighCutoff_LowPass,1) = 1;
FrequencyMask(paddedLength-idxLowCutoff_HighPass+2:-1:paddedLength-idxHighCutoff_LowPass+2,1) = 1;
FrequencySetZero_Index = find(FrequencyMask==0);
%Remove the mean before zero padding
Data = Data - repmat(mean(Data),size(Data,1),1);
Data = [Data;zeros(paddedLength -sampleLength,size(Data,2))]; %padded with zero
Data = fft(Data);
Data(FrequencySetZero_Index,:) = 0;
Data = ifft(Data);
Data_Filtered = Data(1:sampleLength,:);