% Aparece una senal de frecuencia  fr aliada como una frecuencia menor
% dentro del  rango de interes  [0 fN0] debido  al fS0. Esta  senal es
% indeseable para  hacer unas mediciones.   Aqui se simula el  dato de
% interes como una cosinusoide de f1 a fines ilustrativos.
%
% Para solucionar este inconveniente se propone:
%
% a) Subir el fS  = fr. En este caso el ruido se  ve como una senal de
% frecuencia cero. Si  se remueve el valor de continua  de la senal se
% remueve el ruido.
%
% b) Subir fS = 2*fr. En este caso el ruido entra en el ancho de banda
% justo sobre Nyquist.   Como en el dato original  se utilizaba fs=100
% Hz,  quiere decir  que  se interpretaban  frecuencias  hasta 50  Hz.
% Ahora el ruido  se va por fuera  del ancho de banda  de interes.  Al
% aplicar un  pasa bajas  a la  senal puede  atenuarse la  energia del
% ruido en el ancho de banda de interes original [0-50] Hz.
%
%  fr = fo + k fS, observo fr  como fo debido a la existencia de fS =>
%  fS =  fr-fo, k=1  => fN  = (fr-fo)/2. Doy  un fo  y obtengo  fN. Es
%  decir, pido  como quisiera  ver a  fr y puedo  jugar con  ello para
%  buscar atenuar la presencia del ruido en la senal.
%
% Casos considerados
%
% fN = (fr-2*fN0)/2 => observo fr como de 2*fN0. En este caso el ruido
% cae por encima del ancho de banda de interes original [0 fN0]. De
% pasar un pasa bajas atenuo el ruido. El dt cae a la mitad del
% original, lo que supone que tomo al menos el doble de muestras,
% lo cual puede no ser tan drastico.
%
% fN  =  fr/2  =>  fS  =  fr, observo  el  ruido  como  una  senal  de
% continua. Quitando el  promedio a la  senal observada  atenuo el
% ruido. el dt es mucho mas chico que el original. Tomo mas muestras
%
% fN  = fr  => el  ruido  cae sobre  Nyquist, lejos  del intervalo  [0
% fN0]. De pasar un pasa bajas se logra atenuar el ruido en la senal.
%
clear all;close all;

f1  = 10   % [Hz] frecuencia del "dato"
fr  = 1020 % [Hz] frecuencia del "ruido"
fN0 = 50   % frecuencia de Nyquist del dato medido originalmente
T  = 1;    % [s]  tiempo maximo de observacion

% filtro de pasa bajas rectangular

nh = 31;
mh = (nh-1)/2;


for fN = [ fN0 (fr-2*fN0)/2 fr/2 fr ]    % Hz
    fN
    fS = 2*fN;  % frecuencia de muestreo
    dt = 1/fS; % intervalo de muestreo

    t=(0:dt:T);
    
    r  = 0.5*cos(2*pi*fr*t);  % ruido
    s1 = cos(2*pi*f1*t);      % senal 
    s  = s1+r;                % dato
    N  = length(s);
    
figure('Name',[' fN = ' num2str(fN) ' Hz'])

 if      fS==fN0
     sf  = s;
 elseif fS == fr
     sf  = s-mean(r);  % remuevo la amplitud de "continua"
 else
     % filtrar con pasa bajas rectangular
     wc = 2*pi*fN0*dt;
     h  = (wc/pi) * sinc(wc.*(-mh:mh)/pi);
     sf = conv(h,s)(mh+1:mh+1+N-1);
 end

subplot(4,1,1)
plot(t,s1,'k;senal;')
    axis([0 T -1.5 +1.5])
subplot(4,1,2)
plot(t,s,'b;medido;')
    axis([0 T -1.5 +1.5])
    subplot(4,1,3)
plot(t,r,'r;ruido;')
    axis([0 T -1.5 +1.5])
    subplot(4,1,4)
plot(t,sf,'g;filtrado;')
    axis([0 T -1.5 +1.5])
xlabel('t (s)')

Npad=5*length(s);
f=(0:Npad-1)*fS/Npad;

figure('Name',[' Espectros  fN = ' num2str(fN) ' Hz'])
subplot(2,1,1)
hold on
plot(f,(2/Npad)*abs(fft(s,Npad)),'b;medido;')
plot(f,(2/Npad)*abs(fft(r,Npad)),'r;ruido;')
hold off
axis([0 fN])
subplot(2,1,2)
hold on
if (fN == fN0)
plot(f,(2/Npad)*abs(fft(s,Npad)),'b;medido;')
plot(f,(2/Npad)*abs(fft(r,Npad)),'r;ruido;')
else
plot(f,(2/Npad)*abs(fft(s,Npad)),'b;medido;')
plot(f,(2/Npad)*abs(fft(r,Npad)),'r;ruido;')
plot(f,(2/Npad)*abs(fft(sf,Npad)),'g;filtrado;')
end
hold off
xlabel('f (Hz)')
axis([0 50])
end
