CODEWITHSUNDEEP
matlabfilter
Graham
Graham

Reputation: 581

Implementing a Sum of Sincs filter in Matlab

I am trying to implement a filter that looks like

SoS Filter

in Matlab. I have:

omegas = (-length(t)/2:length(t)/2)*2*pi/tau/10;
SOS_freq_sum = zeros(1,length(omegas));
for i = 1:length(K_set)
   k = K_set(i);
   SOS_freq_sum = SOS_freq_sum + b_k(i)*sinc(omegas/(2*pi/tau)-k);
end
SOS_filter_in_frequency = (tau/sqrt(2*pi))*((SOS_freq_sum)');

How do I use this to filter my data? Matlab's filter function defines numerator and denom coefficients for a transfer function, but not for the SOS form. Is there any way to do it without transforming the input signal into the frequency domain?

Many thanks

Upvotes: 0

Views: 190

Answers (1)

nivag
nivag

Reputation: 573

Filtering a signal can be viewed as convolving your time domain signal with with a (time domain) filter kernel. If you had the filter kernel of your signal you could perform the convolution with conv().

Similarly convolution in the time domain is multiplication in the frequency domain so you can also multiply the frequency response of your filter with the fourier transform of the signal. And then inverse fourier transform the result to get your filtered signal.

You appear to have calulated the frequency response of your filter. So you can either fourier transform your signal, multiply (with suitable zero adding if necessary) and then inverse fourier transform, or inverse fourier transform your filter to get the kernel and convolve it with the signal.

For large datasets and filters performing the convolution in the frequency domain can be significantly faster but I doubt this will be noticeable for 1D signals.

Outline code:

ftsignal=fft(signal);
ftsignal=fftshift(signal); %I think this makes it the same as your frequency response.
ftfiltsignal=ftsignal.*frequecyresponse; %if these are not the same size you will need to zero-pad
filtsignal=ifft(ftfiltsignal);

If you really want to avoid convolution in the frequency domain for some reason:

filterkernel=ifft(frequencyresponse);
filteredsignal=conv(signal, filterkernel, 'same');

Note this approach will produce edge effects which can be significant if you filter size is reasonaby large compared to your signal length.

Upvotes: 1

Related Questions