diff --git a/matlab/UnivariateSpectralDensity.m b/matlab/UnivariateSpectralDensity.m
index 95e8dab69..22fee8739 100644
--- a/matlab/UnivariateSpectralDensity.m
+++ b/matlab/UnivariateSpectralDensity.m
@@ -1,11 +1,25 @@
-function [omega,f] = UnivariateSpectralDensity(dr,var_list)
+function [oo_] = UnivariateSpectralDensity(M_,oo_,options_,var_list)
 % This function computes the theoretical spectral density of each
 % endogenous variable declared in var_list. Results are stored in 
-% oo_ and may be plotted. Plots are saved into the graphs-folder.
+% oo_.SpectralDensity and may be plotted. Plots are saved into the 
+% graphs-folder. 
 % 
+% INPUTS
+%   M_                  [structure]    Dynare's model structure
+%   oo_                 [structure]    Dynare's results structure
+%   options_            [structure]    Dynare's options structure
+%   var_list            [integer]      Vector of indices for a subset of variables.
+%    
+% OUTPUTS
+%   oo_                 [structure]    Dynare's results structure,
+%                                       containing the subfield
+%                                       SpectralDensity with fields freqs
+%                                       and density, which are of size nvar*ngrid.
+%
+
 % Adapted from th_autocovariances.m.  
 
-% Copyright (C) 2006-2013 Dynare Team
+% Copyright (C) 2006-2015 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -22,10 +36,6 @@ function [omega,f] = UnivariateSpectralDensity(dr,var_list)
 % You should have received a copy of the GNU General Public License
 % along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 
-global options_ oo_ M_
-
-
-omega = []; f = [];
 
 if options_.order > 1
     disp('UnivariateSpectralDensity :: I Cannot compute the theoretical spectral density') 
@@ -33,12 +43,7 @@ if options_.order > 1
     disp('Please set order = 1. I abort')
     return
 end
-pltinfo  = options_.SpectralDensity.plot; 
-cutoff   = options_.SpectralDensity.cutoff; 
-sdl      = options_.SpectralDensity.sdl; 
-omega    = (0:sdl:pi)';
-GridSize = length(omega);
-exo_names_orig_ord  = M_.exo_names_orig_ord;
+
 if isoctave
     warning('off', 'Octave:divide-by-zero')
 else
@@ -60,22 +65,19 @@ for i=1:nvar
         ivar(i) = i_tmp;
     end
 end
-f = zeros(nvar,GridSize);
-ghx = dr.ghx;
-ghu = dr.ghu;
+
+ghx = oo_.dr.ghx;
+ghu = oo_.dr.ghu;
 nspred = M_.nspred;
 nstatic = M_.nstatic;
-kstate = dr.kstate;
-order = dr.order_var;
+kstate = oo_.dr.kstate;
+order = oo_.dr.order_var;
 iv(order) = [1:length(order)];
 nx = size(ghx,2);
 ikx = [nstatic+1:nstatic+nspred];
-A = zeros(nx,nx);
 k0 = kstate(find(kstate(:,2) <= M_.maximum_lag+1),:);
 i0 = find(k0(:,2) == M_.maximum_lag+1);
 i00 = i0;
-n0 = length(i0);
-A(i0,:) = ghx(ikx,:);
 AS = ghx(:,i0);
 ghu1 = zeros(nx,M_.exo_nbr);
 ghu1(i0,:) = ghu(ikx,:);
@@ -91,87 +93,72 @@ for i=M_.maximum_lag:-1:2
     AS(:,j1) = AS(:,j1)+ghx(:,i1);
     i0 = i1;
 end
-Gamma = zeros(nvar,cutoff+1);
-[A,B] = kalman_transition_matrix(dr,ikx',1:nx,M_.exo_nbr);
+
+[A,B] = kalman_transition_matrix(oo_.dr,ikx',1:nx,M_.exo_nbr);
 [vx, u] =  lyapunov_symm(A,B*M_.Sigma_e*B',options_.lyapunov_fixed_point_tol,options_.qz_criterium,options_.lyapunov_complex_threshold,[],[],options_.debug);
 iky = iv(ivar);
 if ~isempty(u)
     iky = iky(find(any(abs(ghx(iky,:)*u) < options_.Schur_vec_tol,2)));
-    ivar = dr.order_var(iky);
+    ivar = oo_.dr.order_var(iky);
 end
+
+iky = iv(ivar);  
 aa = ghx(iky,:);
 bb = ghu(iky,:);
-
-if options_.hp_filter == 0
-    tmp = aa*vx*aa'+ bb*M_.Sigma_e*bb';
-    tmp(abs(tmp) < 1e-12) = 0;
-    Gamma(:,1) = diag(tmp);
-    vxy = (A*vx*aa'+ghu1*M_.Sigma_e*bb');
-    tmp = aa*vxy;
-    tmp(abs(tmp) < 1e-12) = 0;
-    Gamma(:,2) = diag(tmp);
-    for i=2:cutoff
-        vxy = A*vxy;
-        tmp = aa*vxy;
-        tmp(abs(tmp) < 1e-12) = 0;
-        Gamma(:,i+1) = diag(tmp);
-    end
-else
-    iky = iv(ivar);  
-    aa = ghx(iky,:);
-    bb = ghu(iky,:);
+ngrid = options_.hp_ngrid; %number of grid points
+freqs = (0 : pi/(ngrid-1):pi)'; % grid on which to compute
+tpos  = exp( sqrt(-1)*freqs); %positive frequencies
+tneg  = exp(-sqrt(-1)*freqs); %negative frequencies
+if options_.hp_filter == 0 && ~options_.bandpass.indicator %do not filter
+   filter_gain=ones(ngrid,1);
+elseif ~(options_.hp_filter == 0 && ~options_.bandpass.indicator) && options_.bandpass.indicator %filter with bandpass
+    filter_gain = zeros(1,ngrid);
+    lowest_periodicity=options_.bandpass.passband(2);
+    highest_periodicity=options_.bandpass.passband(1);
+    highest_periodicity=max(2,highest_periodicity); % restrict to upper bound of pi
+    filter_gain(freqs>=2*pi/lowest_periodicity & freqs<=2*pi/highest_periodicity)=1;
+    filter_gain(freqs<=-2*pi/lowest_periodicity+2*pi & freqs>=-2*pi/highest_periodicity+2*pi)=1;
+elseif ~(options_.hp_filter == 0 && ~options_.bandpass.indicator) && ~options_.bandpass.indicator %filter with HP-filter
     lambda = options_.hp_filter;
-    ngrid = options_.hp_ngrid;
-    freqs = 0 : ((2*pi)/ngrid) : (2*pi*(1 - .5/ngrid)); 
-    tpos  = exp( sqrt(-1)*freqs);
-    tneg  =  exp(-sqrt(-1)*freqs);
-    hp1 = 4*lambda*(1 - cos(freqs)).^2 ./ (1 + 4*lambda*(1 - cos(freqs)).^2);
-    mathp_col = NaN(ngrid,length(ivar)^2);
-    IA = eye(size(A,1));
-    IE = eye(M_.exo_nbr);
-    for ig = 1:ngrid
-        f_omega  =(1/(2*pi))*( [(IA-A*tneg(ig))\ghu1;IE]...
-                               *M_.Sigma_e*[ghu1'/(IA-A'*tpos(ig)) IE]); % state variables
-        g_omega = [aa*tneg(ig) bb]*f_omega*[aa'*tpos(ig); bb']; % selected variables
-        f_hp = hp1(ig)^2*g_omega; % spectral density of selected filtered series
-        mathp_col(ig,:) = (f_hp(:))';    % store as matrix row
-                                                 % for ifft
-    end;
-    imathp_col = real(ifft(mathp_col))*(2*pi);
-    tmp = reshape(imathp_col(1,:),nvar,nvar);
-    tmp(abs(tmp)<1e-12) = 0;
-    Gamma(:,1) = diag(tmp);
-    sy = sqrt(Gamma(:,1));
-    sy = sy *sy';
-    for i=1:cutoff-1
-        tmp = reshape(imathp_col(i+1,:),nvar,nvar)./sy;
-        tmp(abs(tmp) < 1e-12) = 0;
-        Gamma(:,i+1) = diag(tmp);
-    end
+    filter_gain = 4*lambda*(1 - cos(freqs)).^2 ./ (1 + 4*lambda*(1 - cos(freqs)).^2);    
 end
 
-H = 1:cutoff;
+mathp_col = NaN(ngrid,length(ivar)^2);
+IA = eye(size(A,1));
+IE = eye(M_.exo_nbr);
+for ig = 1:ngrid
+    f_omega  =(1/(2*pi))*( [(IA-A*tneg(ig))\ghu1;IE]...
+                           *M_.Sigma_e*[ghu1'/(IA-A'*tpos(ig)) IE]); % state variables
+    g_omega = [aa*tneg(ig) bb]*f_omega*[aa'*tpos(ig); bb']; % selected variables
+    f_hp = filter_gain(ig)^2*g_omega; % spectral density of selected filtered series
+    mathp_col(ig,:) = (f_hp(:))';    % store as matrix row
+end;
+
+f = zeros(nvar,ngrid);
 for i=1:nvar
-    f(i,:) = Gamma(i,1)/(2*pi) + Gamma(i,H+1)*cos(H'*omega')/pi;
+    f(i,:) = real(mathp_col(:,(i-1)*nvar+i)); %read out spectral density
 end  
 
+oo_.SpectralDensity.freqs=freqs;
+oo_.SpectralDensity.density=f;
+
 if isoctave
     warning('on', 'Octave:divide-by-zero')
 else
     warning on MATLAB:dividebyzero
 end
 
-if pltinfo
+if options_.nograph == 0
     if ~exist(M_.fname, 'dir')
         mkdir('.',M_.fname);
     end
-    if ~exist([M_.fname '/graphs'])
+    if ~exist([M_.fname '/graphs'],'dir')
         mkdir(M_.fname,'graphs');
     end
 
     for i= 1:nvar
         hh = dyn_figure(options_,'Name',['Spectral Density of ' deblank(M_.endo_names(ivar(i),:)) '.']);
-        plot(omega,f(i,:),'-k','linewidth',2)
+        plot(freqs,f(i,:),'-k','linewidth',2)
         xlabel('0 \leq \omega \leq \pi')
         ylabel('f(\omega)')
         box on
diff --git a/matlab/stoch_simul.m b/matlab/stoch_simul.m
index d6b6ceadf..ff8b53cfd 100644
--- a/matlab/stoch_simul.m
+++ b/matlab/stoch_simul.m
@@ -350,7 +350,7 @@ if options_.irf
 end
 
 if options_.SpectralDensity.trigger == 1
-    [omega,f] = UnivariateSpectralDensity(oo_.dr,var_list);
+    [oo_] = UnivariateSpectralDensity(M_,oo_,options_,var_list);
 end