% puerto_serie.m
% macro de MatLab para trabajar con el puerto serie

clear;                          % borramos todas las variables

fichero_out_t='x_tiempo.txt';   % fichero de salida (la señal en el dominio del tiempo)
fichero_out_f='x_freq.txt';     % fichero de salida (la señal en el dominio de la frecuencia, espectro de potencia)
clave_num_mues=3;               %        0,    1,    2,    3,     4,     5,     6, o      7 
                                % para 500, 1250, 2500, 5000, 12500, 25000, 50000, o 125000 muestras
RBW=3000;                       % resolution-band-width (en Hz) utilizado para el espectro suavizado
timeout=10;                     % timeout en segundos para el puerto serie (subirlo si se quieren transferir muchas muestras)

% PRIMER PASO: configuración del puerto serie
s1=serial('com1','BaudRate',38400,'Parity','even','StopBits',1,'InputBufferSize',262144,'OutputBufferSize',4096,'Timeout',timeout);

% TEST: identificación del dispositivo (abrir,comunicar,cerrar)
%fopen(s1);
%fprintf(s1,'*IDN?');
%idn=fscanf(s1)
%fclose(s1);
%% si no se cierra el dispositivo y se intenta abrir de nuevo puede haber
%% problemas. Para cerrar todos los dispositivos conectados al puerto serie se usa lo siguiente:
obj=instrfind;  %% localiza todos los puertos serie abiertos o cerrados 
fclose(obj);    %% cierra los abiertos
close all;      %% tambien cerramos todas las figuras

% SEGUNDO PASO: configuramos la transmisión y leemos parámetros:
fopen(s1);
fprintf(s1,'*IDN?');                              % identificador del dispositivo
pause(0.20); idn=fscanf(s1);
pause(0.20); fprintf(s1,':RUN');                  % activamos el sistema de adquisición del osciloscopio
pause(0.20); fprintf(s1,':ACQuire:LENGth %d',clave_num_mues); % Configuramos el numero de muestras
pause(0.20); fprintf(s1,':ACQuire:LENGth?');
pause(0.20); longitud=fscanf(s1); 
fprintf('Sistema configurado para la adquisición de %d muestras\n',sscanf(longitud,'%d'));
pause(0.20); fprintf(s1,':ACQuire:MODe 0');        % modo 0 para "muestras"; los otros modos para "average" o "picos"
pause(0.20); fprintf(s1,':CHANnel:DISPlay 1');     % activa canal 1
pause(0.20); fprintf(s1,':MEASure:SOURce 1');      % selecciona canal 1 para medidas
pause(0.20); fprintf(s1,':STOP');                  % paramos el sistema de adquisición
pause(0.20); fprintf(s1,':MEASure:VRMS?');         % pide VRMS
pause(0.20); vrms=fscanf(s1);
pause(0.20); fprintf(s1,':MEASure:VPP?');          % pide VPP
pause(0.20); vpp=fscanf(s1);
fprintf('Medidas de amplitud proporcionadas por osciloscopio:\n');
fprintf('    Amplitud pico-pico:  %s    Amlpitud rms:        %s\n',vpp,vrms);
pause(0.20); fprintf(s1,':CHANnel1:SCALe?');       % leemos la escala de amplitud (necesario para obtener x(t))
pause(0.20); escala=fscanf(s1);
A=sscanf(char(escala),'%f %s');
escala_v=A(1);                     % valor
unidad_v=char(A(2:end));           % unidades (V, mV, uV)
switch unidad_v'
	case 'V',  factor_v=1;         % factor_v es el factor por el que hay que multiplicar para pasar a V
	case 'mV', factor_v=1e-3;
	case 'uV', factor_v=1e-6;
end
escala=escala_v*factor_v;          % obtenemos la escala en Voltios/division

% TERCER PASO: volcado de memoria con las muestras:
pause(0.20); fprintf(s1,':ACQuire1:MEMory?');      % solicitamos transferencia de datos
pause(0.20); dato=fread(s1,1,'uchar');             % leyendo de la cabecera: caracter especial
pause(0.20); dato=fread(s1,1,'uchar');             % leyendo de la cabecera: "digito de tamaño de datos" 
   dig_tam_datos=sscanf(char(dato),'%d');
pause(0.20); dato=fread(s1,dig_tam_datos,'uchar'); % leyendo de la cabecera: "tamaño de datos" que permite calcular nmues
   tam_datos=sscanf(char(dato),'%d');
   nmues=(tam_datos-8)/2;                          
pause(0.20); dato=fread(s1,4,'uchar');             % leyendo de la cabecera: la frecuencia de muestreo
   tmp_file=fopen('tmp_file','wb');
   fwrite(tmp_file,dato);
   fclose(tmp_file);
   tmp_file=fopen('tmp_file','rb');                % por problemas con ordenación de bytes, hay que escribir en fichero temporal
   freq_muestreo=fread(tmp_file,'float',1,'ieee-be');
   fclose(tmp_file);
pause(0.20); dato=fread(s1,4,'uchar');             % leyendo de la cabecera: canal y tres bytes de relleno
canal=char(dato(1));
pause(0.20); muestras=fread(s1,nmues,'int16');     % leyendo las muestras: se vuelve a usar el fichero temporal
   tmp_file=fopen('tmp_file','wb');
   fwrite(tmp_file,muestras,'int16');
   fclose(tmp_file);
   tmp_file=fopen('tmp_file','rb');
   muestras=fread(tmp_file,'int16','ieee-be');
   fclose(tmp_file);
fclose(s1);

% CUARTO PASO: asignamos valores de tiempo y amlitud a las muestras:
x=muestras*escala/25;                              % estableciendo la amplitud (cada division son 25 cuantos)
t=(1:length(x))/freq_muestreo;                     % estableciendo el tiempo para cada muestra
VRMS=sqrt(mean(x.*x.*hanning(length(x)))/mean(hanning(length(x))));  % calculamos amplitud rms
VDBM=20*log10(VRMS/sqrt(0.05));                    % calculamos amplitud dBm
fprintf('Medidas de amplitud (sobre señal adquirida):\n');
fprintf('    Amplitud pico-pico:  %f Vpp\n    Amlpitud rms:        %f Vrms\n',max(x)-min(x),VRMS);
fprintf('    Amplitud dBm:        %f dBm\n',VDBM);

% QUINTO PASO: representación gráfica de la señal, el espectro de potencia,
% y el espectro de potencia suavizado

figure(1)   % la señal
plot(t*1000,x*1000,'r');
xlabel('tiempo (ms)');
ylabel('Amplitud (mV)');
ylim([-escala*4000 escala*4000]);
title('Señal en el tiempo')
grid on;                            % pinta el grid
zoom on;                            % activa el zoom

figure(2);  % espectro          
L0=length(x);
Nbins=ceil(log(L0)/log(2));
L=2^(Nbins+2);                      % 4 veces potencia de 2 mayor que el número de muestras (para suavizar)
y=x.*hanning(L0);                   % previo a FFT se le aplica ventana de Hanning (para evitar rizado)
Y0=abs(fft(y,L))/sqrt(L);           % FFT
vrms_y0=sqrt(sum(Y0.*Y0)/L/2);      % amplitud rms sobre el espectro
alpha=VRMS/vrms_y0;                 % factor de escala para compensar ventana de hamming y tamaño de ventana fft
Y=alpha*Y0;                         % escalando
YY=Y.^2;                            % densidad de potencia espectral escalada
vrms_y=sqrt(sum(YY(1:L/2))/L);      % amplitud de la señal en voltios rms (igual a vrms_x)
PY=10*log10(YY/L0/0.05);            % densidad de potencia espectral en dBm (dividimos por L0 en lugar de L para que la altura
                                    %    de los picos sea la misma que si hubieramos hecho la FFT con LO puntos)
Ts=t(2)-t(1);                       % periodo de muestreo
freq=(0:(L/2))/L/Ts;                % vector de frecuencias (en Hz)
PY=PY(1:(L/2+1));                   % sólo nos interesan del espectro la muestra 1 (freq=0) hasta la L/2+1 (freq=fs/2)
                                    % el resto de puntos del espectro son redundantes por ser x real
plot(freq/1000,PY,'r');             % y dibujamos el espectro
xlabel('frecuencia (kHz)');
ylabel('potencia espectral (dBm)');
title('Espectro de potencia');
limites_x=xlim();
grid on;                            % pinta el grid
zoom on;                            % activa el zoom

figure(3); % espectro suavizado de acuerdo con el "resolution-band-width" RBW establecido
num_freq=sum(freq<RBW);
ventana=hanning(num_freq)*1.5;
ventana=ventana.*(ventana<1)+1*(ventana>=1);
YY_smooth=filter(ventana,1,YY);
PY_smooth=10*log10(YY_smooth/L/0.05);
PY_smooth=PY_smooth(1:(L/2+1));
freq_smooth=freq-0.5*RBW;
plot(freq_smooth/1000,PY_smooth,'r');
cadena1=sprintf('frecuencia (kHz)');
xlabel(cadena1);
xlim(limites_x);
cadena1=sprintf('potencia espectral (dBm)');
ylabel(cadena1);
cadena1=sprintf('Espectro de potencia suavizado: RBW=%.1fkHz',RBW/1000');
title(cadena1);
grid on;                            % pinta el grid
zoom on;                            % activa el zoom

% SEXTO PASO: escribiendo la señal y el espectro en sendos ficheros de salida
% Señal en el tiempo: 
t_out=t*1000000;            %    t_out es el tiempo en us  
x_out=x*1000;               %    x_out  es la amplitud en mV (sin aplicar ventana de Hanning)  
% Espectro de potencia (densidad de potencia espectral):
freq_out=freq/1000000;      %    freq_out es la frecuencia en MHz
PX_out=PY;                  %    PX_out es la densidad de potencia espectral en dBm (dB de mW a 50 Ohm)
% Escribiendo datos:
A=[t_out' x_out];A=A';      % datos preparados: tiempo en us - tension en mV
B=[freq_out' PX_out];B=B';  % datos preparados: freq. en MHz - dens. pot. espectral en dBm
fid=fopen(fichero_out_t,'wt'); fprintf(fid,'%f\t%f\n',A); fclose(fid);
fid=fopen(fichero_out_f,'wt'); fprintf(fid,'%f\t%f\n',B); fclose(fid);

%%%% FIN DE LA MACRO
