%% Lowpass Fourth-Order Section (FOS) Filter
% Design a lowpass fourth-order section (FOS) filter using the |fdesign| 
% function. Using this filter, filter a noisy sinusoidal signal with two 
% tones, one at 3 kHz, and the other at 12 kHz.

%%
% Design a fifth-order filter using the elliptic method in the |'df2tsos'|
% structure. Use L-infinity norm scaling in the frequency domain. Specify the
% passband frequency to be 0.15pi rad/sample and the stopband frequency to be
% 0.25pi rad/sample. Specify 1 dB of allowable passband ripple and a stopband 
% attenuation of 60 dB. 

Fp = 0.15; 
Fst = 0.25;
Ap = 1; 
Ast = 60;

%%
% The filter coefficients are scaled using an |fdopts.sosscaling| object.
% The scaling object is defined to have no numerator constraints, and the
% |ScaleValueConstraint| is set to |'unit'|, specifying the scaling to be 
% unity scaling.

fdo = fdopts.sosscaling;
fdo.NumeratorConstraint='none';
fdo.ScaleValueConstraint='unit';

f = fdesign.lowpass('Fp,Fst,Ap,Ast',Fp,Fst,Ap,Ast);
hFilter = design(f,'ellip','SystemObject',true,...
    'FilterStructure','df2tsos','SOSScaleNorm','Linf',...
    'SOSScaleOpts',fdo);

%%
% Visualize the lowpass frequency response of the designed filter using
% |fvtool|.
fvtool(hFilter)

%%
% Extract the SOS matrix (second-order section representation) of the filter.
sosV =  hFilter.SOSMatrix;

%%
% Extract the numerator and denominator coefficients from the SOS Matrix.
num = zeros(size(sosV,1),5);
den = zeros(size(sosV,1),5);

for i = 1:size(sosV,1)
    [num0,den0] = iirlp2bp(sosV(i,1:3),sosV(i,4:6),Fp,[0.25,0.75]);
    num(i,1:length(num0)) = num0;
    den(i,1:length(num0)) = den0;
end

%%
% Create a fourth-order section filter using the extracted numerator and
% denominator coefficients.
fos = dsp.FourthOrderSectionFilter(num,den);

%%
% Visualize the frequency response of the fourth-order section filter using
% |fvtool|.
fvtool(fos);

%% 
% The input is a sum of two sine waves with frequencies 3 kHz and 12 kHz,
% respectively. The input sample rate is 44.1 kHz, and the frame size is set 
% to 1024 samples.

fs = 44100;
FrameLength = 1024;

SINE1 = dsp.SineWave('SamplesPerFrame',FrameLength,'SampleRate',fs,'Frequency',3000);
SINE2 = dsp.SineWave('SamplesPerFrame',FrameLength,'SampleRate',fs,'Frequency',12000);

%%
% Initialize a spectrum analyzer to visualize the signal spectra.
scope = dsp.SpectrumAnalyzer(...
    'SampleRate',fs,...
    'PlotAsTwoSidedSpectrum',false,...
    'Method','Filter bank',...
    'Title','Original and Filtered Signals',...
    'ShowLegend',true,...
    'YLimits',[-180 50],...
    'ChannelNames',{'Original Signal','Filtered Signal'});



%%
% Filter the noisy input signal with the fourth-order section filter.
% Visualize the spectrum of the original signal and the filtered signal
% using the spectrum analyzer.

for index = 1:1000    
    x = SINE1() + SINE2()+ 0.001*randn(FrameLength,1);    
    y = fos(x);
    scope([x,y]);    
end


%%
% Copyright 2018 The MathWorks, Inc.