%navrh filtru FIR
%PSP '99
clear;

%parametry filtru
as=-40
ap=-2
AP=-10^(ap/20)
AS=-10^(as/20)

fs=600
fd=1000 		%dolni
fh=2000 		%horni
fq=2580

%pro spravne normovani
df=fh-fd;
fm1=sqrt(fd*fh);
fm2=sqrt(fq*fs);
if fm1<fm2 fm=fm1 ;
 else fm=fm2;
end

%vzorkovaci kmitocet
fv=8000

%transformace na NDP
%z ni potrebujeme v podstate jen delta omega pro spravny vypocet radu
omh=(fh^2-fm1^2)/(fh-fd)/fh;
omd=(fd^2-fm2^2)/(fh-fd)/fd;
oms=(fs^2-fm1^2)/df/fs;
omq=(fq^2-fm2^2)/df/fq;
if abs(oms)>=abs(omq);
 omegas=abs(omq)
 delom=(omegas-1)*1/((fq^2-fm2^2)/df/fq);
else
 omegas=abs(oms)
 delom=(omegas-1)*1/((fs^2-fm1^2)/df/fs);
end

%odhad radu filtru
N=6.2*pi/delom;
N=round(N)
%pro dany utlum vybereme podle literatury vhodne okno
b=fir1(N,[(fs+fd)/2/fv*2 (fh+fq)/2/fv*2],hanning(N+1));
a=1;
[H w]=freqz(b,a,2048);
ft1=(fs+fd)/2;
ft2=(fh+fq)/2;

%iterace - specificka pro tento filtr
iter=0;

while abs(H(round(fd/fv*2*2048)))>abs(AP),
	N=6.2*pi/delom;
	N=round(N);
	ft1=ft1+10;
   ft2=ft2-10;
   b=fir1(N,[(fs+fd)/2/fv*2 (fh+fq)/2/fv*2],hanning(N+1));
	[H w]=freqz(b,a,2048);
	while abs(H(round(fs/fv*2*2048)))>=abs(AS),
   	N=N+1;
		b=fir1(N,[ft1/fv*2 ft2/fv*2],hanning(N+1));
      [H w]=freqz(b,a,2048);
 	end
   iter=iter+1;
 	%rekneme ze 10 kroku staci, stejne to uz lepsi nebude
   if iter>10 break
   end
end

%konecny rad filtru a mezni frekvence
N
ft1
ft2

%vykresleni frekvencni charakteristiky jak ji poskytuje matlab
figure(1);
freqz(b,a,512);

%vykresleni |H(jOMEGA)|
figure(2);
%fv*w/(2*pi)
plot(w,abs(H));
xlabel('OMEGA')
ylabel('|H(jOMEGA)|')
title(['Tolerancni schema a |H(jOMEGA)|, FIR radu ',int2str(N)]);

grid;
%tolerancni pole
hold on;
FS=fs/fv*2*pi;
FD=fd/fv*2*pi;
FH=fh/fv*2*pi;
FQ=fq/fv*2*pi;
FV=pi;
plot([0 FS],[-AS -AS],'r');
plot([FS FS],[-AS 1],'r');
plot([FD FD],[-AS -AP],'r');
plot([FD FH],[-AP -AP],'r');
plot([FH FH],[-AP -AS],'r');
plot([FQ FQ],[1 -AS],'r');
plot([FQ FV],[-AS -AS],'r');
hold off;


%uprava koeficientu FIRu - scaling Q15
koefb=round(b.*32768);

%vypsani koeficientu ve tvaru pro assembler - da se i stahnout v okne matlabu
%zde ulozena pomoci fwrite
skoefb=sprintf(' .word %d\n',koefb)
fid=fopen('c:\xampl\skola\psp\sem\asm\noncirc\koefs','wb');
fwrite(fid,skoefb,'char');
fclose(fid);

%impulsova odezva
figure(4);
drc=[1 zeros(1,59)];

%jeji ulozeni
fid=fopen('c:\xampl\skola\psp\sem\asm\noncirc\input.hex','wb');
%+32767 ~ +1, -32768 ~ -1; lze pouze .*32767 protoze v drc je pouze jedna 1 jinak nuly
fprintf(fid,'%04x\n',drc.*32767);
fclose(fid);

%odezva filtru vypocitaneho
h=filter(b,1,drc);
subplot(211);
plot(h,'b');
title('Impulsova odezva vypocitana');
%hold on;

%nahrani vystupu FIRu emulatoru TMS
fid=fopen('c:\xampl\skola\psp\sem\asm\noncirc\output.hex','r');
outputh=fread(fid);
fclose(fid);
outputd=[];
f=1;
%pri save jeden prikaz, pri load cele pocitat 'matlabove' :(
while f<length(outputh),
   cislo=[outputh(f) outputh(f+1) outputh(f+2) outputh(f+3)];
   cislo=hex2dec(setstr(cislo));
   if cislo>32767 cislo=cislo-65536;
   end;
   outputd=[outputd cislo/32768];
   f=f+6;
end
%odezva filtru ziskana emulaci
subplot(212);
plot(outputd(1:60),'r');
title('Impulsova odezva ziskana ze simulatoru');
%hold off;