// Testing the function timevarconv()
// 
// -----------------------
// Federico Miyara
// 2025-06-04
// 2025-06-17
// 2025-06-27

// Clean start
clear
close(winsid())


// ------------------------------------------
// Create signal
// ------------------------------------------



// Signal
caso = 1
select caso
case 1
   // Independent variable and its length
   K = 0:1000;
   N = length(K)
   x = 1 + 0.3*sin(2*%pi*3*K/N) + 0.3*sin(2*%pi*30*K/N);
case 2
   // Independent variable and its length
   K = 0:1000;
   N = length(K)
   x = 1 + 0.3*sin(2*%pi*40*K/N);
end


// ------------------------------------------
// Time-variable impulse response
// ------------------------------------------

// Empty cell-array to contain an impulse
// response for each value of the independent 
// variable 
h = cell(1,N);

caso = 2
select caso
case 1
   // Linear decay impulse responses
   tic
   for k = 1:N
      // Length of the k-th impulse response
      Jo(1,k) = ceil(k/5);
      // Create the k-th impulse response as
      // a linear decay to 0 with a slope that
      // is inversely proportional to k. The 
      // sample at which 0 value is reached 
      // is excluded.   
      h{k} = 1 - (0:Jo(k)-1)*5/k;
      // Compute the k-th normalization constant
      hnorm = sum(h{k});
      // Normalize impulse response
      h{k} =  h{k} / hnorm; 
   end
   toc
case 2
   // Exponential decay impulse responses
   // Initial and final attenuation constant
   mo1 = 2;
   mo2 = 60;
   tic
   for k = 1:N
      // Exponential attenuation constant 
      // (increases linearly from mo1 to mo2)
      mo = mo1 + (mo2 - mo1)/(N - 1) * (k - 1);
      // Impulse response length (ensures 
      // vanishing to 0.67 %)
      Jo(1,k) = 5*mo;
      // Exponential impulse response
      h{k} = exp(-(0:Jo(k)-1)/mo); 
      // Compute the k-th normalization constant
      hnorm = sum(h{k});
      // Normalize impulse response
      h{k} =  h{k} / hnorm; 
   end
   toc
case 3
   // Constant impulse response
   for k = 1:N
      h{k} = exp(-(0:20)/4);
      // Compute the k-th normalization constant
      hnorm = sum(h{k});
      // Normalize impulse response
      h{k} =  h{k} / hnorm; 
   end
end



// ------------------------------------------
// Time-variable convolution 
// ------------------------------------------

// Load function
exec timevarconv.sci;

// Compute the convolution
y = timevarconv(x, h);

// In the case of constant impulse response
// it can be computed using the standard 
// convolution function
yy = convol(x, h{50});

// ------------------------------------------
// Plot the signal before and after the
// time-variable convolution 
// ------------------------------------------

scf(1);
clf(1);
gcf().figure_position = [200,100];
subplot(2,1,1)
plot(K, x)
title("Input signal")
xlabel("k")
ylabel("x")  
gca().grid = [1 1];  
gca().data_bounds = [0 0; 1000 2]; 
subplot(2,1,2)
plot(0:length(y)-1, y, 0:length(yy)-1, yy, "r")  
title("Output signal after convolution (red: constant; blue: variable)")
xlabel("k")
ylabel("y = x*h")
gca().grid = [1 1];  
gca().data_bounds = [0 0; 1000 2]; 

// --------------------------------------
// Test delay and anticipation
// --------------------------------------
K = 1:100;
x = zeros(K);
Lx = length(x);
x([1 25 75]) = 1;  
h = cell(1,Lx);
Jo = zeros(1,Lx);
d1 = zeros(1,Lx);
d2 = zeros(1,Lx);

for m=1:length(x)
   h{m} = ones(1, 4 + floor(m/2));
   Jo(m) = length(h{m});
   d1(m) = 0;   
   d2(m) = -round(Jo(m)/2);
end

exec timevarconv.sci;
A1 = max(1 - min((1:Lx) + d1), 0);
A2 = max(1 - min((1:Lx) + d2), 0);
y1 = timevarconv(x, h, d1); 
y2 = timevarconv(x, h, d2); 

scf(2);
clf(2);
gcf().figure_position = [220,80];
subplot(2,1,1)
plot(K, x)
gca().data_bounds = [0,0;Lx,1.1];
subplot(2,1,2)
plot(K, y1(K), K, y2(K), "r")
gca().data_bounds = [0,0;Lx,1.1];

