https://en.wikipedia.org/wiki/Ringing_artifacts#Introduction
Brickwall filtering in the frequency domain can introduce ringing in
the time domain.
On Mon, Mar 4, 2024 at 5:47 PM Mayank Jog <[log in to unmask]> wrote:
>
> Dear experts,
> Just following up on my query regarding the spm implementation of high pass filtering, in case anyone had insights,
> Thank you,
> Mayank
>
> On Sun, Feb 25, 2024 at 12:49 PM Mayank Jog <[log in to unmask]> wrote:
>>
>> Hello experts!
>> I was trying to understand an oddity I observed with high-pass filtering in spm.
>>
>> Basically, I constructed a signal = y1+ a*y2;
>> y1 = sinusoid whose freq > hpf_cutoff, ie. it shouldn't be filtered out
>> y2 = sinusoid whose freq < hpf_cutoff, ie. it should be filtered out.
>>
>> The issue I'm having is that the filter gives different results based on "a" above (MATLAB code @ end of this email). Thinking from a brick wall** -type filtering POV, this shouldn't happen... the result of filtering "signal" above should be y1 everytime, independent of "a".
>>
>> 1. Reading the documentation, I realized that SPM implements high-pass filtering using DCT.... why do we prefer filtering fMRI data with a DCT filter, since as the above case shows, a brick wall filter seems to be more accurate?
>> 2. Thinking of y2 as "noise", it's almost as if the output is dependent on the scale of noise (captured by the scaling factor "a" above). Is this the right way to think about it/ Am I missing something here?
>>
>> Thank you!
>> Mayank
>>
>>
>> **By brick wall, I mean doing an fft, and nulling all frequencies above hpf_cutoff, followed by an inverse fft.
>>
>> MATLAB Code: ===================
>> L = 1024; %<length of signal
>> filter_100s = ... %< filter with hpf_cutoff = 100s
>> spm_filter( struct('RT', 1, 'HParam', 100, 'row', 1:L) );
>>
>>
>> y1 = sin(2*pi*[1:L]/50 )'; %< sinusoid with period = 50s, shouldn't be filtered
>> y2 = sin(2*pi*[ylim1:L]/350)'; %< sinusoid with period = 350s, should be filtered
>>
>>
>> y_filter_test1 = spm_filter(filter_100s, y1+50*y2); %< 'a' = 50
>> y_filter_test2 = spm_filter(filter_100s, y1+.1*y2); %< 'a'=0.1
>>
>> figure; subplot(3,1,1); plot([y1 y2]);
>> subplot(3,1,2); plot(y_filter_test1);
>> subplot(3,1,3); plot(y_filter_test2);
>> %============================
|