clear 
clc
clear a
clear M
clear b

// IR MOSTRANDO A MATRIZ M EM CADA CASO ABAIXO

//  Exemplo 5 ("muitas" amostras):
nx = 1000;  // numero de amostras
inti=-1;
ints=1;
Tinv = (1/(ints-inti));  // eh a frequencia fundamental

n=5;     //  "grau" do Pol. trigonometrico, max multiplo da freq. fund.

//  Exemplo 6 (poucas amostras):
// nx = 10; inti=-1; ints=1;  n=5; 

//  Exemplo 7 (poucas amostras, e grau muito alto, "degenerando" 
//  a matriz M. Pedir o valor de M  e observar que eh singular):
// nx = 10; inti=-1; ints=1;  n=8; 

// DEFININDO A FUNCAO A SER APROXIMADA E A ORDEM DE TRUNCAMENTO
deff('y=f(x)','y=0*(x<0)+1*(x>=0)');
//deff('y=f(x)','y=-x.*(x<0)+x.*(x>=0)');
//deff('y=f(x)','y=cos(2*%pi*Tinv*x)');
//deff('y=f(x)','y=sin(2*%pi*11*x)+sin(2*%pi*20*x)');
// deff('y=f(x)','y=cos(18*x)');
gridx = inti:(ints-inti)/nx:ints;
gridx = gridx(1:size(gridx,2)-1); // Eliminando a última amostra!

tic;

// DEFININDO A BASE (PHIs)
deff('y=f1(x)','y=1.*x.^0');
for k=1:n-1
deff(strcat(['y=f',string(k+1),'(x)']),strcat(['y=cos(2*%pi*Tinv*',string(k),'*x)']));
end
for k=1:n-1    //n+1:2*n-1
deff(strcat(['y=f',string(k+n),'(x)']),strcat(['y=sin(2*%pi*Tinv*',string(k),'*x)']));
end


// CALCULANDO O SISTEMA DE EQS NORMAIS
for i=1:2*n-1   
for j=1:2*n-1
 deff('y=fmult(x)',strcat(['y=f',string(i),'(x).*f',string(j),'(x)']));
//H(i,j)=intg(inti,ints,fmult,err,err);
  H(i,j)=sum(fmult(gridx));
end
end
//disp('Veja como eh a estrutura da matriz de Hilbert:')
//H


// CALCULANDO O LADO DIREITO DAS EQS. NORMAIS Ha=b
 deff('y=fmult(x)',strcat(['y=f(x).*f1(x)']));
 b(1)=sum(fmult(gridx));  
for i=2:2*n-1
 deff('y=fmult(x)',strcat(['y=f(x).*f',string(i),'(x)']));
 b(i)=sum(fmult(gridx));  
end

a = inv(H)*b;

toc
disp(toc())


// TRACANDO CURVAS:
cf=[];   // dados para curva de f
cp=[];   // dados para curva de P
for x = inti:(ints-inti)/(10*nx):ints;
cf = [ cf f(x) ];
// CALCULANDO P(x) = a(1)/2 + a(2)cos(2*%pi*x) + a(3)cos(2*%pi*2x)...
  soma = 0;
  soma = soma + a(1);
  for k=1:n-1
   soma = soma + a(k+1)*cos(2*%pi*Tinv*k*x);
  end
  for k=1:n-1
   soma = soma + a(k+n)*sin(2*%pi*Tinv*k*x);
  end
cp = [ cp soma ];
end
subplot(411)
plot(inti:(ints-inti)/(10*nx):ints,cf,'r+')
plot(inti:(ints-inti)/(10*nx):ints,cp,'b')

// trancando a curva apenas nos pontos do grid (onde ela 
//  foi projetada para ser aproximadora)
cp=[];
for x=inti : (ints-inti)/(nx) : ints
// CALCULANDO P(x) = a(1)/2 + a(2)cos(2*%pi*x) + a(3)cos(2*%pi*2x)...
  soma = 0;
  soma = soma + a(1);
  for k=1:n-1
   soma = soma + a(k+1)*cos(2*%pi*Tinv*k*x);
  end
  for k=1:n-1
   soma = soma + a(k+n)*sin(2*%pi*Tinv*k*x);
  end
cp = [ cp soma ];
end
subplot(412)
plot(inti : (ints-inti)/(nx) : ints,cp,'*')


subplot(413)
plot(a,'*') 

subplot(414)
c = - [nx*a(1); (nx/2)*(a(2:n)+%i*a(n+1:2*n-1))];  // o zero eh o coef de sin(0x)
                           // pois a(1) eh o coef de cos(0x)
[phi,db]=phasemag(c);
plot(db,'r') 


  // t = -100:1/sample_rate:100;
  // N=size(gridx,'*'); //number of samples
  // s=1*(gridx>=0);
  s = f(gridx);
 
tic
  y=fft(s);
 toc
disp(toc())

//  c
//  y'

