Accelerating the pace of engineering and science

# Documentation Center

• Trial Software

## Creating New Types of System objects: a Teager-Kaiser DESA-2 Operator

This example shows how to design a custom System object™ to estimate the instantaneous frequency of an FM-modulated tone using a specific type of Discrete Energy Separation Algorithm known as DESA-2, based on the Teager-Kaiser Energy Operator

Introduce the Discrete Energy Separation Algorithms (DESA)

The Discrete Energy Separation Algorithms (DESA) provide instantaneous estimates of amplitude and frequency for sinusoidal signals. The basic building block of DESA is the nonlinear Teager-Kaiser Energy Operator (TKEO)

In particular, the DESA-2 algorithm provides instantaneous estimates for amplitude (or AM function) and frequency (or FM function) through the following expressions

Develop with an Object Oriented approach

This example takes the perspective of a MATLAB developer willing to author a DESA-2 operator. It illustrates how to use an Object-Oriented Programming (OOP) approach to designing both the core algorithm and a simple test bench around it. More specifically, all discrete-time algorithmic components are designed as System objects™. System objects simplify the modeling of discrete dynamic systems by providing a coherent programming interface for time-domain simulations. The step method of all System objects is responsible for executing the core algorithm, by producing the required output (or other effect, e.g. some visualization or data storage) based on the input signals provided. As all objects, once created System objects store their own properties and internal states. This allows to efficiently separate initialization tasks (typically executed only once) from the core runtime algorithm. The latter is indeed called iteratively within a simulation loop so it should be as lightweight as possible in the interest of efficiency.

Introduce the simulation setup - test signal and estimation operator

This simulation uses a DESA-2 operator to estimate the instantaneous frequency of an FM-modulated sinusoid with constant amplitude and randomly varying frequency. From a high-level perspective, this is realized by constructing an object to generate a suitable test signal, followed by another object in charge of the actual frequency estimation.

Define system parameters

% Discrete time parameters for the simulation
fs = 2000;
totalTime = 4;
frameLength = 64;
nFrames = floor(totalTime/(frameLength/fs));

% Parameters describing the constant-amplitude FM-modulated sinusoid
amplitude = 2;
frequencyOffset = 200;
frequencyStandardDeviation = 20;
frequencyResolution = 0.1;
frequencyNoiseCutoffFreq = 10;
frequencyNoiseBandRatio = frequencyNoiseCutoffFreq/fs;
frequencyNoiseSeed = 1;


Generate test signal

An instance of a RandomFMToneGenerator is used to generate the test signal. After creating it and setting its properties, its step method is called in a loop to generate consecutive frames of the test signal. dsp.TimeScope is used to visualize the signal dynamically as it is generated. Visually, one can see the period of the sinusoidal signal varying randomly over time.

hFMTone = dspdemo.RandomFMToneGenerator(...
'Amplitude', amplitude,...
'FrequencyOffset', frequencyOffset,...
'FrequencyResolution', frequencyResolution,...
'FrequencyStandardDeviation', frequencyStandardDeviation,...
'FrequencyRandomSeed', 1,...
'FrequencyVariationNormalizedBandwidth',frequencyNoiseBandRatio,...
'SampleRate', fs,...
'SamplesPerFrame', frameLength);

hScope = dsp.TimeScope(...
'YLimits', 2 * [-1.2 1.2],...
'YLabel','Oscillator output',...
'TimeSpan',totalTime/20,...
'BufferLength', 2 * totalTime/20 * fs,...
'FrameBasedProcessing', true,...
'SampleRate', fs,...
'Title', 'Frequency-modulated signal with constant amplitude',...
'ShowLegend', true,...
'ShowGrid', true);

for k = 1:nFrames/2
% FM test tone generation
[testTone, ~] = step(hFMTone);

% Visualization
step(hScope, testTone)
end


Inspect architecture of RandomFMToneGenerator

RandomFMTonegenerator is a custom, MATLAB-authored System object that in turn uses a number of other classes. The latter are a mixture of built-in and custom authored System objects. The following is a UML diagram that summarizes the architecture of RandomFMToneGenerator:

Inspect code of RandomFMToneGenerator

Inspecting the code of RandomFMToneGenerator helps to clarify the roles of the lower-level components. The stepImpl method implements the core functionality of the public step() method. Notice how the separation of reponsibilities helps to simplify the code within stepImpl, which in turn calls the step() methods of the two inner System objects respectively

type(fullfile(matlabroot,'toolbox','dsp','dspdemos', '+dspdemo',...
'RandomFMToneGenerator.m'))

classdef RandomFMToneGenerator < matlab.System
%FMToneGenerator Generates a randomly FM-modulated tone
%
% Copyright 2013 The MathWorks, Inc.

properties
Amplitude = 1
FrequencyOffset = 200
FrequencyResolution = 0.1
FrequencyStandardDeviation = 20
FrequencyRandomSeed = 1
FrequencyVariationNormalizedBandwidth = 0.0025
SampleRate = 2000
SamplesPerFrame = 1
end

properties(Access = private)
pNoiseGenerator
pNCO
pFrequencyGain
end

methods
function obj = RandomFMToneGenerator(varargin)
% Support name-value pair arguments
setProperties(obj,nargin,varargin{:});
end
end

methods (Access=protected)
function setupImpl(obj)
% Initialize noise generation
obj.pNoiseGenerator = dspdemo.BandlimitedNoiseGenerator(...
'NoiseStandardDeviation',obj.FrequencyStandardDeviation,...
'NoiseNormalizedBandwidth',...
obj.FrequencyVariationNormalizedBandwidth,...
'NoiseRandomSeed', obj.FrequencyRandomSeed,...
'SamplesPerFrame', obj.SamplesPerFrame);

% Initialize NCO
fs = obj.SampleRate;
nNCO = ceil(log2(fs/obj.FrequencyResolution));
obj.pNCO = dsp.NCO('SamplesPerFrame', 1, ...
'CustomAccumulatorDataType', numerictype([],nNCO),...
'OutputDataType', 'double');

obj.pFrequencyGain = 2^nNCO/fs;

end

function resetImpl(obj)
reset(obj.pNoiseGenerator)
reset(obj.pNCO)
end

function [fmtone, frequencyActual] = stepImpl(obj)
frequencyActual = step(obj.pNoiseGenerator) + ...
obj.FrequencyOffset;
L = obj.SamplesPerFrame;
fgain = obj.pFrequencyGain;
fmtone = zeros(L, 1);
for k = 1:obj.SamplesPerFrame
fmtone(k) = obj.Amplitude * step(obj.pNCO, ...
int32(frequencyActual(k)*fgain));
end
end

function N = getNumInputsImpl(~)
% Specify number of System inputs
N = 0; % Because stepImpl has no input arguments beyond obj
end

function N = getNumOutputsImpl(~)
% Specify number of System outputs
N = 2; % Because stepImpl has one output
end
end
end



Test BandlimitedNoiseGenerator, building block of RandomFMToneGenerator

A key component of RandomFMToneGenerator is BandlimitedNoiseGenerator. BandlimitedNoiseGenerator generates a zero-mean random signal with a prescribed bandwidth and RMS amplitude. The signal generated by BandlinitedNoiseGenerator defines the frequency variations of the test signal around its offset value. BandlimitedNoiseGenerator can be tested in isolation by creating an istance, setting its relevant parameters and calling its step method in a loop. The same instance of dsp.TimeScope (hScope) is used to visualize the resulting signal over time.

hNoiseGen = dspdemo.BandlimitedNoiseGenerator(...
'NoiseStandardDeviation', frequencyStandardDeviation,...
'NoiseNormalizedBandwidth', frequencyNoiseBandRatio,...
'NoiseRandomSeed', frequencyNoiseSeed,...
'SamplesPerFrame', frameLength);

% Change visualization settings
release(hScope)
hScope.YLimits = 6*frequencyStandardDeviation*[-1 1];
hScope.TimeSpan = totalTime/2;
hScope.BufferLength = totalTime/2 * fs;
hScope.Title = 'Frequency-modulating signal';
hScope.YLabel = 'Frequency variation (Hz)';

% Simulate for totalTime seconds
for k = 1:nFrames
% Random frequency variation signal
frequencyRandZeroMean = step(hNoiseGen);

% Visualization
step(hScope, frequencyRandZeroMean)
end


Inspect architecture and code for the DESA-2 operator

The DESA-2 frequency-estimation operator shown here below uses the same System object interface encountered before. This time though its step method accepts an array of signal samples as input and returns local estimates of tone amplitude and frequency. Comparing the equations for the operator on one hand and the class diagram on the other, notice how the DESA-2 operator is composed of two Teager-Kaiser energy operators of type TeagerKeiserEnergyOperator as expected.

type(fullfile(matlabroot,'toolbox','dsp','dspdemos', '+dspdemo',...
'DesaTwo.m'))

classdef DesaTwo < matlab.System
% DESATWO Discrete Energy Separation Algorithm (DESA-2)
%
% Copyright 2013 The MathWorks, Inc.

properties(Access = public)
SampleRate = 1
end

properties(SetAccess = private, GetAccess = private)
% Two Teager-Kaiser enery operators (with internal states)
pDiscreteEnergyOperatorOne  % TeagerKaiserEnergyOperator
pDiscreteEnergyOperatorTwo  % TeagerKaiserEnergyOperator
% Two delay elements
pDelayOne                   % dsp.Delay
pDelayTwo                   % dsp.Delay
end

methods
% Constuctor
function obj = DesaTwo(varargin)
% Support name-value pair arguments
setProperties(obj,nargin,varargin{:});
end
end

% System object "services" (all protected)
methods(Access=protected)
function [amp, freq] = stepImpl(obj, inSignal)
% Used to avoid division by zero and asin argument > 1
DELTA = 1e-5;
MAX_ENRATIO = 4;

% Use delay element to form input to two Teager Keiser
% Energy Operators
tkeoOneInput = step(obj.pDelayOne, inSignal);
tkeoTwoInput = inSignal - step(obj.pDelayTwo, inSignal);
% Get energy outputs from two energy operator instances
energyOne = step(obj.pDiscreteEnergyOperatorOne,tkeoOneInput);
energyTwo = step(obj.pDiscreteEnergyOperatorTwo,tkeoTwoInput);
% Constrain energy estimates to be positive and non-zero
energyOne = max(DELTA, energyOne);
energyTwo = max(DELTA, energyTwo);
% Constrain energy ratio to be less than 4 - this makes
% argument of asin less than 1
energyRatio = min(energyTwo ./ energyOne, MAX_ENRATIO);
% Compute amplitude and frequency estimates
amp = 2*energyOne./(sqrt(energyTwo));
freq = obj.SampleRate * 1/(2*pi) * asin(1/2 * sqrt(energyRatio) );
end

function numIn = getNumInputsImpl(~)
numIn = 1;
end

function numOut = getNumOutputsImpl(~)
numOut = 2;
end

function setupImpl(obj, inSignal) %#ok<INUSD>
% Instantiate two internal Teager Kaiser Energy Operators
% with Kaiser order = 1 (default)
obj.pDiscreteEnergyOperatorOne = ...
dspdemo.TeagerKaiserEnergyOperator;
obj.pDiscreteEnergyOperatorTwo = ...
dspdemo.TeagerKaiserEnergyOperator;
% Delay #1 has only one tap, Delay #2 has two
obj.pDelayOne = dsp.Delay(1);
obj.pDelayTwo = dsp.Delay(2);

end

function resetImpl(obj)
% Reset all system object properties
reset(obj.pDiscreteEnergyOperatorOne);
reset(obj.pDiscreteEnergyOperatorTwo);
reset(obj.pDelayOne);
reset(obj.pDelayTwo);
end

function releaseImpl(obj)
% Release all system object properties
release(obj.pDiscreteEnergyOperatorOne);
release(obj.pDiscreteEnergyOperatorTwo);
release(obj.pDelayOne);
release(obj.pDelayTwo);
end

end

end


Verify the DESA-2 operator on the test FM signal

DesaTwo can be tested with a signal generated with RandomFMToneGenerator. Beyond the actual FM tone, RandomFMToneGenerator can also return the time-varying frequency values used to generate it. The latter are expected to match the estimate returned by DesaTwo within a certain tolerance. The code here below simulates the two components together and uses again dsp.TimeScope for visual verification.

reset(hFMTone)

hDESA2 = dspdemo.DesaTwo('SampleRate', fs);

% Change visualization settings
release(hScope)
hScope.Title = 'Actual (1) and estimated (2) frequencies';
hScope.YLabel = 'Frequency (Hz)';
hScope.YLimits = 6*frequencyStandardDeviation*[-1 1] + frequencyOffset;
hScope.PlotType = 'Stairs';
hScope.BufferLength = 2 * totalTime/2 * fs;

for k = 1:nFrames
% FM test tone generation
[testTone, frequencyActual] = step(hFMTone);

% Instantaneous frequency estimation (the instantaneous amplitude is
% also estimated but not plotted for simplicity)
[estimatedAmplitude, estimatedNormFreq] = step(hDESA2, testTone);

% Visualization
step(hScope, [frequencyActual, estimatedNormFreq])
end


Reference

P. Maragos, J.F. Kaiser, T.F. Quartieri, Energy Separation in Signal Modulations with Application to Speech Analysis, IEEE Transactions on Signal Processing, vol. 41, No. 10, October 1993