function xrec=muestreouniforme(N,M,d);
% Funcion muestreouniforme(N,M,d)
% Este programa genera, de forma aleatoria, una señal x(n) de tamaño N y ancho de banda M y, tras
% muestrearla uniformemente con paso d, la recupera utilizando el teorema del muestreo digital uniforme.
%
% Evidentemente, este programa se puede modificar de forma muy simple para que, en vez de generar una señal
% de modo aleatorio, tome dicha señal como dato y a dicha señal el teorema del muestreo (Ejercicio!).

% N=longitud de la señal
% M=ancho de banda
% d=espaciado de la muestra
% xrec = la salida
% %%%%%%%%%%%
% Generamos una señal aleatoria (con distribución normal) de tamaño N.
close all;
 y=randn(1,N);
% Usamos dicha señal para convertirla en una señal de banda limitada con ancho de banda <=M.
filt1=ones(1,M+1);
filt2=zeros(1,N-2*M-1);
filt3=ones(1,M);
filt=[filt1 filt2 filt3];
y=y.*filt;
x=1000*ifft(y);
% Calculamos las muestras que nos haran falta para aplicar el teorema del muestreo digital uniforme
r=(N/d)-1;
% Aplicamos el teorema del muestreo digital uniforme (Para ello, haremos uso de la función de Matlab
% discretsincvector.m, que mostramos implementada a continuación de este algoritmo).
sincvector=discretsincvector(M,N);
for t=0:1:N-1
    for j=0:1:r
    tau(j+1)=d*x((d*j)+1)*sincvector(t-d*j+N);
end
xrec(t+1)=sum(tau);
end

% Dibujamos los resultados logrados.

% Dibujamos la parte real de: la señal original, los valores muestrales, y la señal recuperada
z=1:1:N;
xp=0:d:N-1;
figure;
subplot(3,1,1); plot(z,real(x)); xlabel('Parte real de la señal');
subplot(3,1,2); stem(xp,real(x(xp+1)));xlabel('Valores muestrales');
subplot(3,1,3); plot(z,real(xrec),'r');xlabel('Parte real de la señal recuperada');
figure;
% Dibujamos la parte imaginaria de: la señal original, los valores muestrales, y la señal recuperada
subplot(3,1,1); plot(z,imag(x)); xlabel('Parte imaginaria de la señal');
subplot(3,1,2); stem(xp,imag(x(xp+1))); xlabel('Valores muestrales');
subplot(3,1,3); plot(z,imag(xrec),'r'); xlabel('Parte imaginaria de la señal recuperada');
% Dibujamos los errores cometidos tanto para la parte real como para la parte imaginaria
figure;
subplot(2,1,1); plot(z,real(x-xrec));xlabel('Parte real del Error');
subplot(2,1,2); plot(z,imag(x-xrec));xlabel('Parte imaginaria del Error');

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

function SM = discretsincvector(M,N)
% function SM = discretsincvector(M,N)
% Calcula el analogo discreto de la funcion seno cardinal para trabajar con señales de banda limitada
% (ancho de banda M) y tamaño N dados. Calcula dicha funcion entre -(N-1) y N-1, que es donde es útil
%
%  Parametros:
%   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;