Patent Application
20150169498
1. Field
The present disclosure relates to digital signal processing. More particularly, this invention is directed toward a fast convolution of a signal with a one-sided exponential function.
2. Description of Related Technology
Many fields of engineering and science utilize convolution of a signal of interest with a kernel comprising a class of symmetric functions. Such a class may comprise, e.g., finite window functions such as cosine window, raised cosine window, Lanczos window with small radius, as well as infinite functions like Gaussian function, Gabor atoms (or functions), specifically designed low pass filters, and other symmetric functions known to a person of ordinary skills in the art.
A convolution of a signal with a kernel comprising Gaussian function finds use in such fields as speech recognition, signal equalization, image processing and other fields known to a person of ordinary skills in the art. By means of an example, an edge detection in a signal comprising an image takes an advantage of the fact that the edges of the objects in the image with a characteristic size or spacing related to the widths of the Gaussian functions appear as zero values. By means of another example, Gabor functions, are frequently used as kernels for time-frequency analysis, e.g., feature extraction, especially in texture-based image analysis (e.g. classification, segmentation or edge detection), iris recognition. By means of yet another example, raised cosine function is frequently used as a kernel in matched filters.
A traditional approach to calculating a convolution is to implement the convolution through direct numerical computation between the signal of interest and the kernel, corresponding to a function from, e.g., the above-discussed class of functions. However, in applications, where computational time is critical, such implementation is not suitable and is desired to be improved.
One approach is disclosed in a document by Seeman, M., Zem{hacek over (e)}ik, P., entitled Histogram Smoothing for Bilateral Filter, Proc. GraVisMa, Plzen, Czech Rep., pp. 145-148, 2009. As disclosed therein, the approach considers a kernel comprising Gaussian function. One half of the Gaussian function is approximated by one or more exponential functions, thus decreasing the computational complexity form (m×n) (where n is the number of signal samples and m is the size of the convolution kernel) to n (where n is the number of the signal samples). The convolution with the exponential function is then calculated in a forward direction, i.e., for signal samples starting with the first sample and then in a reverse direction, i.e., for signal samples starting with the last sample. The result of the forward and reverse calculations are then added together.
A problem with this approach is that due to the need for calculation in a forward and reverse directions, the entire signal must be known in advance; therefore, only post-processing of a signal that has been acquired is feasible.
However, many engineering and scientific processes produce signals, which are not necessary finite, e.g., speech, motion pictures, and other signals known to a person of ordinary skills in the art. Accordingly, there is a need in the art for a fast convolution of a signal with an exponential function, used, e.g., to approximate kernels comprising a class of symmetric functions, thus providing a solution to the above identified problem, enabling calculation of convolution of such kernels with both definite and indefinite signals, as well as providing additional advantages.
In an aspect of the disclosure, an apparatus and a method for fast convolution of a signal with a one-sided exponential function according to appended independent claims is disclosed. Additional aspects are disclosed in the dependent claims.
The foregoing aspects described herein will become more readily apparent by reference to the following description when taken in conjunction with the accompanying drawings wherein:
Unless otherwise defined, all terms (including technical and scientific terms) used herein have the same meaning as commonly understood by a person having ordinary skill in the art to which this invention belongs. It will be further understood that terms, such as those defined in commonly used dictionaries, should be interpreted as having a meaning that is consistent with their meaning in the context of the relevant art and this disclosure.
As used herein, the singular forms “a,” an and the are intended to include the plural forms as well, unless the context clearly indicates otherwise. It will be further understood that the terms “comprise,” “comprises,” and/or “comprising,” when used in this specification, specify the presence of stated features, integers, steps, operations, elements, and/or components, but do not preclude the presence or addition of one or more other features, integers, steps, operations, elements, components, and/or groups thereof. The term “and/or” includes any and all combinations of one or more of the associated listed items.
Various disclosed aspects may be illustrated with reference to one or more exemplary configurations. As used herein, the term “exemplary” means “serving as an example, instance, or illustration,” and should not necessarily be construed as preferred or advantageous over other configurations disclosed herein.
Various aspects of the present invention will be described herein with reference to drawings that are schematic illustrations of conceptual configurations of the present invention, unless explicitly noted. The various aspects of this disclosure are provided to enable a person having ordinary skill in the art to practice the present invention. Modifications to various aspects of a presented throughout this disclosure will be readily apparent to a person having ordinary skill in the art, and the concepts disclosed herein may be extended to other applications.
One-dimensional symmetric functions that can be approximated by an exponential function are denoted as w(t), where the t is an independent variable. By means of an example, the well-known one-dimensional Gaussian function G(t) is defined as:
where σ is standard deviation. As well known the equation evaluated over the independent variable can be depicted as the well known symmetrical, bell-shaped curve. By means of another example(s) a raised cosine window C(t) is defined as:
where R is radius. By means of yet another example, Lanczos window is used for signal re-sampling with a radius equal to two signal samples. Lanczos window L(t) is defined as
where S is distance between the samples of the signal to which the original signal is being re-sampled.
As understood by a person of ordinary skills in the art, not every function has an analytical expression. By means of an example, a process may possess a functional relationship between an independent variable and the observed variable. As long as the relationship can be approximated by a sum of one or more exponential functions, which possesses an analytical expression, the function is within the scope of this disclosure.
One half of any of the symmetrical functions w(t) to be used as a kernel, can be approximated as a sum of one or more exponential functions, wherein the approximation by one exponential curve is defined as:
w(t)∝(1−α)−n Eq. (4)
where ∝ indicates approximation; α is a constant determining approximation via the decay constant (1−α)<1; and n is an exponent. The constant α may be determined by any approximation method known to a person of ordinary skills in the art, e.g., the least squares method.
Thus, referring now to
A person of ordinary skills in the art will understand that a single exponential is used for clarity of explanation without a loss of any generality, and extension to a plurality of exponential functions, should the approximation error be unacceptable, is an application of a principle of superposition. Thus by means of an example,
A convolution of a first half of the symmetrical function w(t) approximated by Eq. (4) with a signal of interest ƒ(m) yields an Exponential Moving Average (EMA):
E
α(n,ƒ)=Σm=−∞(1−α)n−m·ƒ(m) Eq. (5)
Instead of calculating the EMA directly from Eq. (5), a recursive relationship may be used:
E
α
n,ƒ=ƒ(n)+(1−α)·Eα(n−1,ƒ) Eq. (6)
Referring to
In block 302, a number of input signal ƒ(m) samples (404) are accumulated into a memory structure, e.g., a register (402). A person of ordinary skills in the art will understand that since the recursive relationship of Eq. (6) operates on only one input signal sample, only a register position (402_1) must contain the input signal sample before the recursion may start. Consequently, the size of the register (402) must comprise at least one position (402_x). The process continues in block 304.
In block 304, the recursive calculation of a sample of the output signal is carried out. The first sample of the signal in the register position (402_1) is provided to a first input of a block (404). A previously calculated sample of the output signal in a position (406_0) of an output register (406) is provided to a block (408) where the sample is multiplied by (1−α) and the result of the multiplication is provided to a second input of the block (404). In the block (404), the two input signals are summed together and provided as an output sample (406_1) into the register (406). A person of ordinary skills in the art will understand that since the recursive relationship of Eq. (6) operates on a previously calculated output signal sample in register position (406_0) and outputs to a different register position (406_1), the register (406) must comprise at least two positions (406_x). The process continues in block 306.
In block 306, the signal samples in both the input register (402) and the output register (406) are shifted one position to the right, so that the register position (402_1) holds the second sample of the input signal, the register position (406_0) holds the previously calculated sample of the output signal, and the register position (406_1) is empty. The process continues in block 304.
To calculate the other half of the symmetrical function w(t) approximated by Eq. (4) with the signal of interest ƒ(m), the function (1−α)−n used in Eq. (5) must be mirrored, i.e., a function (1−α)+n must be used. Considering Eq. (4) for such mirrored function, a recursive relationship may be used:
Examining Eq. (7) two observations can be made. The first one is that the relationship is non-causal, since the output, i.e., the Eα(n, ƒ) depends only on one sample of the signal ƒ(m) and one Eα(n+1, ƒ) in the future. Furthermore, the initial condition for Eq. (7) is determined as:
E
α(N,ƒ)=Em=0N(1−α)N−m·ƒ(n) Eq. (8)
where N is the last sample of the signal ƒ(m). However, since as discussed supra, the signal ƒ(m) has, at least in theory, infinite length practical realization of Eq. (8) is impossible.
The second observation concerns the coefficient
a power of which grows with a number of samples of the signal ƒ(m) up to the power of N. The number of significant digits in a result of a multiplication is a sum of the number of significant digits representing the factors in the multiplication. Since a calculation of a power of a number is implemented as multiplication, it follows that the required precision, i.e., the size of an accumulator holding the significant digits of the multiplication grows with the number of samples of the signal N. However, it is clear from Eq. (7) that all the significant digits must be retained. Thus for an infinite signal length, the precision would have to be infinite, which means that practical realization of Eq. (7) would be impossible in any hardware.
To resolve these issues, the mirrored function (1−α)+n is truncated to a pre-determined number of samples W. Consequently, the initial condition given by the Eq. (8), can now be evaluated, because the contribution of the signal ƒ(m) to the initial condition is calculated over a limited number of samples N=W. However, the truncation causes an approximation error Ek in the kernel summation, given as:
Because the approximation error Ek decreases with increasing number W, the pre-determined number W may be determined in accordance with an approximation error Ek acceptable for a particular application.
Once the pre-determined number W is determined, the required number of the significant digits can be calculated. After W multiplications by (1−α), the number of significant digits needed for infinite precision is exactly the number of significant digits needed to express (1−α)w, which in the worst case is:
W·d Eq. (10)
where d is the number of significant digits needed to express the decay constant (1−α).
Referring now to
The process starts in block 502, where a number W, determining the length of the mirrored function (1−α)+n is determined as disclosed supra. The process continues in block 504.
In block 504, a number of input signal ƒ(m) samples (602_1-602-W) equal to the number W is accumulated to a memory structure, e.g., a register (602). The process continues in block 506.
In block 506, the recursive calculation of a sample of the output signal (606_1) is carried out. The first sample of the signal in a register position (602_1) is multiplied by −1 in multiplier (604) and is provided on the first input of a block (608). The W—the sample of the signal in a register position (602_W) is multiplied by the last sample of the truncated mirrored function (1−α)+w in a multiplier (610) and is provided on the second input of the block (608). The content of the accumulator (612), is multiplied by
in a block (614) and provided on the third input of the block (608). In the block (608), the three input signals are summed together and provided both to the accumulator (612) and as the output sample (606_1) in a register (606). The process continues in block 508
In block 508, the signal samples in the register (602) are shifted one position to the right, so that the register position (602_1) holds the second sample of the signal, and new sample of the signal is written into the now empty register position (602_W). Concurrently, the signal samples in the register (606) are shifted one position to the right, so that the register position (606_1) is empty. The process continues in block 506.
Because as disclosed supra, both calculations use only a shift of input signal to the right, the process disclosed in
Considering a process for calculation of the convolution of signal ƒ(m) with the two-sided exponential function, the calculation of the convolution of signal ƒ(m) with the mirrored exponential function is identical to the process disclosed in
When the process of calculating a sample of the convolution of signal ƒ(m) with the mirrored exponential function is finished, and the samples in the register 702 are shifted one position to the right, thus the sample of the signal ƒ(m) is shifted from the register position 702_2 to the register position 702_1. Consequently, the calculation of the convolution of signal ƒ(m) with the exponential function may start by providing the sample of the signal in the register position 702_1 to a first input of a block 716. Concurrently with the shift of the samples in the register 702, the samples in the register 706 are shifted one position to the right, so that the register position 7062 is empty and the register position 706_1 contains the sample of a convolution of the signal ƒ(m) with the mirrored exponential function.
Since the sample in the register position 706_1 is to be summed with a sample of a convolution of signal ƒ(m) with the exponential function, the sample in the register position 706_1 is removed from the register and provided to a first input of an block 722. Block 722 sums the sample provided from the register position 706_1 with the sample of a convolution of signal ƒ(m) with the exponential function provided to a second input of an block 722. The output of the block 722 is then written back into the register position 706_1.
It thus follows that the sample in the register position 706_0 is likewise a sum of a sample of a convolution of the signal ƒ(m) with the mirrored exponential function and a sample of a convolution of the signal ƒ(m) with the exponential function. Consequently, the sample in the register position 706_1 cannot be used as the previous term Eα(n−1, ƒ) of the recursion given by Eq. (6). Therefore, block 720 serves as an accumulator, which holds the previous term Eα(n−1, ƒ) provided to the block 722 by the output of block 716, which together with block 718 calculates convolution of signal ƒ(m) with the exponential function. Because the convolution calculation by block 716 and block 718 is identical to the calculation as disclosed in
A person of ordinary skills in the art will understand that the above-disclosed structure may in one aspect be implemented as individual hardware blocks, e.g., a multiplier, an accumulator, a summer. In another aspect, the functionality of the blocks, e.g., multiplication, accumulation, summation, may be implemented as software entities on a processing unit. In yet another aspect, the implementation may comprise a mixture of hardware and software entities.
The above disclosed calculation of fast convolution may be extended to a complex exponential function e−jωt localized around an origin by any of the symmetric functions that can be approximated by an exponential function envelope w(t):
{circumflex over (g)}(t)=e−jωt·w(t) Eq. (11)
where the ĝ(t) is a complex function.
As discussed supra, of a particular interest is a Gabor function, comprising the complex exponential function e−jωt localized around an origin by Gaussian function envelope:
However, since the fast convolution as disclosed herein is by no means limited to Gabor function, the general function expressed by Eq. (11) is used.
Convolution {circumflex over (F)}(n, {circumflex over (ƒ)}) of input signal {circumflex over (ƒ)}(t) with the function ĝ(t) can be expressed as:
{circumflex over (F)}(n,{circumflex over (ƒ)})=Σ{circumflex over (g)}(m−n){circumflex over (ƒ)}(m) Eq. (13)
Substituting Eq. (11) into Eq. (13) yields:
{circumflex over (F)}(n,{circumflex over (ƒ)})=Σme−jω(m−n)w(m−n){circumflex over (ƒ)}(m) Eq. (14)
The term e−jω(m−n) can be rewritten as e−jωme+jωn. The e+jωn can be factored out of the summation:
{circumflex over (F)}(n,{circumflex over (ƒ)})=e+jωnΣme−jωmw(m−n){circumflex over (ƒ)}(m) Eq. (15)
The term {circumflex over (ƒ)}(m)·e−jωm can be recognized as modulation of the signal of interest {circumflex over (ƒ)}(m) by the complex exponential e−jωm:
{circumflex over (ƒ)}n(m)=e−jωm{circumflex over (ƒ)}(m) Eq. (16)
Substituting Eq. (16) into Eq. (15) yields:
{circumflex over (F)}(n,{circumflex over (ƒ)})=e+jωnΣmw(m−n){circumflex over (ƒ)}n(m) Eq. (17)
As with Eq. (15) the right-hand side of the Eq. (17) may be recognized as a demodulation of the term Σmw(m−n){circumflex over (ƒ)}n(m) by the complex exponential e−jωm. The convolution may be separated to the real and imaginary part
{circumflex over (F)}(n,{circumflex over (ƒ)})=e+jωtΣmw(m−n)(ƒnr(m)+jƒni(m)) Eq. (18)
{circumflex over (F)}(n,{circumflex over (ƒ)})=e+jωt(Σmw(m−n)ƒnr(m)+jΣmw(m−n)ƒni(m)) Eq. (19)
Considering Eq. (13) through Eq. (19), the convolution may be separated into several stages as illustrated in
A source 802 provides an input signal {circumflex over (ƒ)}(m) to be convolved with a kernel comprising the function ĝ(t). In block 804, the input signal {circumflex over (ƒ)}(m) is modulated by a complex exponential e−jωm provided by an oscillator 806, yielding on the output of the block 804 complex signal {circumflex over (ƒ)}n(m). The complex signal {circumflex over (ƒ)}n(m) is split in block 808 into a real part 810 and an imaginary part 812. A convolution of two-sided exponential function, approximating the function w(t) with a real part 810 of the signal is carried out in block 814, in accordance with the process implemented by a conceptual structure as disclosed supra. Similarly, a convolution of two-sided exponential function, approximating the function w(t) with an imaginary part 812 of the signal is carried out in block 816. Output signals produced by the blocks 814 and 816 are merged in block 818. To transform the complex signal at the output of block 818, a demodulation, with the complex exponential e+jωn provided by an oscillator 820 is carried out in block 822. The output signal {circumflex over (ƒ)}(n, {circumflex over (ƒ)}) is provided to a sink 824.
A person of ordinary skills in the art will understand that the above-disclosed structure may in one aspect be implemented as individual hardware blocks, e.g., a modulator, a splitter, a demodulator. In another aspect, the functionality of the blocks, e.g., modulation, splitting, summation, may be implemented as software entities on a processing unit. In yet another aspect, the implementation may comprise a mixture of hardware and software entities.
The various aspects of this disclosure are provided to enable a person having ordinary skill in the art to practice the present invention. Various modifications to these aspects will be readily apparent to persons of ordinary skill in the art, and the concepts disclosed therein may be applied to other aspects without departing from the spirit or scope of the invention. Thus, the present invention is not intended to be limited to the aspects shown herein but is to be accorded the widest scope consistent with the principles and novel features disclosed herein. By means of an example, although the processes and structures were disclosed in regards to calculation of fast convolution with an indefinite real-valued signal, a person skilled in the art will understand the processes and structures are equally applicable to complex-valued signals both for indefinite and finite signals. It is noted that a person skilled in the art would use the disclosed calculation even for finite signals because the traffic to memory, i.e., writing and reading samples is halved as compared to the previously used calculated in a forward direction, i.e., for signal samples starting with the first sample and then in a reverse direction, i.e., for signal samples starting with the last sample.
All structural and functional equivalents to the various illustrative logical blocks, modules, circuits, and algorithm steps described in connection with the various aspects described throughout this disclosure that are known or later come to be known to those of ordinary skill in the art are expressly incorporated herein by reference and are intended to be encompassed by the claims. Such illustrative logical blocks, modules, circuits, and algorithm steps may be implemented as electronic hardware, computer software, or combinations of both. Thus a processing unit may be implemented within one or more application specific integrated circuits (ASICs), digital signal processors (DSPs), digital signal processing devices (DSPDs), programmable logic devices (PLDs), field programmable gate arrays (FPGAs), processors, controllers, micro-controllers, microprocessors, electronic devices, other devices units designed to perform the functions described herein, and/or combinations thereof.
Those of skill in the art would understand that information and signals may be represented using any of a variety of different technologies and techniques. For example, data, instructions, commands, information, signals, bits, symbols, and chips that may be referenced throughout the above description may be represented by voltages, currents, electromagnetic waves, magnetic fields or particles, optical fields or particles, or any combination thereof.
Moreover, nothing disclosed herein is intended to be dedicated to the public regardless of whether such disclosure is explicitly recited in the claims. No claim element is to be construed under the provisions of 35 U.S.C. §112, sixth paragraph, unless the element is expressly recited using the phrase “means for” or, in the case of a method claim, the element is recited using the phrase “step for.”
