%% Bandstop Filtering of Sinusoids
% Construct a bandstop filter to reject the discrete frequency band between 3π/8 and 5π/8 rad/sample. With a sampling frequency of 48 kHz, these values translate to a frequency range of [9 15] kHz. Apply the filter to a discrete-time signal consisting of the superposition of three discrete-time sinusoids.
% The filter is designed by first creating a bandstop filter design specifications object, and then passing the object as an input to the design function.
%
%% Design Bandstop Filter
% Create a bandstop filter design specifications object using fdesign.bandstop.

bandstopSpecs = fdesign.bandstop(1/4,3/8,5/8,6/8,1,60,1)

% List the available design methods for this object.

designmethods(bandstopSpecs)

% To design an equiripple filter, pick 'equiripple'.

bsFilter = design(bandstopSpecs,'equiripple','Systemobject',true)

% Visualize the frequency response of the designed filter.

fvtool(bsFilter,'Fs',48000)

%
%% Create Sinusoidal Signal
% Create a signal that is a sum of three sinusoids with frequencies at 1 kHz, 12 kHz, and 16 kHz. Initialize Spectrum Analyzer to view the original signal and the filtered signal.

Sine1 = dsp.SineWave('Frequency',1e3,'SampleRate',44.1e3,'SamplesPerFrame',4000);
Sine2 = dsp.SineWave('Frequency',12e3,'SampleRate',44.1e3,'SamplesPerFrame',4000);
Sine3 = dsp.SineWave('Frequency',16e3,'SampleRate',44.1e3,'SamplesPerFrame',4000);

SpecAna = dsp.SpectrumAnalyzer('PlotAsTwoSidedSpectrum',false, ...
    'SampleRate',Sine1.SampleRate, ...
    'NumInputPorts',2,...
    'ShowLegend',true, ...
    'YLimits',[-240,45]);

SpecAna.ChannelNames = {'Original noisy signal','Filtered signal'};

%
%% Filter Sinusoidal Signal
% Filter the sinusoidal signal using the bandstop filter that has been designed. View the original signal and the filtered signal in the Spectrum Analyzer. The tone at 1 kHz is unaffected. The tone at 12 kHz is filtered out and attenuated, and the tone at 16 kHz is mildly attenuated because it appears in the transition band of the filter.

for i = 1 : 1000
    x = Sine1()+Sine2()+Sine3();
    y = bsFilter(x);
    SpecAna(x,y);
end
release(SpecAna)

% Copyright 2012 The MathWorks, Inc.