% W matlab wstępny program do filtrowania muzyki *rewertynPL*
clear all;
%genercja sygnału
N=2^5 % liczba próbek potęga 2

nmod=0:N-1; % modyfikator
f1=8000; % 8000 próbek/s dla przykladu 1
f2=400; % 400 próbek/s dla przykladu 2
t1=1/f1; % o ile sekund przeskok pomiędzy próbkami
t2=1/f2; % o ile sekund przeskok pomiędzy próbkami

ile_rund_w_sek=ceil(f1/N) % zaokrąglenie w górębraki usupełnić zerami floor i ceil , round i fix
ile_rund_w_sek=2

for n3=1:90
tmp30(n3,1)=n3;
tmp31(n3,1)=n3;
end

nn1=0:(N*ile_rund_w_sek-1); % oś x1 ile probek przed FFT
nn2=0:(N*ile_rund_w_sek-1); % oś x2 wielkośc częstotliwości po FFT
nn3=0:(N*ile_rund_w_sek-1); % oś x2 wielkośc częstotliwości po FFT
x=1:(N*ile_rund_w_sek); % oś y
for n3=0:(N*ile_rund_w_sek-1)
x(n3+1)=0;
end

x
length(ile_rund_w_sek)
length(x)
length(nn1)
length(nn2)
pause

for nr_rundy=1:ile_rund_w_sek
for n3=0:(N-1)
nmmn(n3+((N)*(nr_rundy-1))+1)=n3+((N)*(nr_rundy-1))+1;
end
end
nmmn;
length(nmmn)
%pause

for nr_rundy=1:ile_rund_w_sek
for n3=0:(N-1)
%x(n3+((N)*(nr_rundy-1))+1) =cos(2*pi*2*n3*t2); %przykładowe wartości próbek2 genercja sygnału dla cosinusa
n3
x(n3+((N)*(nr_rundy-1))+1) =(sin(2*pi*1000*(n3+N*nr_rundy+0)*t1)+0.5*sin(2*pi*2000*(n3+N*nr_rundy+0)*t1+3/4*pi)); %przykładowe wartości próbek1
% n3 nr próbki, 1000 Hz i 2000 Hz, n3*t1 z której czeci sekundy próbka wzor z ksiązki L.Lyonsa str 64
end
end
%pause
x

% zmienic 1 lub 2
tmp10=f1/N % czestotliwosc 1 skoku
tmp11=f1/(N*1) % czestotliwosc 1 skoku
tmp2=f2/N % czestotliwosc 1 skoku

for n3=1:N;
nn2(n3)=tmp10*(n3-1); %oś X2 numeracja indeksów od 1, przeskok o czestotliwosc tmp10
end
for n3=1:N*ile_rund_w_sek;
nn3(n3)=tmp11*(n3-1); %oś X2 numeracja indeksów od 1, przeskok o czestotliwosc tmp11
end
plot(nn1,x,'--rs','LineWidth',2,...
'MarkerEdgeColor','k',...
'MarkerFaceColor','g',...
'MarkerSize',10)
pause


%.............
typBitReverse=2; % 1 lub 2 -wybór algorytmu przestawianai próbek
typFFT=2; % 1 lub 2 wybór właściwej pętli FFT
xc=x; % kopia sygnału x do obliczeni z gotowej funkcji FFT
%.............


%przestawienie kolejnośći próbek: v1
if(typBitReverse==1)
tekst = '....................................................'
for n3=1:99
tmp33(n3,1)=n3;
ii=1;
end
y=1:(N*ile_rund_w_sek); % wielkość y
for nr_rundy=1:ile_rund_w_sek
blok=N*(nr_rundy-1)
MSB=log2(N); % liczba bitów numerów próbek
for n=0:N-1; % kolejne próbki
ncopy=n;
nr=0; % stary numer próbki
for m=1:MSB % po wszystkich bitach
if (rem(n,2)==0) % czy jedynka na LSB
n=n/2; % jeśli nie przesuń w prawo
else
nr=nr+2^(MSB-m);% dodaj 2^(MSB-m)
n=(n-1)/2; % odejmni jedynkę przesuń w prawo
end
end
tmp33(ii,2)=nr+1+blok;
tmp33(ii,3)=ncopy+1+blok;
ii=ii+1;

y(nr+1+blok)=x(ncopy+1+blok); % skopiuj we właśiwe miejsce
end
x(1+blok:N+blok)=y(1+blok:N+blok); % podstaw wynik pod x tylko seriami po 1 rundzie bo nadpisze wszystkie orginalne dane
tmp33 % wynik kolejnosc sortowania
x % wynik po przeniesieniu bitowym
end
%pause
end

%przestawienie kolejnośći próbek: v2
if(typBitReverse==2)
tekst = '....................................................'
for nr_rundy=1:ile_rund_w_sek
blok=N*(nr_rundy-1);
a=1;
for b=1:N-1

if(b T=x(a+blok);x(a+blok)=x(b+blok);x(b+blok)=T;
end
c=N/2;
while (c a=a-c; c=c/2;
end
a=a+c;
end
end
x % wynik po przeniesieniu bitowym
%pause
end


% właściwe FFt v1 z ksiązki Tomasza Zielińskiego rozdział 9 algorytmy wyznaczania dyskretnej transformacji Fouriera str 249
if(typFFT==1)
tekst = '....................................................'
for nr_rundy=1:ile_rund_w_sek
blok=N*(nr_rundy-1);
for e=1:log2(N) % kolejne etapy
SM=2^(e-1); % szerokość motylka
LB=N/(2^e); % liczba bloków
LMB= 2^(e-1); % liczba motylków w bloku
OMB=2^e; % odległośc mędzy blokami
W=exp(-j*2*pi/2^e); % podstawa bazy fouriera
for b=1:LB % kolejne bloki
for m=1:LMB % kolejne motylki
g=(b-1)*OMB+m; % indeks górnej próbki motylka
d=g+SM; % indeks dolnej próbki motylka
xgora=x(g+blok); % skopiowanie górnej próbki
xdol=x(d+blok)*W^(m-1); % korekta dolnej próbki
x(g+blok)=xgora+xdol; % nowa górna próbka:górna plus dolna po korekcie
x(d+blok)=xgora-xdol; % nowa dolna próbka:górna minus dolna po korekcie
end
end
end
end
end

% właściwe FFt v2 z ksiązki Tomasza Zielińskiego rozdział 9 algorytmy wyznaczania dyskretnej transformacji Fouriera str 249
if(typFFT==2)
tekst = '....................................................'
ii=1
bb=1
for nr_rundy=1:ile_rund_w_sek
blok=N*(nr_rundy-1);
for e=1:log2(N) % kolejne etapy
L=2^e; % długośc bloków DFT, przesunięcie bloków
M=2^(e-1); % liczba motylków w bloku.szerokośc każdego motylka
Wi=1; % Wl=W.=1 startowa wartoś wsp bazy w etapie
W=cos(2*pi/L)-j*sin(2*pi/L); % mnożnik bazy fouriera Wl=W.
for m=1:M % kolejne motylki
for g=m:L:N % w kolejnych blokach od:krok:do
d=g+M; % g-górny d - dolny indeks próbki mtylka

length(x)
if blok == 0
tmp30(ii,2)=g+blok;
tmp30(ii,3)=d+blok;
tmp30(ii,4)=m;
tmp30(ii,5)=x(g+blok);
tmp30(ii,6)=x(d+blok);
if d+blok+2<(N*ile_rund_w_sek-1)
tmp30(ii,7)=x(g+blok+2);
tmp30(ii,8)=x(d+blok+2);
end
ii=ii+1;
else
tmp31(bb,2)=g+blok;
tmp31(bb,3)=d+blok;
tmp31(bb,4)=m;
tmp31(bb,5)=x(g+blok);
tmp31(bb,6)=x(d+blok);
if d+blok+2<(N*ile_rund_w_sek-1)
tmp31(bb,7)=x(g+blok+2);
tmp31(bb,8)=x(d+blok+2);
end
bb=bb+1;
end

T=x(d+blok)*Wi; % "serce" FFT
x(d+blok)=x(g+blok)-T; % nowa dolna próbka:górna minus "serce"
x(g+blok)=x(g+blok)+T; % nowa górna próbka:górna plus "serce"
end
Wi=Wi*W; % kolejna wartośc bazy fouriera Wl=W.
end % koniec pętli motylków
end % koniec pętli etapów
nr_rundy
blok
end
end


% porównaj w matlabem
fragment=N/ile_rund_w_sek % do potegi 2
xc(1:N)=fft(xc(1:N));
%xc=fft(xc);
tekst = '....................................................'
%format short g
blok;
tmp30;
tmp31;

pause
x
xc

blad_real = max(abs(real(x-xc))) % odejmowanie wyniku naszej funkcji z funkcja matlaba
blad_imag = max(abs(imag(x-xc))) % odejmowanie wyniku naszej funkcji z funkcja matlaba

% wykres magnitudy w widmie
for n3=1:N*ile_rund_w_sek;
Xmag1(n3)=sqrt(real(x(n3))^2+imag(x(n3))^2);

Xmag2(n3)=sqrt(real(xc(n3))^2+imag(xc(n3))^2);
end

Xmag1
Xmag2
length(Xmag1)
length(Xmag2)
pause

bar(nn2,Xmag1,'r') % plot,set,http://www.mathworks.com/help/matlab/ref/bar.html#bthxceu
title('Group1')
% plot(nn2,Xmag1,'--rs','LineWidth',2,...
% 'MarkerEdgeColor','k',...
% 'MarkerFaceColor','g',...
% 'MarkerSize',10)

pause

bar(nn3(1:N),Xmag2(1:N),'g') % plot,set,http://www.mathworks.com/help/matlab/ref/bar.html#bthxceu
title('Group2')
% plot(nn2,Xmag2,'--rs','LineWidth',2,...
% 'MarkerEdgeColor','k',...
% 'MarkerFaceColor','g',...
% 'MarkerSize',10)
%http://michalbereta.pl/dydaktyka/KPO/FFT_Lab.pdf
pause
close all