multiratewebinar_oct2005_spanish.m

Design of Multistage Decimators

%% Diseno de Decimadores Multietapa
% Esta demostracion ilustra como podemos obtener decimadores
% progresivamente mas eficientes dado un conjunto de especificaciones de
% diseno. Comenzamos con un diseno polifase de una sola etapa y luego
% pasamos a disenos multietapa. Posteriormente, permitimos que ocurra
% "aliasing" en la banda de transicion del filtro y disenamos un filtro
% Nyquist multietapa. Finalmente, usamos un filtro CIC para la primera
% etapa de un diseno multietapa para asi obtener la implementacion mas
% eficiente.
%
%   Author(s): Ricardo A. Losada
%   Copyright 2005 The MathWorks, Inc.


%% Especificaciones de Diseno
% Estas son las especificaciones que utilizaremos a traves de la
% demostracion:
M   = 15;      % Factor de decimacion
Fs  = 100e6;   % Frecuencia de muestreo de entrada 100 MHz 
Fc  = 3.125e6; % Frecuencia de corte 3.125 MHz
Fp  = 2.925e6; % Maxima frecuencia de paso  2.925 MHz
Fst = 3.325e6; % Minima frecuencia de rechazo 3.325 MHz
TW  = Fst-Fp;  % Ancho de banda de region de transicion 0.4 MHz
Ap  = 1;       % Ripple en banda de paso 1 dB (pico a pico)
Ast = 80;      % Atenuacion minima en banda de rechazo 80 dB

%% Creacion de un Objeto para Agrupar las Especificaciones
% Podemos generar un objeto que agrupe todas las especificaciones con el
% siguiente comando:
Hfd = fdesign.decimator(M,'lowpass','Fp,Fst,Ap,Ast',...
      Fp,Fst,Ap,Ast,Fs);
  
%% Diseno del Filtro de Decimacion
% Para disenar el filtro con una sola etapa, invocamos el comando 'design'
% y especificamos que queremos un diseno optimo de ripples iguales.
Hm_eq = design(Hfd,'equiripple');
info(Hm_eq)

%%
% El filtro que obtenemos tiene 642 coeficientes. Dado que utilizamos una
% implementacion polifase, el costo de implementar el decimador es 43
% multiplicaciones por cada muestra de entrada.
fvtool(Hm_eq)

%% Espectro de Senal Decimada
% Asumiendo que la senal original ocupa totalmente el intervalo de Nyquist
% de manera uniforme, el espectro de la senal filtrada adquiere la forma de
% la respuesta en frecuencia del filtro. Al submuestrear dicha senal
% filtrada, las replicas espectrales se posicionan de manera adyacente pero
% sin que haya superposicion entre ellas por encima de los 80 dB
% permitidos. De esta manera se evita un efecto "aliasing" significativo.
[H,f]=freqz(Hm_eq,8192,Fs,'whole');
figure
plot(f-Fs/2,20*log10(abs(fftshift(H))))
hold on,
plot(f-Fs/2+Fs/M,20*log10(abs(fftshift(H))),'r-')
plot(f-Fs/2-Fs/M,20*log10(abs(fftshift(H))),'k-')
grid on
xlabel(' Frecuencia (Hz)')
ylabel('dB')
axis([-1e7 1e7 -100 5])
title('Espectro de la senal decimada')
legend('Espectro de banda base','Replica positiva','Replica Negativa')

%% Creacion de un Modelo para Simulink
% Con un simple comando, podemos generar un modelo para Simulink del
% decimador que hemos disenado:
block(Hm_eq,'Destination','new')

%% Diseno en Multiples Etapas
% Ahora veremos como aumentar la eficiencia computacional disenando un
% filtro multietapa para las mismas especificaciones de diseno.
Hm_multi = design(Hfd,'multistage');
fvtool(Hm_multi)

%%
% El resultado de este diseno son dos filtros de decimacion conectados en
% serie. La particion del diseno en varias etapas reduce el costo
% computacional. Para este diseno, el costo es de 14.6 multiplicaciones por
% cada muestra de entrada (en promedio).

%% Visualizacion del Modelo en Simulink
% Para generar un modelo para Simulink de este nuevo diseno, utilizamos el
% mismo comando que hemos invocado anteriormente:
block(Hm_multi,'Destination','new')

%% Diseno de Filtros de Nyquist para Decimacion
% Dadas las especificaciones de diseno, es posible incrementar el factor de
% decimacion de manera que las replicas de la senal filtrada, luego del
% submuestreo, se superpongan unicamente en la banda de transicion del
% filtro. Con esto podemos incrementar la eficiencia del diseno, utilizando
% filtros de Nyquist, sin incurrir en "aliasing" en la banda de interes. A
% manera de ejemplo, si tuvieramos una senal de audio muestreada a una
% frecuencia muy alta y quisieramos reducir la frecuencia de muestreo a
% 44.1 kHz, podriamos permitir superposicion de replicas espectrales en la
% banda de frecuencia entre 20 kHz y 22.05 kHz. El "aliasing" que ocurra en
% dicha banda no es de importancia ya que esta no forma parte de la banda
% auditiva del oido humano.

M2   = 16; % Factor de decimacion incrementado
Band = M2; % La banda debe ser igual al factor de decimacion
Hfn  = fdesign.decimator(M2,'nyquist',Band,...
      'TW,Ast',TW,Ast,Fs);
  
%%
% Notese que la frecuencia de corte de estas especificaciones esta dada por
% Fs/(2*Band) que corresponde a la misma frecuencia de corte que utilizamos
% originalmente. De esta manera aseguramos que las especificaciones de
% diseno no han cambiado. Tambien podemos ver que en este caso no estamos
% especificando el ripple permitido en la banda de paso. Nos es posible con
% disenos de Nyquist especificar dicho ripple. Sin embargo, como veremos,
% el ripple resultante sera mucho menor al ripple requerido originalmente
% (1 dB pico a pico).
  
%% Diseno Nyquist Multietapa!
% Cuando invocamos el comando "multistage" en este caso, cada una de las
% etapas del filtro multietapa sera un filtro Nyquist (debido a que Hfn
% especifica que queremos disenar un decimador con filtros Nyquist).
Hn_multi = design(Hfn,'multistage');
fvtool(Hn_multi)

%%
% El resultado son 4 filtros de media banda conectados en serie. Cada
% filtro media banda es muy eficiente ya que (aproximadamente) la mitad de
% sus coeficientes tienen valor zero, luego solamente es necesario utilizar
% la mitad de la longitud del filtro. En este ejemplo el costo
% computacional de este diseno es de 10.31 multiplicaciones por cada
% muestra de entrada.

%% Espectro de Senal Decimada con Filtros Nyquist
% De manera similar al caso anterior, podemos ver la superposicion de las
% replicas espectrales de la senal decimada, en este caso utilizando filtros
% Nyquist. Como ya hemos dicho, en este caso, las replicas se superponen
% mas que en el caso anterior. Sin embargo, no existe superposicion por
% encima de los 80 dB permitidos en regiones distintas a las bandas de
% transicion. Esto demuestra que la banda de interes no incurre en
% "aliasing" significativo.

[H,f]=freqz(Hn_multi,8192,Fs,'whole');
figure
plot(f-Fs/2,20*log10(abs(fftshift(H))))
hold on,
plot(f-Fs/2+Fs/M2,20*log10(abs(fftshift(H))),'r-')
plot(f-Fs/2-Fs/M2,20*log10(abs(fftshift(H))),'k-')
grid on
xlabel(' Frecuencia (Hz)')
ylabel('dB')
axis([-1e7 1e7 -100 5])
title('Espectro de la senal decimada con filtros Nyquist')
legend('Espectro de banda base','Replica positiva','Replica Negativa')

%% Diseno en multiples etapas utilizando filtros CIC
% Finalmente, para obtener un decimador extremadamente eficiente, vamos a
% utilizar un filtro CIC en la primera etapa de un diseno multietapa.
% Primero disenaremos un decimador con dos etapas y veremos que obtendremos
% una eficiencia comparable al diseno con filtros Nyquist. Luego
% cambiaremos a un diseno con tres etapas y asi obtendremos el diseno mas
% eficiente de todos los que hemos realizado en esta demostracion.

%% Diseno del Filtro CIC
% El factor global de decimacion es 16. Escogemos un factor de decimacion
% de 4 para el filtro CIC dejando un factor de 4 para las etapas restantes.

M1   = 4;
D    = 1; % Differential delay
Hd1  = fdesign.decimator(M1,'cic',D,Fp,Ast,Fs);
Hcic = design(Hd1,'multisection');

%% Diseno de la Segunda Etapa: Filtro de Compensacion
% El filtro de compensacion se debe disenar de acuerdo a las
% especificaciones del filtro CIC a compensar. Por ello, aqui utilizamos el
% numero de secciones del filtro CIC asi como el differential delay como
% parametros de diseno.

M2    = 4;
Nsecs = Hcic.NumberOfSections;
Hd2   = fdesign.decimator(M2,'ciccomp',D,Nsecs,Fp,Fst,Ap,Ast,Fs/M1);
Hcomp = design(Hd2,'equiripple');
Hcas = cascade(Hcic,Hcomp); % conexion en serie de los dos filtros

%%
% Este diseno requiere de 10.44 multiplicaciones por muestra de entrada. Lo
% que no representa una mejora en comparacion al diseno Nyquist multietapa.
% Sin embargo, podemos mejorar esto aun mas utilizando tres etapas en total
% en lugar de dos.

%% Diseno Alterno, Segunda Etapa: Filtro de Compensacion
% Decidimos disenar un filtro de compensacion con un factor de decimacion
% de dos.

M3     = 2;
Hd3    = fdesign.decimator(M3,'ciccomp',D,...
    Nsecs,Fp,Fs/(M1*M3)-Fp-TW,Ap,Ast,Fs/M1);
Hcomp2 = design(Hd3,'equiripple');

%% Diseno de la Tercera Etapa
% Fijense que hemos disenado un filtro de decimacion con un factor de dos
% en lugar de un factor de decimacion de 4. Luego tenemos que disenar un
% filtro mas para completar el diseno. Elegimos utilizar un filtro de media
% banda para la ultima etapa (generalmente la que contiene mas
% coeficientes) dada su eficiencia mencionada anteriormente:

M4    = 2;
Hd4   = fdesign.decimator(M4,'halfband',...
      TW,Ast,Fs/(M1*M3));
Hhalf = design(Hd4,'equiripple');
Hcas2 = cascade(Hcic,Hcomp2,Hhalf);

%%
% El costo computacional del diseno en tres etapas es solamente de 6.06
% multiplicaciones por muestra de entrada.

%% Analisis de las Etapas
% Si superponemos la magnitud de la respuesta en frecuencia de cada una de
% las tres etapas que hemos disenado, podemos ver como logramos que la
% respuesta global del decimador cumpla con las especificaciones de diseno.
% El filtro CIC opera a la frecuencia de la senal original, 100 MHz. El
% filtro de compensacion opera a un cuarto de esa frecuencia ya que el CIC
% tiene un factor de decimacion de 4. Finalmente, el filtro de media banda
% opera a la mitad del filtro de compensacion, debido a que el filtro de
% compensacion proporciona decimacion por un factor de dos. 
Hcic_norm  = cascade(Hcic,dfilt.scalar(1/gain(Hcic)));
hfvt = fvtool(Hcic_norm,Hcomp2,Hhalf,'ShowReference','off','Fs',[100,25,12.5]);
legend(hfvt,'Primera etapa: CIC','Segunda etapa: compensacion',...
    'Tercera etapa: halfband');

%%
% Notese como es la ultima etapa la que proporciona la transicion requerida
% entre 2.925 MHz y 3.325 MHz. Sin embargo, como este filtro unicamente
% opera a 12.5 MHz, no puede eliminar bandas de frecuencia en la senal
% alrededor de los multiplos de su frecuencia de operacion. Algunas de
% estas bandas son eliminadas por el filtro de compensacion que opera a 25
% MHz. Sin embargo, este filtro a su vez no puede remover las bandas
% centradas alrededor de multiplos de 25 MHz. Ahi es donde el filtro CIC
% proporciona la atenuacion deseada de manera que en la banda de rechazo,
% la atenuacion minima proporcionada por el conjunto total de filtros es 80
% dB, cumpliendo asi las especificaciones de diseno.

%% Analisis del Diseno Total
% Podemos verificar que el diseno cumple las especificaciones, mostrando la
% magnitud de la respuesta en frecuencia del diseno global.

hfvt = fvtool(Hcas2,'ShowReference','off','Fs',Fs, ...
 'NormalizeMagnitudeto1','on');

%% Generacion de un Modelo Simulink a partir de los Filtros Disenados
% De la misma manera en que lo hemos hecho anteriormente, con un solo
% comando podemos generar un modelo de Simulink para nuestro diseno en tres
% etapas.
block(Hcas2,'Destination','new')

%% Generacion de codigo VHDL 
% Para generar codigo VHDL correspondiente a nuestro diseno, tenemos que
% primero representar cada filtro con aritmetica de coma fija (o punto
% fijo). 
set(Hcas2.stage(:),'Arithmetic','fixed');
generatehdl(Hcas2)