clear all
close all

load sndb_prbs;
u=u.';
y=y.';

figure;
subplot(211); plot(u); title('u(n)'); grid;
subplot(212); plot(y); title('y(n)'); grid;

M = 40;          % dlugosc odpowedzi impulsowej
L = length(u);   % d�ugo�� sygna��w
L = 2*M;         % liczba r�wna�, minimum M, maksimum L-M+1

% Metoda parametryczna model FIR o dlugosci M
%U=[];
%for i=0:L-1
%    U=[U; u(M+i:-1:i+1)];
%end
U = toeplitz(u(M:M+L-1),u(M:-1:1));

Y = y(M:M+(L-1))';
b = U\Y;            %  pinv(X)*y;  % inv(X'*X) * X'*y = Rxx * rxy

fpr = 44100; dt=1/fpr; t=dt*(0:M-1);

% Odpowied� impulsowa
figure; stem(t,b); title('h(t)'); xlabel('t[s]'); grid;

% Jej widmo
Nfft=1000; Nf=Nfft/2+1;
df=fpr/Nfft; f1=0:df:fpr-df; f3=f1;
% FFT
H1=fft(b,Nfft);
% FREQZ
[ H2, f2 ] = freqz(b(1:M),1,Nf,fpr);

% Metoda nieparametryczna dla porownania
M = L;                 % od L do Nfft/2
ryu=xcorr(y,u,M);      % d�ugie y i u! przesuni�te wektory, je�li nie trafisz w pobli�e to nie zrobisz z kr�tkich u(M:L)
ruu=xcorr(u,u,Nfft/2);
H3 = fft( ryu, Nfft ) ./ fft( ruu, Nfft );

figure
plot(f1,20*log10(abs(H1)),'ro-',f2,20*log10(abs(H2)),'kx-',f3,20*log10(abs(H3)),'b-')
title('|H(f)| [dB]'); xlabel('f[Hz]'); grid;

% Metoda dla coraz wi�kszej liczby r�wna�
LL = [50, 100, 200, 1000 ];
y = y + 0.1*randn(1,length(y));
for L=LL
    U=[];
    for i=0:L-1
        U=[U; u(M+i:-1:i+1)];
    end
    Y=y(M:M+(L-1))';
    b=U\Y;
    [ H2, f2 ] = freqz(b(1:M),1,Nf,fpr);
    figure
    subplot(211)
    stem(1:M,b)
    subplot(212)
    plot(f2,20*log10(abs(H2)),'r-')
end




/////////////// zad 3 /////////////////



clear all; close all;
i1=imread('Krata1.jpg');
figure; imshow(i1),
i1g=double(rgb2gray(i1));
[sx, sy]=size(i1g);
cx=sx/2;
cy=sy/2;
% wyszukanie przeci�� linii
wzorzec=[
0 0 0 1 1 1 0 0 0
0 0 0 1 1 1 0 0 0
0 0 0 1 1 1 0 0 0
1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1
0 0 0 1 1 1 0 0 0
0 0 0 1 1 1 0 0 0
0 0 0 1 1 1 0 0 0
];
ic=normxcorr2(wzorzec,i1g);
ig=ic<-0.7;
ig=bwmorph(ig,'shrink',Inf);
[ix,iy]=find(ig>0);
r=sqrt((ix-cx).^2+(iy-cy).^2);
[rs,i]=sort(r);
% ustalenie punkt�w odniesienia (siatki regularnej)
x0=ix(i(1)); y0=iy(i(1));
dx=abs(ix(i(1))-ix(i(4)))-1;
dy=abs(iy(i(1))-iy(i(4)))-3;
jx=x0+dx*[-4:1:5];
jy=y0+dy*[-6:1:7];
[jX, jY]=meshgrid(jx,jy);
jx=jX'; jx=jx(:);
jy=jY'; jy=jy(:);
rm=sqrt((jx-cx).^2+(jy-cy).^2);
[rms,j]=sort(rm);
figure; imshow(ig), hold on
plot(cy,cx,'r+') % punkt centralny
plot(jY,jX,'bo') % punkty odniesienia
% identyfikacja modelu algorytmem LS
% uzupe�nij wyliczanie modelu stopnia 3 zale�no�ci rs od rms
% warto�ci wyliczone z modelu zapisz do zmiennej R
%???
%???
%???
figure(1)
plot(rms,rs,'b.', rms,R,'g-', [0 3e2],[0 3e2],'r-')
% prostowanie obrazka
xx=1:sx;
yy=1:sy;
[XX,YY]=meshgrid(xx,yy);
RR=sqrt((XX-cx).^2+(YY-cy).^2);
RRr=RR*A(1)+RR.^2*A(2)+RR.^3*A(3);
RRu=RRr./RR;
i1r=interp2(i1g,(YY-cy).*RRu+cy,(XX-cx).*RRu+cx);
figure,
    subplot(2,1,1),imagesc(i1g)
    subplot(2,1,2),imagesc(i1r)
    colormap gray