function grochenig(func,dt,L,r,tol)
% function grochenig(func,dt,L,r,tol)
% Este programa, a partir de una señal analógica  f(t) (implementada en el fichero de nombre func.m)
% genera una señal digital tomando las muestras x(n)=f(n*dt), para n=0,1,...,N-1 y (N-1)dt<=L<Ndt.
% A continuación construye, de forma aleatoria, un  conjunto de puntos
%              0<=n_0<n_1<...<n_r<N
% distribuidos de forma irregular en {0,1,...,N-1} y calcula M con la propiedad de que
%                2M(2[delta/2]+1)<N
% y, por tanto, nuestro vector es valido para el muestreo de señales digitales
% de tamaño N y ancho de banda M
%
% Finalmente, utiliza los valores muestrales x(n_0),....,x(n_r) para recuperar la señal x(n) para todo n
% mediante la aplicación del Teorema de Grochenig (Teorema 5.4.10 del texto).
%
% Autores: J. M. Almira, J. A. Cid, A. E. Romero (2005)
%
close all;
h=0:dt:L;
x=feval(func,h);
N=length(x);
% elegimos r indices para muestra aleatoria (Para ello se usa una función de Matlab que implementamos a
% continuación de este algoritmo: la función muestrairreg.m)
xp=muestrairreg(N,r);
% calculamos los pesos correspondientes
wp=pesos(xp,N);
%calculamos el maximo ancho de banda M para ese r
t=2:1:r; xp2=[xp(t),xp(1)+N] ;
xp1=[xp(t-1),xp(r)];
delta = max(xp2-xp1)
Mmax =floor(N/(2*(2*floor(delta/2)+1)))
M=input('Introduzca el valor de M: ');
% Cálculo de la matriz matrix=lambda*wp*P
gamma=((sin(pi*M/N))/(sin(pi/(4*floor(delta/2) + 2))));

%%%% Cálculo del parámetro óptimo de relajación %%%%

lambda=2/((1+gamma)^2+(1-gamma)^2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% A continuación hacemos uso de la función discretsincvector.m, que ya habíamos hecho explícita en el
% programa de Matlab que implementa el teorema del muestreo digital uniforme. De todos modos, incluimos
% también ahora dicho fichero, por completitud (ver el final del algoritmo)
sinc=discretsincvector(M,N);
for i=1:r
    for j=1:N
        matrix(i,j)=lambda*wp(i)*sinc(j-1-xp(i)+N);
    end
end
    x0=x(xp)*matrix;
    x1=x0;

% Calculamos el número de máximo iteraciones (usando la estimación a priori del teorema)

gammabar=max(abs(1-lambda*(1+gamma)^2),abs(1-lambda*(1-gamma)^2))
nitermax=floor(log(tol*(1-gammabar)/(norm(x0)))/log(gammabar)+1)+1

% Realizamos las iteraciones pero introducimos como criterio de parada que
% ||x_{n+1}-x_{n}||.(gamma/(1-gamma))<=toletancia. Al final, se verá
% cuántas iteraciones (comparadas con lo estimado a priori) se han tenido
% que realizar.
    for i=1:1:nitermax
        x2=x1+x0-x1(xp)*matrix;
        if (norm(x2-x1) < tol*(1-gammabar)/gammabar )
           i
           x1=x2;
           break
        end
            x1=x2;
    end
figure(1);
subplot(4,1,1);plot(1:N,real(x),'r');zoom on; xlabel('Señal original');
subplot(4,1,2);stem(xp,real(x(xp)));zoom on; xlabel('Muestras tomadas');
subplot(4,1,3);plot(1:N,real(x1),'b');zoom on; xlabel('Señal recuperada');
subplot(4,1,4);plot(1:N,real(x - x1),'r');zoom on; xlabel('Error');


%%%%%%%%%%%%%%%%%%%%%%%
%
% Función auxiliar muestrairreg
%
%%%%%%%%%%%%%%%%%%%%%%%


function xp=muestrairreg(N,r);
% Construye de forma aleatoria, un  conjunto de puntos
%              0<=n_0<n_1<...<n_r<N
% distribuidos de forma irregular en {0,1,...,N-1} y calcula M con la propiedad de que
%                2M(2[delta/2]+1)<N
% y, por tanto, nuestro vector es valido para el muestreo de señales digitales de tamaño N
% y ancho de banda M
%
%  Parametros:
%   N: tamaño de la señal a muestrear
%   M: ancho de banda de la señal x.
%   r+1: numero de puntos de la muestra
%   xp: señal de salida
%
%  Autores: J. M. Almira, A. E. Romero (2005).
%

 d=floor(N/(r+1));
 xp=d:d:r*d;
 modif=-1+2*rand(1,r);
 modif=sign(modif);
 xp=xp+floor(d/3)*modif;


%%%%%%%%%%%%%%%%%%%%%%%
%
% Función auxiliar discretsincvector
%
%%%%%%%%%%%%%%%%%%%%%%%

function SM = discretsincvector(M,N)
% function SM = discretsincvector(M,N)
% Calcula el análogo discreto de la función seno cardinal para trabajar con señales de banda limitada
% (ancho de banda M) y tamaño N dados. Calcula dicha función entre -(N-1) y N-1, que es donde es útil
%
%  Parámetros:
%   M: ancho de banda
%   N: tamaño de la señal
%   SM: vector salida
%
%  Autores: J. M. Almira, A. E. Romero (2005).
%

% Vector de abcisas
t = (-N+1):(N-1);

% Reservamos memoria: la señal SM (salida) será un vector 1xN (todo a cero, inicialmente)
  SM = repmat(double(0), size(t));

% Calculamos la función muestreo para valores distintos de cero
% Ojo al direccionamiento con vectores

  SM( t>0 ) =  sin(pi*(2*M+1)*t(t>0)/N)./(N*sin(pi*t(t>0)/N));
  SM( t<0 ) =  sin(pi*(2*M+1)*(t(t<0)+N)/N)./(N*sin(pi*(t(t<0)+N)/N));

% Asignamos el valor (en el límite) de la funcion para t = 0

  SM( t==0 ) = (2*M+1)/N;
  
  %%%%%%%%%%%%%%%%%%%%%%%
%
% Función auxiliar pesos
%
%%%%%%%%%%%%%%%%%%%%%%%

function W = pesos(x,N);

r=length(x);
x=[x, x(1)+N];
L=0.*(1:1:(r+1));
W=0.*(1:1:r);
for k=1:1:r
    L(k)=floor((x(k)+x(k+1))./2)
end
L(r+1)=L(1)+N;
for k=1:1:r
    W(k)=L(k+1)-L(k)
end;