ARL-TR-8272 ● JAN 2018

US Army Research Laboratory

An Automated Energy Detection Algorithm 

Based on Morphological and Statistical 

Processing Techniques

by Kwok F Tom

Approved for public release; distribution is unlimited.
NOTICES

Disclaimers

The findings in this report are not to be construed as an official Department of the 

Army position unless so designated by other authorized documents.

Citation of manufacturer’s or trade names does not constitute an official 

endorsement or approval of the use thereof.

Destroy this report when it is no longer needed. Do not return it to the originator.
ARL-TR-8272 ● JAN 2018

US Army Research Laboratory

An Automated Energy Detection Algorithm 

Based on Morphological and Statistical 

Processing Techniques

by Kwok F Tom

Sensors and Electron Devices Directorate, ARL

Approved for public release; distribution is unlimited.
ii

REPORT DOCUMENTATION PAGE Form Approved

OMB No. 0704-0188

Public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and maintaining the 

data needed, and completing and reviewing the collection information. Send comments regarding this burden estimate or any other aspect of this collection of information, including suggestions for reducing the 

burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and Reports (0704-0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302. 

Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to any penalty for failing to comply with a collection of information if it does not display a currently 

valid OMB control number.

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS.

1. REPORT DATE (DD-MM-YYYY)

January 2018

2. REPORT TYPE

Technical Report

3. DATES COVERED (From - To)

1 October 2016–30 September 2017

4. TITLE AND SUBTITLE

An Automated Energy Detection Algorithm Based on Morphological and 

Statistical Processing Techniques

5a. CONTRACT NUMBER

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

6. AUTHOR(S)

Kwok F Tom

5d. PROJECT NUMBER

5e. TASK NUMBER

5f. WORK UNIT NUMBER

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

US Army Research Laboratory

ATTN: RDRL-SER-E

2800 Powder Mill Road

Adelphi, MD 20783-1138

8. PERFORMING ORGANIZATION REPORT NUMBER

ARL-TR-8272

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSOR/MONITOR'S ACRONYM(S)

11. SPONSOR/MONITOR'S REPORT NUMBER(S)

12. DISTRIBUTION/AVAILABILITY STATEMENT

Approved for public release; distribution is unlimited.

13. SUPPLEMENTARY NOTES

14. ABSTRACT

This report is the result of applying morphological image and statistical processing techniques to the energy detection scenario 

of signals in the RF spectrum domain. The algorithm automatically establishes a detection threshold based on the data of the 

RF spectrum measurement. It is able to determine the threshold level even when there are signals present in the RF spectrum.

15. SUBJECT TERMS

RF spectrum, morphological imaging processing, detection threshold algorithm, erosion, dilation, opening

16. SECURITY CLASSIFICATION OF:

17. LIMITATION

 OF

 ABSTRACT

UU

18. NUMBER

 OF 

 PAGES

70

19a. NAME OF RESPONSIBLE PERSON

Kwok F Tom

a. REPORT

Unclassified

b. ABSTRACT

Unclassified

c. THIS PAGE

Unclassified

19b. TELEPHONE NUMBER (Include area code)

301-394-2612

Standard Form 298 (Rev. 8/98)

Prescribed by ANSI Std. Z39.18
Approved for public release; distribution is unlimited.

iii

Contents

List of Tables v

Preface vi

1. Introduction 1

2. Data Collection and Statistical Summary 1

3. Statistical Processing 2

3.1 Statistical Analysis 2

3.1.1 Moments 3

3.1.2 Mean 3

3.1.3 Variance 3

3.1.4 Standard Deviation 3

3.1.5 Kurtosis 4

3.1.6 Maximum 4

3.1.7 Minimum 4

3.1.8 Median 4

3.1.9 Rank Order Filter 5

3.1.10 Crest Factor (CF) 6

3.1.11 Running Median 6

3.2 Statistical Summary 7

4. Morphological Image Processing 8

5. Algorithm 10

6. Conclusion 11

7. References 12

Appendix A. MATLAB Code 13
Approved for public release; distribution is unlimited.

iv

Appendix B. Graphs of Morphological Processed RF Spectrum Files 25

Appendix C. Graphs of RF Spectrum Files Calculated Detection 

Threshold 43

List of Symbols, Abbreviations, and Acronyms 60

Distribution List 61
Approved for public release; distribution is unlimited.

v

List of Tables

Table 1 Summary of the RF spectrum collection .............................................. 2

Table 2 Statistical analysis of the RF spectrum measurements......................... 8
Approved for public release; distribution is unlimited.

vi

Preface

Energy detection in the RF spectrum is the most basic technique for signal 

detection. Typically, this requires establishing an energy detection threshold based 

on a noise-only condition (i.e., no signal present in the RF spectrum). Initially, the 

area of exploration for this report was to examine the potential of an automatic 

energy detection thresholding algorithm based on the RF measurement. It would 

not require the preliminary RF spectrum noise-only measurements to establish the 

energy detection threshold.

A series of RF spectrum measurements were collected at the US Army Research 

Laboratory in 2013. The data files contained different segments of the RF spectrum 

with various resolution bandwidths. In addition, the data files represented RF 

spectral conditions with and without RF signals. Initially, morphological imaging 

processing techniques were investigated to determine the potential for an automated 

energy detection algorithm. 

The technique of “opening” was very successful in determining an energy detection 

threshold for the RF spectrum with signal and noise-only environments. This led 

to the investigation of other techniques that could be used in automated energy 

detection algorithms. This exploration has resulted in a series of techniques for 

automatic energy detection thresholding determination. This is the first of 5 reports 

that detail the energy detection techniques examined with the recorded RF spectrum 

measurements.1–4

 1 Tom K. An automated energy detection algorithm based on morphological filter processing with a semidisk structure. Adelphi (MD): Army Research Laboratory (US); 2018 Jan. Report No.: ARL-TR-8271. 2

Tom K. An automated energy detection algorithm based on morphological filter processing with a modified 

watershed transform. Adelphi (MD): Army Research Laboratory (US); 2018 Jan. Report No.: ARL-TR-8270. 3 Tom K. An automated energy detection algorithm based on kurtosis-histogram excision. Adelphi (MD): 

Army Research Laboratory (US); 2018 Jan. Report No.: ARL-TR-8269. 4 Tom K. An automated energy detection algorithm based on consecutive mean excision. Adelphi (MD): 

Army Research Laboratory (US); 2018 Jan. Report No.: ARL-TR-8268.
Approved for public release; distribution is unlimited.

1

1. Introduction

Energy detection is the simplest method of detecting a signal in the frequency 

spectrum. When doing so, a comparison is made between the frequency spectral 

component energy and a detection threshold level. Knowledge of the frequency 

spectrum is usually necessary to establish this detection threshold level. Various 

methods can be used to determine the detection threshold level, such as establishing 

the system noise statistics offline and setting the threshold for a given probability 

of detection versus probability of false alarm. 

The goal of this study was to develop an algorithm that can analyze the frequency 

spectrum and establish a threshold value based on the spectral components. The 

automated processing employs techniques from image processing and statistical 

analysis. A combination of techniques from these 2 areas resulted in an algorithm 

for determining a threshold detection level. From image processing, we applied 

morphological filtering to the spectral data to establish a noise level. The 2 

fundamental operations in this threshold-establishing algorithm are called “Erode 

and Dilate”. This series of operations is performed on the spectral data until certain 

statistical criteria are met. An offset was then added to the morphological filtered 

data to form the threshold for energy detection.

2. Data Collection and Statistical Summary

The local RF spectrum was measured in 2013 on the rooftop of building 204 at the 

US Army Research Laboratory’s (ARL’s) Adelphi location. An Agilent N9342CN 

spectrum analyzer and a Discone antenna were used to collect RF spectrum data. 

This spectrum analyzer was operated under the control of a LabVIEW software 

program to acquire and store data with different resolution bandwidth (RBW) from 

1 kHz to 1 MHz. 

These data files represent various sizes of RF spectral coverage from 10 MHz up 

to 4 GHz. The number of data files for each spectral band varied. The larger RBW 

files were acquired over seconds of data acquisition time versus the small RBW 

files. Small RBW provides fine spectral resolution, but it impacts data acquisition 

time for spectral coverage and data size. Depending on the RBW and spectral 

coverage, a data file could require a few hours of acquisition.

Table 1 summarizes the data collection measurement files. The RF spectral bands 

covered the spectrum from 10 MHz to 4 GHz. In general, various spectral bands 

were measured with 4 RBW configuration. Data file size is inversely proportional 
Approved for public release; distribution is unlimited.

2

to the RBW. Data file size is proportional to the spectral band coverage. The data 

size varied from approximately 1 KSample to 4 MSample data points per file. 

Table 1 Summary of the RF spectrum collection

Spectral band coverage

Number of RBW 

measurement RBW

1.1–1.6 GHz 4 1 kHz, 10 kHz, 100 kHz, 

1 MHz

2–3 GHz 4 1 kHz, 10 kHz, 100 kHz, 

1 MHz

3–4 GHz 3 10 kHz, 100 kHz, 1 MHz

4–6 GHz 3 10 kHz, 100 kHz, 1 MHz

10 MHz–1 GHz 4 1 kHz, 10 kHz, 100 kHz, 

1 MHz

10 MHz–2 GHz 4 1 kHz, 10 kHz, 100 kHz, 

1 MHz

10 MHz–3 GHz 4 1 kHz, 10 kHz, 100 kHz, 

1 MHz

10 MHz–4 GHz 4 1 kHz, 10 kHz, 100 kHz, 

1 MHz

100 MHz–1 GHz 1 100 kHz

3. Statistical Processing

3.1 Statistical Analysis

Statistical analysis is the mathematical science dealing with the analysis or 

interpretation of data. The data analyst uses a few straightforward statistical 

techniques as a means of summarizing the collected data. These statistical 

techniques are under the area of descriptive statistics, which is a methodology to 

condense the data in quantitative terms.

In commercial prognostics and diagnostic vibrational monitoring applications, 

statistical techniques that are mainly used for alarm purposes in industrial plants 

are the statistical moments of order 2, 3, and 4. The probability density function 

(PDF) of the vibrational time series of a good bearing has a Gaussian distribution 

(also known as a normal distribution), whereas a damaged bearing results in a nonGaussian distribution with dominant tails because of a relative increase in the 

number of high levels of acceleration. These techniques can be applied to the RF 

spectral data with a different interpretation of the results.1
Approved for public release; distribution is unlimited.

3

3.1.1 Moments

If these moments are calculated about the mean, they are called central statistical 

moments. The first and second moments are well known, being the mean and the 

variance, respectively. These are analogous to the first and second area moments of 

inertia with the area shape defined by the PDF. The third moment is termed 

skewness and the fourth moment is termed kurtosis. The general equation for the 

order of moment is as follows:

MMpp = 1

NN�(xxii − xx�)pp

NN

ii=1

,

wwheeeeee pp iiii tthee oo oo oooo tthee mm mm mm ,

NN iiii tthee nn nn oooo dd dd vv vv ,

ii iiii tthee ii oooo tthee dd dd vv vv ,

xx� iiii tthee mm mm vv vv oooo tthee dd dd .

3.1.2 Mean

Mean is the most common measure of a statistical distribution. In this case, mean 

is the arithmetic average for a set of measurements. 

xx� = μμ = 1

NN�xxii

NN

ii=1

.

3.1.3 Variance

Variance is a measure of the dispersion of a waveform about its mean—also called

the second moment of the measurements.

σσ2 = 1

NN�(xxii − xx�)2

NN

ii=1

.

3.1.4 Standard Deviation

Standard deviation is a measure of the variation of a set of data values. The standard 

deviation is defined as the square root of the variance moment.
Approved for public release; distribution is unlimited.

4

σσ = �1

NN�(xxii − xx�)2

NN

ii=1

.

3.1.5 Kurtosis

Kurtosis is the fourth statistical moment, normalized by the standard deviation to 

the fourth power. It is a measure of whether the data are peaked or flat relative to a 

normal distribution. The noise in the RF spectrum is typically considered to have a 

normal distribution. The normal distribution has a value of 3.

κκ = MM4

σσ4 .

κκ =

1

NN ∑ (xxii − xx�) NN 4

ii=1

σσ4 .

κκ = 1

NNσσ4�(xxii − xx�)4

NN

ii=1

.

3.1.6 Maximum

“Max” is the largest value of a set of numbers.

yy = max[xx( )].

3.1.7 Minimum

“Min” is the smallest value of a set of numbers.

yy = min[xx( )].

3.1.8 Median

The statistical median is an order statistic that gives the “middle” value of a set of 

samples. Median is the middle value of a set of data values that divides the set into 

2 groups. Half the groups exist below and half exist above this value.
Approved for public release; distribution is unlimited.

5

3.1.9 Rank Order Filter 

Rank order filter is a sorting process by which a set of numbers is ordered from the 

smallest to the largest value. Rank order filtering is a nonlinear filtering technique 

that orders the contents of a filter kernel and selects the sample indexed by rank 

from the magnitude ordered samples. 

Rank order filtering can be summarized as follows:

yy = ∑ iixx�(ii)

NN

ii=1 ,

where xx�(ii), = 1, ⋯ , NN is the result of sorting the data in ascending order. With this 

definition, the max, min, and median can be obtained from a rank order filter as 

follows2:

Min:

yy = ∑ iixx�(ii)

NN

ii=1

aaii = 11 ffffff ii = 11

00 oooooo oo oooo

Max:

yy = ∑ iixx�(ii)

NN

ii=1

ii = 1 ii = NN

0 oo heeee ee

Median:

yy = ∑ iixx�(ii)

NN

ii=1

ii = 1 ii = NN

2

0 oo heeee ee
Approved for public release; distribution is unlimited.

6

3.1.10 Crest Factor (CF)

Crest factor (CF) is a measure of a waveform showing the ratio of peak values to 

the effective value. In other words, CF indicates how extreme the peaks are in a 

waveform.

CC FF FF = |xx|pppppppp

xxrrrrrr

.

Noise sources are characterized by their CF, which is the peak to average ratio of 

the noise. In a technical bulletin, XiTRON reported CF values between 5 and 7 for 

random noise.3 For example, a 5:1 CF of the noise voltage is 20log5 = 14 dB. This 

is a measure of the quality of the noise distributions and one way to measure its 

Gaussian nature.4 For the purpose of algorithm development, the CF equation was 

modified as follows:

CC FF FF = MM

MM .

3.1.11 Running Median

To define the running median filter, let {x} be a discrete time series. The running 

median passes a window over the sequence {x} that selects, at each instant m, an 

odd number of samples to comprise the observation vector x(n). The observation 

window is typically symmetric and centered at n, resulting in 

xx( ) = [xx( − NN1), ⋯ , xx( ), ⋯ , xx( + NN1)]

TT ,

where NN1 may range in value over the nonnegative integers, and NN = 2NN1 + 1 is

the (odd valued) window size. While processing such noncausal observation 

vectors has traditionally been referred to as smoothing, we loosen the terminology 

somewhat and refer to the processing of both causal and noncausal observations as 

simply filtering. The median filter operating on the input sequence {xx} produces 

the output sequence{yy}, where at time index ,

yy( ) = MMEEEE[xx( )]

= MM vv vv oooo [xx( − NN1), ⋯ , xx( ), ⋯ , xx( + NN1)].

That is, the samples in the observation window are sorted, and the middle, or 

median, value is taken as the output.5
Approved for public release; distribution is unlimited.

7

3.2 Statistical Summary 

There were 31 different groupings of the RF spectrum data measurements. The 

number of data files under each of the main groupings was varied. For the purpose 

of developing the algorithm, only a single data file was selected from each group. 

Each data file was processed to obtain the following characteristics: RBW, mean, 

standard deviation, median, max, min, kurtosis, and CF. A summary of the results 

is in Table 2. 

The following results were noted:

• The smaller the RBW, the lower the noise floor. The equation for thermal 

noise power is PP = kkkkkk, where B is bandwidth. In this case, the RBW of 1 

kHz has the lowest noise value. The RBW of 1 MHz has the highest noise 

value.

• Each 10-fold increase in bandwidth results in a 10-dB increase in noise 

power. This relationship is illustrated in this data set.

• The calculated mean and median values are very close for a given RF 

measurement configuration.

• The CF for noise-only data files was on the order of approximately 10 to 

13. Noise-only data files were estimated by visually inspecting the spectrum 

plot.

• The RF spectrum is not necessarily flat across the band.
Approved for public release; distribution is unlimited.

8

Table 2 Statistical analysis of the RF spectrum measurements

4. Morphological Image Processing

A description of a technique for automatically estimating the noise floor spectrum 

was given in a conference back in 1997.6 This technique is based on applying the 

morphological binary image processing operators to the RF spectrum. The 

technique works well in both flat and nonflat noise floor spectra. Morphology 

image processing is a set of nonlinear operations. Ready et al. note that “humans 

are good at estimating the noise floor spectrum by ‘eyeballing’ a spectral plot. 

Intuitively, we separate the spectral humps from the noise floor spectrum by 

eliminating those parts of the spectrum shape that are due to signals and visually 

draw in the noise floor spectrum.”6

Morphology is a broad set of processing techniques that process images based on 

shapes. Morphological operations apply a structuring element to an input image, 

creating an output image of the same size. In a morphological operation, the value 

of each pixel in the output image is based on a comparison of the corresponding 

Filename Spectral Band RBW Mean SD Median Max Min Kurtosis CF

Air_test_1.1GHz_1.6GHzb_03_28_14_06_33_22 1.1 - 1.6 GHz 1000 -112.694 3.8444 -112.3 -101.2 -138.4 3.5571 11.1

Air_test_1.1GHz_1.6GHza_03_27_14_07_14_05 1.1 - 1.6 GHz 10000 -102.12 3.7226 -101.8 -91.69 -124.5 3.5134 10.11

Air_test_1.1GHz_1.6GHzc_03_31_14_06_47_26 1.1 - 1.6 GHz 100000 -90.9109 2.597 -90.75 -83.13 -101.2 3.0217 7.62

Air_test_1.1GHz_1.6GHzd_04_01_14_06_54_03 1.1 - 1.6 GHz 1000000 -80.5634 3.5741 -80.745 -55.36 -89.8 14.6861 25.385

Air_test_2GHz_3GHza_06_05_14_07_15_58 2 - 3 GHz 1000 -111.219 4.4991 -111.1 -77.55 -139 6.8918 33.55

Air_test_2GHz_3GHzb_05_29_14_06_28_28 2 - 3 GHz 10000 -101.006 4.1808 -100.8 -74.67 -122 6.1995 26.13

Air_test_2GHz_3GHzc_05_29_14_04_09_27 2 - 3 GHz 100000 -89.5666 3.4154 -89.68 -64.27 -100.3 13.2282 25.41

Air_test_2GHz_3GHzd_05_28_14_00_06_18 2 - 3 GHz 1000000 -79.6014 3.0753 -79.6 -58.74 -88.88 8.4466 20.86

Air_test_3GHz_4GHzb_06_12_14_06_30_17 3 - 4 GHz 10000 -100.51 3.7072 -100.2 -89.68 -123.7 3.5901 10.52

Air_test_3GHz_4GHzC_06_11_14_06_44_40 3 - 4 GHz 100000 -89.1625 2.5513 -89.075 -80.68 -102.3 3.0774 8.395

Air_test_3GHz_4GHzd_06_09_14_07_09_51 3 - 4 GHz 1000000 -79.0828 2.445 -78.94 -72.76 -88.5 2.9651 6.18

Air_test_4GHz_6GHza_04_24_14_07_05_26 4 - 6 GHz 1000 -107.559 4.1793 -107.2 -94.31 -136.3 3.3696 12.89

Air_test_4GHz_6GHzb_04_28_14_06_23_12 4 - 6 GHz 10000 -97.245 4.0398 -96.95 -84.05 -119.7 3.3249 12.9

Air_test_4GHz_6GHzc_04_29_14_06_48_42 4 - 6 GHz 100000 -85.9125 3.082 -85.77 -75.84 -98.78 2.9068 9.93

Air_test_10MHz_1GHza_04_17_14_07_09_30 10 MHz - 1 GHz 1000 -110.348 8.6126 -112.1 -28.37 -140.3 6.4788 83.73

Air_test_10MHz_1GHzb_04_21_14_07_17_44 10 MHz - 1 GHz 10000 -100.168 8.4775 -102 -27.95 -126.7 6.7602 74.05

Air_test_10MHz_1GHzc_04_22_14_06_39_58 10 MHz - 1 GHz 100000 -87.948 8.8113 -90.74 -27.28 -103.6 6.513 63.46

Air_test_10MHz_1GHzd_04_23_14_06_36_16 10 MHz - 1 GHz 1000000 -78.1054 8.3327 -80.59 -27.82 -89.88 7.5282 52.77

Air_test_10MHz_2GHzd_03_06_14_08_52_04 10 MHz - 2 GHz 1000 -111.275 6.9798 -112.1 -27.34 -138.3 9.5085 84.76

Air_test_10MHz_2GHzC_02_27_14_06_43_55 10 MHz - 2 GHz 10000 -100.894 7.0157 -101.8 -26.58 -125.9 10.4224 75.22

Air_test_10MHz_2GHzB_02_25_14_06_35_53 10 MHz - 2 GHz 100000 -89.242 6.7487 -90.59 -26.09 -102.4 12.7961 64.5

Air_test_10MHz_2GHz_02_20_14_06_55_36 10 MHz - 2 GHz 1000000 -78.8448 7.3203 -80.6 -26.54 -92.41 11.3819 54.06

Air_test_10MHz_3GHz_03_11_14_06_32_51 10 MHz - 3 GHz 1000 -111.373 6.1356 -111.8 -29.42 -140.5 11.0474 82.38

Air_test_10MHz_3GHzB_03_13_14_06_32_45 10 MHz - 3 GHz 10000 -100.787 6.3591 -101.4 -25.36 -125.9 11.7254 76.04

Air_test_10MHz_3GHzC_03_18_14_06_32_08 10 MHz - 3 GHz 100000 -89.1614 6.0257 -90.14 -24.16 -102.4 15.3376 65.98

Air_test_10MHz_3GHzd_03_20_14_08_09_07 10 Mhz - 3 GHz 1000000 -78.6492 6.7418 -80.8 -25.81 -90.38 13.0673 54.99

Air_test_10MHz_4GHza_04_10_14_07_01_22 10 MHz - 4 GHz 1000 -111.158 5.6813 -111.4 -28.22 -143.1 10.9573 83.18

Air_test_10MHz_4GHzb_04_11_14_06_10_18 10 MHz - 4 GHz 10000 -100.838 5.5851 -101.1 -26.85 -128 12.2576 74.25

Air_test_10MHz_4GHzc_04_15_14_07_01_09 10 MHz - 4 GHz 100000 -89.4431 5.1076 -90.03 -25.05 -103.5 18.317 64.98

Air_test_10MHz_4GHzd_04_16_14_06_42_27 10 MHz - 4 GHz 1000000 -79.0396 5.5441 -79.89 -26.72 -91.64 16.9073 53.17

Air_test_100MHz_1GHz_02_19_14_07_20_46 100 MHz - 1 GHz 100000 -88.0772 9.0701 -90.82 -23.91 -102.2 6.994 66.91
Approved for public release; distribution is unlimited.

9

pixel in the input image with its neighbors. By choosing the size and shape of the 

neighborhood, you can construct a morphological operation that is sensitive to 

specific shapes in the input image.

The most basic morphological operations are dilation and erosion. Dilation adds 

pixels to the boundaries of the objects in an image, while erosion removes pixels 

on object boundaries. The number of pixels added or removed from the objects in 

an image depends on the size and shape of the structuring element used to process 

the image. In the morphological dilation and erosion operations, the state of any 

given pixel in the output image is determined by applying a rule to the 

corresponding pixel and its neighbors in the input image. These rules are known as

the erosion and dilation:

Erosion:

AA ⊖ = ⋂

∈

AA−bb.

Dilation:

AA⨁ = ⋃

∈ AA .

The opening of A by B is obtained by the dilation of A by B, followed by erosion of 

the resulting structure by B.

Opening:

AA ∘ = (AA ⊖ )⨁ .

The application of the opening technique on the RF spectrum is the basis for noise 

floor estimation. The estimation of the noise floor is used as a reference level to set 

a threshold detection level.7

In signal, statistical, and image processing, the minimum and maximum operators 

are typically encountered. Gil and Kimmel note, “In mathematical morphology, the 

result of such an operator is referred to as the erosion (or dilation) of the signal with 

a structuring element given by a pulse of width p.”8 

For the 1-D case, this reduces to a simple filter of just providing the max or min 

value of a set of values.

8

1-D max filter: Given a sequence xx0, ⋯ , xxnn−1 and an integer pp > 1, compute 

yyii = max

0≤jj<pp

xxii+jj
Approved for public release; distribution is unlimited.

10

for = 0, ⋯ , − pp .

1-D min filter: Given a sequence xx0, ⋯ , xxnn−1 and an integer pp > 1, compute

yyii = min

0≤jj<pp

xxii+jj

For = 0, ⋯ , − pp .

5. Algorithm

The algorithm was originally evaluated and developed in LabVIEW. This 

development process was then executed in MATLAB. Appendix A is the code used 

to process and generate the enclosed RF spectrum signatures with the 

corresponding results of the morphological processing and resultant threshold 

detection level. The following is a description of the detection thresholding level 

generation based on morphological filter processing:

1) Determine the RBW of the spectral data.

2) Determine some statistics on the spectral data file: median, max, and CF.

3) Perform a sliding running median on data for a window size of 5. The start 

and end of the spectral data are somewhat an issue. In this case, the first 5 

points used the median value determined at the first valid value. The last 5 

points used the median value determined at the last valid point.

4) Perform an erosion operation on the median filtered spectral data as 

computed in step 3. Starting from the lowest-frequency component, create

the eroded spectral array with an initial sliding window of size 2. Basically, 

calculate the minimum value for the given window, fill in that position, and 

shift the processing window one frequency position to calculate the next 

value. Continue on across the entire median filtered spectral data array. For 

the last value, use the previous valid value that is calculated.

5) Perform a dilation operation on the eroded spectral array as calculated in 

step 4. In this case, the calculations are performed starting at the highest 

frequency and continue to the lowest frequency. Dilation is simply the 

maximum value in the windowed data set. The initial window size is 2. In 

this case, the maximum value is determined from the sliding window data 

and used to create the dilated spectral array.

6) Calculate statistics on the dilated spectral array from step 5: mean, median, 

max, and CF.
Approved for public release; distribution is unlimited.

11

7) Repeat the morphological filtering of erosion and dilation on the spectral 

array by increasing the window size by 1 for every pass. For example, on 

the second pass, the window size is 3 for erosion and dilation operations. 

The erosion operation is applied to the dilated spectral array from step 5. 

Repeat this process as long as the CF is greater than 10 and the mean power 

greater than 0.01 on each pass of the morphological filtering.

8) If the CF is not greater than 10, then perform the morphological filtering 

operations until the average power is greater than 0.01.

9) Once the convergence has been met, use the resulting morphological 

filtered data array (i.e., the dilated spectral array) to form the threshold as 

follows:

h hoo oo =

MM hoo oooo FF SS AA + 25 −

2.9 x 10 log10 �

RRRRRR

1000�

6. Conclusion

In general, the morphological filtering for establishing the thresholding does well. 

Appendix B displays the graphs of the morphological processing for each of the 

selected RF spectrum data files. The red curve is the result of the morphological 

processing overlaid on top of the RF spectrum signature. 

Appendix C displays the results of the threshold overlaid on top of the RF spectrum 

signature. The graphs are intended to provide the reader with a qualitative sense of 

the effectiveness of the automatic threshold generation algorithm. Actual 

processing time and the number of iterations through the morphological filtering 

routine are highly dependent on the RBW and spectral coverage band. A small 

RBW over a large coverage band requires more time to complete the process. The 

threshold detection level obtained works well for signals present or noise-only RF 

spectrum data files. The red curve is the threshold when added to the morphological 

processed array as shown in Appendix B. 

Overall, a visual inspection of the graphs shows that the algorithm works well. The 

threshold parameters can be modified to obtain the desired level of detection versus 

the false alarm. The algorithm worked well for the RF spectrum where there are 

signals present as well as for the RF spectrum that was judged to be noise only.
Approved for public release; distribution is unlimited.

12

7. References

1. Tom KF. A primer on vibrational ball bearing feature generation for 

prognostics and diagnostics algorithms. Adelphi (MD): Army Research 

Laboratory (US); 2015 Mar. Report No.: ARL-TR-7230.

2. Mitra AK, Lewis TL, Paul AS, Shaw AK. Ultra-wideband radar data models 

and target detection with adaptive rank-order filters. Proc SPIE 4727; 2002.

3. Crest factors and power analyzers. San Diego (CA): XiTRON Technologies; 

nd. Technical Brief TB 103 [accessed 2017 Nov 8]. 

http://www.xitrontech.com/assets/002/5788.pdf.

4. Noise: frequently asked questions. East Hanover (NJ): Noisewave; nd 

[accessed 2017 Nov 8]. http://www.noisewave.com/faq.pdf.

5. Barner KE, Arce GR. Order–statistic filtering and smoothing of time series: 

part II. In: Rao CR, Balakrishnan N, editors. Handbook of statistics. 

Amsterdam (Netherlands): Elsevier Science; 1998. p 555–602; (Order 

statistics and their applications; vol. 16).

6. Ready MJ, Downey ML, Corbalis LJ. Automatic noise floor spectrum 

estimation in the presence of signals. Proceedings of Thirty-First Asilomar 

Conference on Signals, Systems & Computers; 1997 Nov. p 877–881.

7. Sequeira S. Energy based spectrum sensing for enabling dynamic spectrum 

access in cognitive radios [master’s thesis]. [Camden (NJ)]: Rutgers 

University; 2011 May.

8. Gil Y, Kimmel R. Efficient dilation, erosion, opening, and closing algorithms.

IEEE Transactions on Pattern Analysis and Machine Intelligence.

2002;24(12): 1606–1617.
Approved for public release; distribution is unlimited.

13

Appendix A. MATLAB Code
Approved for public release; distribution is unlimited.

14

function MorphologicalFilterROF9report()

%Make selection for data file to process

%SELECT = 1 to 31

SELECT = 15;

dir{1} = 'K:\CognitiveRadar\spectrum monitoring\data\building204-

4c085\Air_test_1.1GHz_1.6GHza\';

filename{1} = 'Air_test_1.1GHz_1.6GHza_03_27_14_07_14_05';

dir{2} = 'K:\CognitiveRadar\spectrum monitoring\data\building204-

4c085\Air_test_1.1GHz_1.6GHzb\';

filename{2} = 'Air_test_1.1GHz_1.6GHzb_03_28_14_06_33_22';

dir{3} = 'K:\CognitiveRadar\spectrum monitoring\data\building204-

4c085\Air_test_1.1GHz_1.6GHzc\';

filename{3} = 'Air_test_1.1GHz_1.6GHzc_03_31_14_06_47_26';

dir{4} = 'K:\CognitiveRadar\spectrum monitoring\data\building204-

4c085\Air_test_1.1GHz_1.6GHzd\';

filename{4} = 'Air_test_1.1GHz_1.6GHzd_04_01_14_06_54_03';

dir{5} = 'K:\CognitiveRadar\spectrum monitoring\data\building204-

4c085\Air_test_2GHz_3GHza\';

filename{5} = 'Air_test_2GHz_3GHza_06_05_14_07_15_58';

dir{6} = 'K:\CognitiveRadar\spectrum monitoring\data\building204-

4c085\Air_test_2GHz_3GHzb\';

filename{6} = 'Air_test_2GHz_3GHzb_05_29_14_06_28_28';

dir{7} = 'K:\CognitiveRadar\spectrum monitoring\data\building204-

4c085\Air_test_2GHz_3GHzc\';

filename{7} = 'Air_test_2GHz_3GHzc_05_29_14_04_09_27';

dir{8} = 'K:\CognitiveRadar\spectrum monitoring\data\building204-

4c085\Air_test_2GHz_3GHzd\';

filename{8} = 'Air_test_2GHz_3GHzd_05_28_14_00_06_18';

dir{9} = 'K:\CognitiveRadar\spectrum monitoring\data\building204-

4c085\Air_test_3GHz_4GHzb\';

filename{9} = 'Air_test_3GHz_4GHzb_06_12_14_06_30_17';

dir{10} = 'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_3GHz_4GHzC\';

filename{10} = 'Air_test_3GHz_4GHzC_06_11_14_06_44_40';

dir{11} = 'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_3GHz_4GHzd\';

filename{11} = 'Air_test_3GHz_4GHzd_06_09_14_07_09_51';
Approved for public release; distribution is unlimited.

15

dir{12} = 'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_4GHz_6GHza\';

filename{12} = 'Air_test_4GHz_6GHza_04_24_14_07_05_26';

dir{13} = 'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_4GHz_6GHzb\';

filename{13} = 'Air_test_4GHz_6GHzb_04_28_14_06_23_12';

dir{14} = 'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_4GHz_6GHzc\';

filename{14} = 'Air_test_4GHz_6GHzc_04_29_14_06_48_42';

dir{15} = 'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_10MHz_1GHza\';

filename{15} = 'Air_test_10MHz_1GHza_04_17_14_07_09_30';

dir{16} = 'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_10MHz_1GHzb\';

filename{16} = 'Air_test_10MHz_1GHzb_04_21_14_07_17_44';

dir{17} = 'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_10MHz_1GHzc\';

filename{17} = 'Air_test_10MHz_1GHzc_04_22_14_06_39_58';

dir{18} = 'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_10MHz_1GHzd\';

filename{18} = 'Air_test_10MHz_1GHzd_04_23_14_06_36_16';

dir{19} = 'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_10MHz_2GHz\';

filename{19} = 'Air_test_10MHz_2GHz_02_20_14_06_55_36';

dir{20} = 'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_10MHz_2GHzB\';

filename{20} = 'Air_test_10MHz_2GHzB_02_25_14_06_35_53';

dir{21} = 'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_10MHz_2GHzC\';

filename{21} = 'Air_test_10MHz_2GHzC_02_27_14_06_43_55';

dir{22} = 'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_10MHz_2GHzd\';

filename{22} = 'Air_test_10MHz_2GHzd_03_06_14_08_52_04';

dir{23} = 'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_10MHz_3GHz\';

filename{23} = 'Air_test_10MHz_3GHz_03_11_14_06_32_51';

dir{24} = 'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_10MHz_3GHzB\';

filename{24} = 'Air_test_10MHz_3GHzB_03_13_14_06_32_45';

dir{25} = 'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_10MHz_3GHzC\';
Approved for public release; distribution is unlimited.

16

filename{25} = 'Air_test_10MHz_3GHzC_03_18_14_06_32_08';

dir{26} = 'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_10MHz_3GHzd\';

filename{26} = 'Air_test_10MHz_3GHzd_03_20_14_08_09_07';

dir{27} = 'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_10MHz_4GHza\';

filename{27} = 'Air_test_10MHz_4GHza_04_10_14_07_01_22';

dir{28} = 'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_10MHz_4GHzb\';

filename{28} = 'Air_test_10MHz_4GHzb_04_11_14_06_10_18';

dir{29} = 'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_10MHz_4GHzc\';

filename{29} = 'Air_test_10MHz_4GHzc_04_15_14_07_01_09';

dir{30} = 'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_10MHz_4GHzd\';

filename{30} = 'Air_test_10MHz_4GHzd_04_16_14_06_42_27';

dir{31} = 'K:\CognitiveRadar\spectrum 

monitoring\data\building204-4c085\Air_test_100MHz_1GHz\';

filename{31} = 'Air_test_100MHz_1GHz_02_19_14_07_20_46';

str_meta = 

sprintf('%s%s.mspecrawdata',dir{SELECT},filename{SELECT});

str_data = 

sprintf('%s%s.specrawdata',dir{SELECT},filename{SELECT});

% Get metadata on selected data file

fid_meta = fopen(str_meta);

META = textscan(fid_meta,'%s');

ave = META{1}{1};

ref = META{1}{2};

attn = META{1}{3};

rbw = META{1}{4};

% Read in data file

fid_data = fopen(str_data);

g=0;

f=0;

a=0;

while(g==0)

 ft=fgetl(fid_data);

 f=[f,str2num(ft)];

 at=fgetl(fid_data);

 a=[a,str2num(at)];

 g=feof(fid_data);

end
Approved for public release; distribution is unlimited.

17

f=f(2:end);

a=a(2:end);

% Determine resolution bandwidth

Bandwidth = f(3) - f(2);

% Plot data file

figure(1);

plot(f/1e6,a);

axis tight;

grid;

xlabel('Frequency (MHz)');

ylabel('Power (dBm)');

title(strcat(rbw,'--',ave,'--',ref,'--

',attn),'Interpreter','none');

% Determine some statistics values for data

M = mean(a);

Med = median(a);

S = std(a);

Max = max(a);

Min = min(a);

Kurt = kurtosis(a);

Range = abs(Min - Max);

Number = floor(Range);

Bins = Number * 2;

disp([M S Med Max Min Kurt]);

% Initialize parameters

 

InputIndex = length(a);

ErosionIndex = 1;

InputArray = a;

ErosionArray = a;

DilationArray = a;

Medsize = 5;

Medcount = fix(Medsize);
Approved for public release; distribution is unlimited.

18

Medstart = Medcount + 1;

% Process data file with a running median

for m = Medstart:InputIndex-Medcount

 InputArray(m) = median(a(m-Medcount:m+Medcount)); 

end

for m = 1:Medcount

 InputArray(m) = InputArray(Medstart); 

end

for m = InputIndex-Medcount:InputIndex

 InputArray(m) = InputArray(InputIndex-Medcount);

end

% Plot processed median of data file

figure(2);

plot(InputArray);

PreviousPower = sum(a)/InputIndex;

display (PreviousPower);

% Initialize parameters

j= 1;

DiffPower = 1;

DilationMedian = Med;

DilationMax = Max;

DilationCF = DilationMax-DilationMedian;

 

if DilationCF > 10 && j == 1

 display('CF & Power'); 

% Perform Erosion operation on data file

 while DilationCF > 10

 

 for k = 1:InputIndex-j

 ErosionArray(k) = min(InputArray(k:k+j));

 end

 for k = InputIndex-j:InputIndex
Approved for public release; distribution is unlimited.

19

 

 ErosionArray(k) = InputArray(InputIndex-j);

 

 end

 

% Plot results of Erosion operation 

 figure(5);

 plot(ErosionArray);

% Perform Dilation operation on data file

 DilationArray = ErosionArray;

 for k = InputIndex:j+1

 DilationArray(k) = max(ErosionArray(k-j:k));

 end

 

 for k = j:1

 DilationArray(k) = ErosionArray(j+1);

 end

% Plot results of Dilation operation

 InputArray = DilationArray;

 figure(6);

 plot(DilationArray); 

 

 DilationPower = sum(DilationArray)/InputIndex;

 DilationMedian = median(DilationArray);

 DilationMax = max(DilationArray);

 DiffPower = abs(PreviousPower - DilationPower);

 DilationCF = DilationMax-DilationMedian;

 display ([j DilationPower DilationMedian DilationMax 

DilationCF DiffPower]);

 PreviousPower = DilationPower;

 j= j+1;

 end

 % Check performance of the Erosion & Dilation to see if 

operations 

 % resulted have converged to accept criteria (If power 

difference in 
Approved for public release; distribution is unlimited.

20

 % the RF spectrum is less than 0.1 between iterations previous 

and 

 % current processing of RF spectrum

 

 while DiffPower > 0.1 

 

 display('Power');

 

 for k = 1:InputIndex-j

 ErosionArray(k) = min(InputArray(k:k+j));

 end

 for k = InputIndex-j:InputIndex

 ErosionArray(k) = min(InputArray(k:InputIndex));

 end

 figure(5);

 plot(ErosionArray);

 DilationArray = ErosionArray;

 for k = InputIndex:j+1

 DilationArray(k) = max(ErosionArray(k-j:k));

 end

 

 for k = j:1

 DilationArray(k) = max(ErosionArray(k:1));

 end

 

 InputArray = DilationArray;

 figure(6);

 plot(DilationArray); 

 

 DilationPower = sum(DilationArray)/InputIndex;

 DilationMedian = median(DilationArray);

 DilationMax = max(DilationArray);

 DiffPower = abs(PreviousPower - DilationPower);

 DilationCF = DilationMax-DilationMedian;
Approved for public release; distribution is unlimited.

21

 display ([j DilationPower DilationMedian DilationMax 

DilationCF DiffPower]);

 PreviousPower = DilationPower;

 j= j+1;

 

 end

else

 

 while DiffPower > 0.01

 

 display('Power');

 

 for k = 1:InputIndex-j

 ErosionArray(k) = min(InputArray(k:k+j));

 end

 for k = InputIndex-j:InputIndex

 ErosionArray(k) = min(InputArray(k:InputIndex));

 end

 figure(5);

 plot(ErosionArray);

 DilationArray = ErosionArray;

 for k = InputIndex:j+1

 DilationArray(k) = max(ErosionArray(k-j:k));

 end

 

 for k = j:1

 DilationArray(k) = max(ErosionArray(k:1));

 end

 

 InputArray = DilationArray;

 figure(6);

 plot(DilationArray); 

 

 DilationPower = sum(DilationArray)/InputIndex;

 DilationMedian = median(DilationArray);

 DilationMax = max(DilationArray);
Approved for public release; distribution is unlimited.

22

 DiffPower = abs(PreviousPower - DilationPower);

 DilationCF = DilationMax-DilationMedian;

 display ([j DilationPower DilationMedian DilationMax 

DilationCF DiffPower]);

 PreviousPower = DilationPower;

 j= j+1;

 

 end

 

end

% Plot original RF spectrum and threshold

figure(3);

plot(f/1e6,a);

axis tight;

grid;

xlabel('Frequency (MHz)');

ylabel('Power (dBm)');

title(strcat(filename(SELECT),' -- ',rbw,' 

Hz'),'Interpreter','none');

hold on

% Determine threshold array

Threshold = DilationArray+25-(2.9*log10(Bandwidth));

plot(f/1e6,Threshold,'r');

hold off;

% Plot original RF spectrum and final results of Erosion & 

Dilation

% processing

figure(10);

plot(f/1e6,a);

axis tight;

grid;

xlabel('Frequency (MHz)');

ylabel('Power (dBm)');
Approved for public release; distribution is unlimited.

23

title(strcat(filename(SELECT),' -- ',rbw,' 

Hz'),'Interpreter','none');

hold on

plot(f/1e6,DilationArray,'r');

hold off;

% Save results to data file

lab = num2str(SELECT);

lab1 = strcat('ROF2-Figure ',lab,'.jpg');

saveas(gcf, lab1);

end
Approved for public release; distribution is unlimited.

24

INTENTIONALLY LEFT BLANK.
Approved for public release; distribution is unlimited.

25

Appendix B. Graphs of Morphological Processed RF Spectrum 

Files
Approved for public release; distribution is unlimited.

26
Approved for public release; distribution is unlimited.

27
Approved for public release; distribution is unlimited.

28
Approved for public release; distribution is unlimited.

29
Approved for public release; distribution is unlimited.

30
Approved for public release; distribution is unlimited.

31
Approved for public release; distribution is unlimited.

32
Approved for public release; distribution is unlimited.

33
Approved for public release; distribution is unlimited.

34
Approved for public release; distribution is unlimited.

35
Approved for public release; distribution is unlimited.

36
Approved for public release; distribution is unlimited.

37
Approved for public release; distribution is unlimited.

38
Approved for public release; distribution is unlimited.

39
Approved for public release; distribution is unlimited.

40
Approved for public release; distribution is unlimited.

41
Approved for public release; distribution is unlimited.

42

INTENTIONALLY LEFT BLANK.
Approved for public release; distribution is unlimited.

43

Appendix C. Graphs of RF Spectrum Files Calculated Detection 

Threshold
Approved for public release; distribution is unlimited.

44
Approved for public release; distribution is unlimited.

45
Approved for public release; distribution is unlimited.

46
Approved for public release; distribution is unlimited.

47
Approved for public release; distribution is unlimited.

48
Approved for public release; distribution is unlimited.

49
Approved for public release; distribution is unlimited.

50
Approved for public release; distribution is unlimited.

51
Approved for public release; distribution is unlimited.

52
Approved for public release; distribution is unlimited.

53
Approved for public release; distribution is unlimited.

54
Approved for public release; distribution is unlimited.

55
Approved for public release; distribution is unlimited.

56
Approved for public release; distribution is unlimited.

57
Approved for public release; distribution is unlimited.

58
Approved for public release; distribution is unlimited.

59
Approved for public release; distribution is unlimited.

60

List of Symbols, Abbreviations, and Acronyms

ARL US Army Research Laboratory

CF crest factor

GHz gigahertz

kHz kilohertz

MHz megahertz

PDF probability density function

RBW resolution bandwidth

RF radio frequency
Approved for public release; distribution is unlimited.

61

1 DEFENSE TECHNICAL

(PDF) INFORMATION CTR

DTIC OCA

2 DIR ARL

(PDF) IMAL HRA 

RECORDS MGMT

RDRL DCL

TECH LIB

1 GOVT PRINTG OFC

(PDF) A MALHOTRA

7 ARL

(PDF) RDRL SER E

M CONN

R DEL ROSARIO

K F TOM

D WASHINGTON

RDRL SER M

J SILVIOUS

RDRL SER U

A MARTONE

RDRL SER W

K RANNEY
Approved for public release; distribution is unlimited.

62

INTENTIONALLY LEFT BLANK.