clear all;
close all;

s=load('090908-091024.profile');
t=load('DEM.profile');
model = load('VISCO.profile');
GPS = load('profile_GPS.txt');

%-- remove NaNs
x=s(:,1);
y=s(:,2);
x = x(~isnan(y));
y = y(~isnan(y));

%-- Interpolate SAR signal
s1=interp1(x,y,t(:,1));
t1=t(:,2);

%-- remove NaNs once again
x=s(:,1);
y1=s1;
y2=t1;
x = x(~isnan(y1));
y2 = y2(~isnan(y1));
y1 = y1(~isnan(y1));

%-- standard deviation of DEM
s2=std(y2);

%-- new coherence
%COR = real(fft(y1).*conj(fft(y2)))./abs(fft(y1))./abs(fft(y2));
COR12 = (fft(y1).*conj(fft(y2)));
COR11 = (fft(y1).*conj(fft(y1)));
COR22 = (fft(y2).*conj(fft(y2)));

%-- cut off freq
dx = 2;
dX = 1;
%L = length(COR11);
%COR(1:dx) = 0; COR(L-dx:L) = 0;
%COR(int16(L/2)-dx:int16(L/2)+dX) = 0;
%COR = smooth(COR);

    myFilt = ones(dx) ./ dx;
    COR12 = imfilter(COR12, myFilt);
    COR11 = imfilter(COR11, myFilt);
    COR22 = imfilter(COR22, myFilt);
%    COR12 = smooth(COR12);
%    COR11 = smooth(COR11);
%    COR22 = smooth(COR22);
    RCOR = real(COR12)./sqrt(COR11.*COR22);
    ICOR = imag(COR12)./sqrt(COR11.*COR22);
    COR = RCOR.^2+ICOR.^2;
    COH = ifft((COR12)./sqrt(COR11.*COR22),'symmetric');

filt = 1 - COR;

fy1=real(ifft(fft(y1).*filt));
fy2=real(ifft(fft(y2).*filt));

figure;
subplot(2,1,1);
plot(x,y1,'b',x,fy1,'r',model(:,1),model(:,2),'g',GPS(:,4),GPS(:,3)/10,'ko','MarkerSize',10);
set(gca,'YLim',[-10 10]);
legend('SAR initial','SAR final','VISCO1D','GPS');
subplot(2,1,2);plot(x,y2,'b',x,fy2,'r');
legend('DEM initial','filtered DEM');

figure;
subplot(2,1,1);
plot(x,y1-fy1); legend('SAR difference');
subplot(2,1,2); plot(x,y2-fy2); legend('DEM difference');

figure; plot(y1,y2,'.');hold on; plot(fy1,fy2,'r.');

figure; subplot(2,1,1);plot(RCOR,'k');hold on;plot(ICOR,'r');hold off
subplot(2,1,2); plot(real(COH),'k');hold on;plot(imag(COH),'r');hold off
