SHARE
TWEET

Untitled

a guest Jan 22nd, 2020 74 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1.  
  2. clear;
  3. t=0:0.01:1;% wektor z wartosciami od 0 do 1 rosnace o 0.01
  4. X=0.5*sin(2*pi*15*t)+2*sin(2*pi*t)+sin(2*pi*7*t);  %tworzymy sobie sygnal wyjsciowy - suma trzech sinusoid
  5.  
  6. windowSize = 10;                  
  7. Num = (1/windowSize)*ones(1,windowSize);     %Filtr średniej ruchomej przesuwa okno długości WindowSize wzdłuż danych, obliczając średnie dane zawarte w każdym oknie.
  8. windowSize1 = 5;
  9. Num1 = (1/windowSize1)*ones(1,windowSize1);
  10. load FFR;
  11. Pp1=filter(Num,1,X) ;           %filtr dolnoprzepustowy Fc:10 Hz
  12. load FFR;
  13. Pp2=filter(Num1,1,X);            %filtr dolnoprzepustowy Fc: 5 Hz
  14.  
  15.  
  16. Lmax=15;                        %poczatkowo przyjmowany maksymalny rzad macierzy P
  17. N=size(t,2);                    %N - liczba wykorzystanych probek sygnalu
  18.  
  19.  
  20. Ptp=returnPtP(N,Lmax,Pp1,Pp2);    %Ta funkcja zwraca macierz, ktora powstaje z pomnozenia P transpon. z P (czyli PtP)
  21.  
  22.  
  23. [e_vec,e]=eig(Ptp);               %tu szukamy wartosci wlasnych macierzy PtP
  24.                                 %(zmienna e) i jej wektorów wlasnych (zmienna e_vec)
  25.                                 %smienne e i e_vec są macierzami
  26.  
  27. eMax = max(max(e));             %szukamy najwiekszej wartosci wlasnej macierzy PtP
  28.                                 %pierwsze max(e) daje nam wektor z
  29.                                 %najwiekszymi wartościami w kolumnach
  30.                                 %macierzy e, drugie max(max(e)) daje nam
  31.                                 %wartość maxymalną z wektora max(e)
  32.                                                                                        
  33.  
  34. %%%%%%%ponizsza petla wyznacza ile jest wartosci wlasnych macierzy PtP,
  35. %%%%%%%ktore sa wieksze lub równe 5% najwiekszej wartosci wlasnej
  36. k=0;
  37. for i=1:size(e)
  38.    
  39.     if e(i,i)>=0.05*eMax
  40.         k=k+1;
  41.     end
  42.  
  43. end
  44.  
  45. odd_or_not = mod(k,2);
  46.  
  47. %%%%%%%%%%%%%%%Okreslamy nowt rzad macierzy
  48.  
  49. if odd_or_not==1
  50. k=k+1;          %jesli petla zliczy nieparzysta liczbe to zaokraglamy w gore
  51. L=0.5*k;
  52. else
  53. L=0.5*k;        %polowa wartosci wlasnych zliczonych przez powyzsza petle to nowy rzad filtru
  54. end
  55.  
  56. %%%%%%%%%%%%%%%Ponownie liczymy macierz PtP juz majac wlasciwy oszacowany rzad (moze on wyjsc duzo nizszy niz rzeczywisty rzad filtrow)
  57.  
  58. Ph=returnPtP(N,L,Pp1,Pp2);
  59.  
  60. [V,E]=eig(Ph);  % zwraca diagonalna macierz wartosci wlasnych E i macierz V gdzie kolumny macierzy V sa odpowiednimi wektorami wlasnymi, macierzy PtP
  61.  
  62. E_dgl = diag(E);                            %otrzymujemy wektor zawierajacy wartosci z przekatnej macierzy diagonalnej E
  63. Emin=min(E_dgl);                            %szukamy najmniejszej wartosci wlasnej macierzy PtP
  64. [Emin_row,Emin_col]=find(E==Emin,1);        %sprawdzamy jej indeks
  65. h=V(:,Emin_col);                            %z wektora V wybieramy odpowiadajacy wektor wlasny
  66.                                             %to bedzie nasz wektor h
  67.                                              
  68. m=size(h);
  69. h1=h(1:(m/2));                              %dzielimy go sobie po polowie na 2 tak ze pol wektora to h1 a drugie pol to h2
  70.                                            
  71. h2=h(((m/2)+1):m);
  72. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%transponujemy Pp1 i Pp2
  73. Pp1_vec=Pp1';
  74. Pp2_vec=Pp2';
  75. Pp=[Pp1_vec;Pp2_vec];
  76.  
  77. %%%%%%%%%%%%%%%%%%%%%%%%%%%tutaj nizej z wektorow h tworzymy macierze toeplitza
  78. %%%%%%%%%%%%%%%%%%%%%%%%%%%oznaczone duzym H (wpiszcie w wikipedie toeplitz matrix to zobaczycie rozmieszczenie elementow)
  79. c1=zeros(N,1);
  80. c2=zeros(N,1);
  81. r1=zeros(1,N+L-1);
  82. r2=zeros(1,N+L-1);
  83.  
  84. for i=1:(m/2)
  85.     c1(i)=h1(i);                     %tu tworzymy pierwsze kolumny H1 i H2
  86.     c2(1)=h2(i);
  87. end
  88.  
  89.  
  90. r1(1)=c1(1);                        %tu tworzymy pierwsze wiersze H1 i H2
  91. r2(1)=c2(1);
  92. H1=toeplitz(c1,r1);                 %tworzenie macierzy H1 i H2
  93. H2=toeplitz(c2,r2);
  94.  
  95. H=[H1;H2];
  96. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%Metoda najmniejszych kwadratow
  97.  
  98. Nom=H'*Pp;                           %tworzymy licznik
  99. Denom=(H'*H);                        %tworzymy mianownik
  100. Pa=mldivide(Nom,Denom);              %dzielimy przez siebie
  101. Pa=Pa(1:N);                          %czasem trzeba uciąć sygnal do interesującej nas długości stąd ten zapis
  102. Pa_s=Pa*(mean(Pp1)/mean(Pa));        %skalowanie otrzymanego sygnału, żeby miał porownywalną wartość średnią do pomierzonych (wzór A6)
  103.  
  104. plot(t/2,Pa_s,'r');
  105.                        
  106.  
  107.  subplot(2,2,1);plot(t,Pp1);title('FIR1');
  108.  subplot(2,2,2);plot(t,Pp2);title('FIR2');
  109.  subplot(2,2,3);plot(t,X); title('original');
  110.  subplot(2,2,4);plot(t,Pa_s,'r');title('reconstructed');
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
Top