% First an easy case, to make sure we have normalization and other details correct.

dtheta = 2*pi/20; 
theta = 0:dtheta:20*pi;
signal = cos(theta);
% The above signal should have zero mean over one or several periods:
Nper = 1
N = round(Nper * 2*pi / dtheta);    % filter size = Nper periods.
[integral, kernel] = ifilter(signal*dtheta, N);
   % Could equivalently do ifilter(signal,N) * dtheta
integral2 = cumsum(signal) * dtheta;  % easy way.
integral2 = integral2 - mean(integral2);
answer = sin(theta);  % exact answer, for comparison.

fig1 = figure(1)
hold off
plot(theta, signal, 'g-')
hold on
xlabel('\theta')
plot(theta, answer, 'b-')
plot(theta, integral, 'b+')
plot(theta, integral2, 'r')
plot(theta, kernel(1:numel(theta)), 'm')
title('verification of method', 'fontsize', 20)
legend('signal','exact','ifilter','cumsum','kernel')
print('../plots/ifilter_test.pdf', '-dpdf')

% What happens if we add noise?
rms_noise = 0.5;
% noise = normrnd(0.0,rms_noise,size(signal));
noise = rms_noise * randn( size(signal) );
signal = signal + noise;

[integral, kernel] = ifilter(signal*dtheta, N);
   % Couild equivalently do ifilter(signal,N) * dtheta
integral2 = cumsum(signal) * dtheta;  % easy way.

figure(2)
hold off
plot(theta, signal, 'g-')
hold on
xlabel('\theta')
plot(theta, answer, 'b-')
plot(theta, integral, 'b+')
plot(theta, integral2, 'r')
plot(theta, kernel(1:numel(theta)), 'm')
title('test with noise', 'fontsize', 20)
legend('signal','exact','ifilter','cumsum','kernel')
print('../plots/ifilter_test_noise.pdf', '-dpdf')
