%CALLDE BY APPLY_FILTER
%https://ccrma.stanford.edu/~jos/sasp/Example_1_Low_Pass_Filtering.html
% implement a length L=257 FIR low-pass filter having a cut-off frequency at  f_c = 600 Hz. The filter is tested on an input signal  x(n) consisting of a sum of sinusoidal components at frequencies (440, 880, 1000, 2000) Hz. We'll filter a single input frame of length M=256, which allows the FFT to be N=512 samples
function [y]=conv_D(D,ts,tfilt)
echo off
acha=find(~isnan(D));
M = length(acha); 
if M == 0,y=D;return,end

x = D(acha);
Fs = 1/ts;%7 days
% Signal parameters:
%f = [ 440 880 1000 2000 ];      % frequencies   
%M = 256;                        % signal length 
%Fs = 5000;                      % sampling rate

% Filter parameters:
%L=10;%M+1;
L=M;%round((M+1)/2)*2;
%fc=100;
fc=1/tfilt; %365 days
%L = 257;    % filter length 
%fc = 600;   % cutoff frequency 

% Design the filter using the window method:
hsupp = (-(L-1)/2:(L-1)/2);
hideal = (2*fc/Fs)*sinc(2*fc*hsupp/Fs);
h = hamming(L)' .* hideal; % h is our filter

% Choose the next power of 2 greater than L+M-1 
Nfft = 2^(ceil(log2(L+M-1))); % or 2^nextpow2(L+M-1)

% Zero pad the signal and impulse response:
xzp = [ones(1,floor((Nfft-M+1)/2))*x(1) x(:)' ones(1,floor((Nfft-M)/2))*x(end) ]; %added +1 to avoid odd numbers
hzp = [ h zeros(1,Nfft-L) ];

% Transform the signal and the filter:
X = fft(xzp);
H = fft(hzp);

Freq = -fc/2:fc/Nfft:fc/2-1/Nfft;

%Now we perform cyclic convolution in the time domain using pointwise multiplication in the frequency domain:

Y = X .* H;
y = ifft(Y);
relrmserr = norm(imag(y))/norm(y); % check... should be zero
y = real(y);

%y=y(L/2+1:end-(L/2+1));
y=y((Nfft-M+L)/2+1:M-1+(Nfft-M+L)/2+1);
D(acha)=y;%(2:end-1))=y(2:end-1);  %added to avoid last point problems
y=D;

echo on
