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

nmod=0:N-1; %modyfikator
x=1:N; %oś y
nn1=0:N-1; %oś x1 ile probek przed FFT
nn2=0:N-1; %oś x2 wielkośc częstotliwości po FFT

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
length(nmod)
length(nn1)
pause

for n3=0:(length(nn1)-1)
%x(n3+1) =cos(2*pi*2*n3*t2); %przykładowe wartości próbek2 genercja sygnału dla cosinusa
x(n3+1) =(sin(2*pi*1000*n3*t1)+0.5*sin(2*pi*2000*n3*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
end
pause


% zmienic 1 lub 2
tmp1=f1/N % czestotliwosc 1 skoku
tmp2=f2/N % czestotliwosc 1 skoku

for n3=1:N;
nn2(n3)=tmp1*(n3-1); %oś X2 numeracja indeksów od 1, przeskok o czestotliwosc tmp1
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)
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
y(nr+1)=x(ncopy+1); % skopiuj we właśiwe miejsce
end
x=y; % podstaw wynik pod x
end

%przestawienie kolejnośći próbek: v2
if(typBitReverse==2)
a=1;
for b=1:N-1
if(b T=x(a);x(a)=x(b);x(b)=T;
end
c=N/2;
while (c a=a-c; c=c/2;
end
a=a+c;
end
end


% właściwe FFt v1 z ksiązki Tomasza Zielińskiego rozdział 9 algorytmy wyznaczania dyskretnej transformacji Fouriera str 249
if(typFFT==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); %%%% % skopiowanie górnej próbki
xdol=x(d)*W^(m-1); % korekta dolnej próbki
x(g)=xgora+xdol; % nowa górna próbka:górna plus dolna po korekcie
x(d)=xgora-xdol; % nowa dolna próbka:górna minus dolna po korekcie
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)
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
d=g+M; % g-górny d - dolny indeks próbki mtylka
T=x(d)*Wi; % "serce" FFT
x(d)=x(g)-T; % nowa dolna próbka:górna minus "serce"
x(g)=x(g)+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
end

% porównaj w matlabem
xc=fft(xc);
tekst = '....................................................'
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;
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
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(nn2,Xmag2,'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