image thumbnail

Filter for discrete time signal

by

 

20 Jun 2013 (Updated )

input: signal, Fs, filter: specific frequency, # of harmonics, shape, width output : a clean signal

[y_fl, Y_fl, f, Y, f_fl, fl_function, filter]=f4ts(y,Fs,filter,plot1)
function [y_fl, Y_fl, f, Y, f_fl, fl_function, filter]=f4ts(y,Fs,filter,plot1)
%Itamar Luzon
%update 23.06.13  @HUJI
%Filter for discrete time signal -
%   give you a clean signal, by a filter.
%   the filter width is at FWHM!!!, and the high is 1.
%   [y_fl, Y_fl, f, Y, f_fl, fl_function, filter]=f4ts(y,Fs,filter,plot1)
%   note : 
%   1. this function does not check the input signal, so check it before!!!
%   2. this function halnde better with noise at diffrent times, than
%   without noise.
%   3. when you try it, choose diffrent shapes, sizes of width for the filter,
%   and number of harmonics, to get the best one filter.
%   4. don't forget to clear the filter before you use it agine.
%
%   input :
%   y : the signal in time domain
%   Fs : frequency of sempeling
%   the filter is the type and the shape of the wanted filter at the wanted
%   frequency, (the child field must be at the same names)
%   filter.h : is the first harmonic frequency that you want to filter(real posative).
%   filter.nh : is the number of harmonics that you want to filter (natural number).
%   if "filter.nh" does not exist it's use default of ceil(max(f)/filter.h).
%   filter.shape : is the shape of the filter that you want to use,
%   the avilable shapes are : 'gaussian', 'lorentzian', 'square'.
%   if there is no "filter.shape" field the deafult will be 'square' shape
%   filter.width : is the shapes bandpass width, using Full Width Half Max(real posative),
%   if "filter.width" does not exist it's use default of filter.h/100.
%   plot1 : is for testing, if plot1==1, it's plot a fig with result
%
%   output :
%   y_fl : is a clean time domain signal, after the filter.
%   Y_fl : is the clean frequency domain signal, after the filter.
%   f : is the x-axis of frequency domain, include negative part!!!!
%   Y : is the frequency domain signal (no filter!!!)
%   f_fl : is the filter domain that was used, note that it made from the
%   X-axis (f) , mean that beacouse of the "df" will not sample ALL the
%   dots. the real filter functin shape is at fl_function
%   fl_function : is the filter function
%
% test for :
% signal at frequencys of 15Hz with 8 harmonics, with lifetime of 2sec
% with noise in signal of 10%,also in every harmonic phase of 10%
% also noise of 4 goussians in diffrent times and amplitude
% also noise at the frequency of 13.5,21Hz
% dt=0.001; t=0:dt:10; noise=0.1; l=length(t);
% y=(3+...
%     -1*cos(2*pi*13.5*t+2*pi*(1+randn(1,l)*noise)).*exp(-t/0.1)+(...
%     +13*sin(2*pi*15.*t+2*pi*(1+randn(1,l)*noise))...
%     +16*sin(2*pi*30.*t+2*pi*(1+randn(1,l)*noise))....
%     +10*sin(2*pi*45.*t+2*pi*(1+randn(1,l)*noise))....
%     +8*sin(2*pi*60.*t+2*pi*(1+randn(1,l)*noise))....
%     +5*sin(2*pi*75.*t+2*pi*(1+randn(1,l)*noise))....
%     +3*sin(2*pi*90.*t+2*pi*(1+randn(1,l)*noise))....
%     +2*sin(2*pi*105.*t+2*pi*(1+randn(1,l)*noise))....
%     +0.5*sin(2*pi*120.*t+2*pi*(1+randn(1,l)*noise))....
%     ).*exp(-t/2)+...
%     3*sin(2*pi*21.*t+2*pi*(1+randn(1,l)*noise))+...
%     -12*exp(-((t-0.5)/0.5).^2)... %gaussian at 0.5 width 0.5
%     +10*exp(-((t-1)/1).^2)... %gaussian at 1 width 1
%     -10*exp(-((t-3)/2).^2)... %gaussian at 3 width 2
%     +15*exp(-((t-6)/0.2).^2)... %gaussian at 6 width 0.2
%     -20*exp(-((t-8)/0.1).^2)... %gaussian at 8 width 0.1
%     -15*exp(-((t-9)/3).^2)... %gaussian at 9 width 3
%     ).*(1+2*randn(1,l)*noise);
% 
% Fs=1/dt;
% filter.h=15;
% filter.width=3;
% filter.nh=8;
% filter.shape='lorentzian';
% f4ts(y,Fs,filter,1);
Y=fft(y);
N=length(y); dt=1/Fs; df=1/dt/N;
if mod(N,2)==0
    f=[0:floor(N/2) -floor(N/2-1):-1]*df;
else
    f=[0:floor(N/2) -floor(N/2):-1]*df;
end
[filter, goodFilter]=fixfilter(filter,f);
if goodFilter
    harmonics=[-filter.nh:1:-1 , 1:1:filter.nh]*filter.h;
    fl_function_temp=[];
    for harmonic=harmonics
        switch filter.shape
            case 'gaussian'
                temp=['+exp(-((4*log(2)*(x-',num2str(harmonic),').^2)/((',num2str(filter.width),')^2)))'];
            case 'lorentzian'
                temp=['+1./(1+4*((x-',num2str(harmonic),')/',num2str(filter.width),').^2)'];
            case 'square'
                temp=['+abs(heaviside(x+' num2str((-harmonic+filter.width/2))...
                    ')-heaviside(x+' num2str((-harmonic-filter.width/2)),'))'];
        end
        fl_function_temp=[fl_function_temp temp];
    end;
    fl_function_temp=['fl_function = @(x) ' fl_function_temp ';'];
    eval(fl_function_temp);
    ll=length(f);
    f_fl=zeros(1,ll);
    II=2:ceil(ll/2);
    f_fl(II)=fl_function(f(II));
    f_fl(ll-II+2)=f_fl(II);
    Y_fl=Y.*f_fl;
    if (nargout>1 || (exist('plot1','var') && plot1))
        y_fl=ifft(Y_fl);
    end
    if exist('plot1','var') && plot1 showIt(f,f_fl,Y_fl,Y,filter,y_fl); end;
end
end

%check if the filter is good and add defaults if needed
function [filter, goodFilter]=fixfilter(filter,f)
inputTestFilterFieldsNames={'h','nh','shape','width'};
inputTestFFiledsExist = isfield(filter,inputTestFilterFieldsNames);%for h: check if ok, for shape and width if not ok is set to default
inputTestdefault=zeros(1,4); inputTestFFiledsCheck=zeros(1,4);%check if ok
if inputTestFFiledsExist(1)
    inputTestFFiledsCheck(1) = isa(getfield(filter,'h'),'double') && imag(filter.h)==0 && filter.h>0;
else
    inputTestdefault(1)=1;
end
if inputTestFFiledsExist(2)
    inputTestFFiledsCheck(2) = isa(getfield(filter,'nh'),'double') && filter.nh>0 && mod(filter.nh,1)==0;
elseif ~inputTestFFiledsExist(2) && inputTestFFiledsCheck(1)
    filter.nh=ceil(max(f)/filter.h); inputTestFFiledsCheck(2)=1; inputTestdefault(2)=1;
end

if inputTestFFiledsExist(3)
    inputTestFFiledsCheck(3) =isa(getfield(filter,'shape'),'char') && (...
        (~isempty(strfind(getfield(filter,'shape'), 'lorentzian'))) ||...
        (~isempty(strfind(getfield(filter,'shape'), 'square'))) ||...
        (~isempty(strfind(getfield(filter,'shape'), 'gaussian'))) );
else
    filter.shape='square'; inputTestFFiledsCheck(3)=1; inputTestdefault(3)=1;
end

if inputTestFFiledsExist(4)
    inputTestFFiledsCheck(4) = isa(getfield(filter,'width'),'double') && filter.width>0;
else
    %width default
    filter.width=filter.h/100; inputTestFFiledsCheck(4)=1; inputTestdefault(4)=1;
end
goodFilter=(inputTestFFiledsExist(1)==1 && sum(inputTestFFiledsCheck)==4);


for ii1=1:3
    if inputTestdefault(ii1)
        if 1<ii1
            disp(['the field "filter.' inputTestFilterFieldsNames{ii1} '" does not exist, so default.']);
        else
            disp(['the field "filter.' inputTestFilterFieldsNames{ii1} '" does not exist.']);
        end
    end
    if ~inputTestFFiledsCheck(ii1)
        disp(['problem with the field "filter.' inputTestFilterFieldsNames{ii1} '" value.']);
        if ii1==3 disp('can be : ''gaussian'', ''lorentzian'', ''square''.'); end
    end
end
end

function showIt(f,f_fl,Y_fl,Y,filter,y_fl)
figure;
subplot(1,2,1); plot(f,f_fl*max(abs(Y_fl)),'y',f,abs(Y),f,abs(Y_fl)); xlabel('f [Hz]');
legend(['filter ampX' num2str(max(abs(Y_fl)))],'f','filtered f');
title(['filter width@FWHM ' num2str(filter.width)]);
dt=1/(length(f)*(f(2)-f(1)));
t=(0:(length(f)-1))*dt;
ytemp=ifft(Y);
subplot(1,2,2); plot(t,ytemp,t,y_fl); xlabel('time [sec]');
xlim([t(1) t(end)]);
legend('time','filtered time');
hold off
end

Contact us