%% conv1demo2
%% show operation of convolution on real 
%% signals for small real masks


% define signal, sample rate, sample domain
N = 256;
t = (0:1:(N-1))/N;
sig = 2*cos(5*2*pi*t)+0.5*sin(60*2*pi*t)+ 0.25*sin(100*2*pi*t)+0.5*(2*rand(size(t))-1);
%sig = cos(5*2*pi*t);
%sig = sin(30*2*pi*t);
%sig = sin(60.3*2*pi*t);
%k = 96; sig = [ones(1,k),zeros(1,N-k)]*(N/k);
%sig = rand(1,N);
 %sig = 2*rand(1,N)-1;


% define convolving mask
% mask1 is the causal part (past and present) 
% and mask2 is the anti-causal (future) part
eps = 10^(-5);
 mask1 = [2,1]/3; mask2 = [1]/3;
% mask1 = [2,-1];mask2 = [-1];
%mask1 = [1,-3,3,-1]; mask2 = [ ];

%  mask1 = [2,2,1]/8, mask2 = [1,2]/8
%  mask1 = [1,1]/2; mask2 = [];

mask = [mask1,zeros(1,N-length(mask1)-length(mask2)),mask2];
mask = mask/sum(abs(mask(:)));
posmaskinds = find(mask > eps)-1;
negmaskinds = find(mask < -eps)-1;
maskinds = [(-length(mask2)):(-1),0:(length(mask1)-1)];
 
% compute convolution and dft
fsig =circconv(sig,mask);
sfmask = fftshift(fft(mask));
sfsig = fftshift(fft(sig));



% iteration and graphing parameters
waitforkey = [0,1,2,3,4,5]
dt =0.1;
M = max(abs(sig(:)));
Mf = max(abs(fsig(:)));
SD = 0:(N-1);
FD = SD-(N/2);
axis tight

% display results 
subplot(3,1,2)
plot(FD, abs(sfmask),'r-')
title('amplitude of dft of mask');
axis tight
subplot(3,1,3)
plot(FD, abs(sfsig),'r-')
title('amplitude of dft of signal');
axis tight




for k = 0:(N-1)
    krange = 1:(k+1);
    pmik = mod((k-posmaskinds),N)+1;
    nmik = mod((k-negmaskinds),N)+1;
    subplot(3,1,1)
    cla
    hold on
    plot(SD,sig,'b-')
    plot(SD(pmik), sig(pmik),'k.')
    plot(SD(nmik), sig(nmik),'r.')
    plot(SD(krange), fsig(krange),'k-')
    plot(SD(k+1), fsig(k+1),'g.')
    axis([0,N-1,-M,M]); 
    title('original signal, black > 0, red <0');
    hold off

    pauseforkey = sum(find(waitforkey==k)); 
    if pauseforkey > 0
      disp('Press a key to proceed'), pause,
      else pause(dt)
      end
   end