% correct_signal_example.m

%

% This MATLAB script file demonstrates the

% use of Fourier analysis to correct a

% measured signal for corruption from

% the measuring device.  This example

% uses only deconvolution to remove

% the measurement error.

%

% K. Beers. MIT ChE. 8/15/2002

% First, we generate the input signal.

N = 2^15;

P = 2*pi;

t_val = linspace(0,4*pi,N);

dt = t_val(2) - t_val(1);

t0 = t_val(1);

f = zeros(size(t_val));

for k=1:10

    f = f + (1/sqrt(k))*sin(k*t_val);

end

noise_mag = 0.2;

f = f + noise_mag*randn(size(f));

% From the details of the signal sampling,

% we create a vector of the sampled frequencies

% in the Fourier spectrum.

omega_crit = pi/dt;

d_omega = pi/P;

omega = linspace(0,2*omega_crit-d_omega,N);

% We now generate the Fourier transform

% of the response function.

r_FT = zeros(size(omega));

omega_dev = 5;  % frequency of max. device sensititivy

std_dev = 2.5*omega_dev;  % standard dev. of sensitivity window

r_cut = sqrt(eps);  % cutoff for response magnitude

for m=1:(N/2)

    omega_m = omega(m);

    var1 = exp(-(omega_m-omega_dev)^2/(2*std_dev));

    if(var1 < r_cut)

        var1 = 0;

    end

    % set positive frequency

    r_FT(m) = var1;

    % set negative frequency

    r_FT(N-m) = r_FT(m);

end

r_FT(N) = r_FT(1);

% calculate power spectrum

r_PS = real(conj(r_FT).*r_FT);

% calculate time-domain response function

r = real(ifft(r_FT));

i_r_zero = 1;

% We now plot this synthesized response function

figure;

% frequency domain

subplot(2,2,1);

plot(omega,r_PS);

buffer = 0.05*(max(r_PS)-min(r_PS));

axis([omega(1) omega(N) ...

        (min(r_PS)-buffer) (max(r_PS)+buffer)]);

xlabel('\omega');

ylabel('|R(\omega)|^2');

title('Device response function');

% time domain

subplot(2,2,2);

plot(t_val,r);

buffer = 0.05*(max(r)-min(r));

axis([t_val(1) t_val(N) ...

        (min(r)-buffer) (max(r)+buffer)]);

xlabel('t');

ylabel('r(t)');

% postive low-freq response

subplot(2,2,3);

plot(omega,r_PS);

buffer = 0.05*(max(r_PS)-min(r_PS));

axis([omega(1) 3*omega_dev ...

        (min(r_PS)-buffer) (max(r_PS)+buffer)]);

xlabel('\omega');

ylabel('|R(\omega)|^2');

% negative low-freq response

subplot(2,2,4);

plot(omega,r_PS);

buffer = 0.05*(max(r_PS)-min(r_PS));

axis([(omega(N)-3*omega_dev) omega(N) ...

        (min(r_PS)-buffer) (max(r_PS)+buffer)]);

xlabel('\omega');

ylabel('|R(\omega)|^2');

% We now use convolution to generate the measured

% output signal of the device.

f_FT = fft(f);

d_FT = r_FT.*f_FT;

d = real(ifft(d_FT));

% power spectra of input and output signals

f_PS = real(conj(f_FT).*f_FT);

d_PS = real(conj(d_FT).*d_FT);

% Next, plot synthetic input and output signals in

% Fourier and time domains.

figure;

% input - real time

subplot(2,2,1);

plot(t_val,f);

xlabel('t');

ylabel('f(t)');

title('Measurement device input and output signals');

% output - real time

subplot(2,2,2);

plot(t_val,d);

xlabel('t');

ylabel('d(t)');

% input - frequency

subplot(2,2,3);

plot(omega,f_PS);

buffer = 0.05*(max(f_PS)-min(f_PS));

axis([omega(1) 3*omega_dev ...

        (min(f_PS)-buffer) (max(f_PS)+buffer)]);

xlabel('\omega');

ylabel('|F(\omega)|^2');

% output - frequency

subplot(2,2,4);

plot(omega,d_PS);

buffer = 0.05*(max(d_PS)-min(d_PS));

axis([omega(1) 3*omega_dev ...

        (min(d_PS)-buffer) (max(d_PS)+buffer)]);

xlabel('\omega');

ylabel('|D(\omega)|^2');

% In reality, the measured device output signal would itself

% contain some random noise, but we assume that this could be

% removed by filtered.  Therefore, we will take the calculated

% device output d(t), and use it directly, with the response

% function to reconstruct the signal.

% Also, we usually do not have a priori knowledge of the

% device response function, but we can obtain it through

% calibration, by inputing a known signal, measuring the

% output, and then devolve the two signals to estimate

% the response function.

% We assume now that we only have information about the

% response function and the measured output signal.

% Using this information, we now wish to reconstruct the

% actual signal input to the system.

% We have already calculated the discrete Fourier transforms

% of the response function and the measured output signal.

% From these, we use deconvolution to obtain the estimated

% Fourier transform of the input signal.

f_est_FT = zeros(size(d_FT));

for m=1:length(d_FT)

    if(abs(r_FT(m)) >= sqrt(eps))

        f_est_FT(m) = d_FT(m)/r_FT(m);

    end

end

f_est_PS = real(conj(f_est_FT).*f_est_FT);

% We now obtain the corrected signal by an inverse

% Fourier transform.

f_est = real(ifft(f_est_FT));

% We now plot the corrected signals, and compare the

% Fourier transform to those of the "true" signal

% and the measured device output.

figure;

% corrected signal - time domain

subplot(2,2,1);

plot(t_val,f_est);

xlabel('t');

ylabel('f(t) - corrected estimate');

title('Removal of device measurement bias');

% corrected signal - frequency domain

subplot(2,2,2);

plot(omega,f_est_PS);

buffer = 0.05*(max(f_est_PS)-min(f_est_PS));

axis([omega(1) 3*omega_dev ...

        (min(f_est_PS)-buffer) (max(f_est_PS)+buffer)]);

xlabel('\omega');

ylabel('|F(\omega)|^2 - estimated');

% measured device output - frequency domain

subplot(2,2,3);

plot(omega,d_PS);

buffer = 0.05*(max(d_PS)-min(d_PS));

axis([omega(1) 3*omega_dev ...

        (min(d_PS)-buffer) (max(d_PS)+buffer)]);

xlabel('\omega');

ylabel('|D(\omega)|^2');

% original signal - frequency domain

subplot(2,2,4);

plot(omega,f_PS);

buffer = 0.05*(max(f_PS)-min(f_PS));

axis([omega(1) 3*omega_dev ...

        (min(f_PS)-buffer) (max(f_PS)+buffer)]);

xlabel('\omega');

ylabel('|F(\omega)|^2');