Detecting peaks in a noisy 2D correlation plot I'm doing pattern matching by means of normalized grey-scale correlation. My input image is noisy and so the peaks in the correlation plot cannot be found just by thresholding.

The original correlation plot displayed as a surface, the peaks are visible, but are not the only local maxima in the picture, I need to find only the ones that are cone-shaped. Note that I know how many I'm looking for.

I noticed that if I calculated the local gradient of the correlation plot, the peaks appear clearly on visual inspection, especially on the gradient orientation plot.

Gradient orientation:

I'm looking for a way to automatically detect the exact position of the peaks (2 in this case, but it is sometimes more) with a one pixel accuracy.

What I guess would help me is a way to identify the center pixel towards which the gradient direction points from all around. The center of the radial gradient, otherwise said.

I have tried to generate an ideal radial gradient such as this one:

and fit it to the image by moving correlation but with no success.

Does someone have an idea? I'm certainly not the first one needing this kind of method, but I couldn't find anything in the literature, am I missing something?

Edit: I found someone who had a similar question here but couldn't find a satisfactory answer.

You could try using the peak-to-sidelobe ratio, i.e. how many standard deviations above the mean is each point in the correlation output.

psr = ${p - \mu }\over\sigma  $

Typically you compute the mean and sigma in a window around each point excluding the region nearest to each point.


Z-Transform vs. Fourier Transform convergence Are there signals for which the Fourier transform is known to exist (perhaps including
singularities) and for which the z-transform does not converge?

There are several signals for which the $\mathcal{Z}$-transform doesn't exist but the (discrete-time) Fourier transform does. One important example is the complex exponential $x[n]=e^{jn\omega_0}$ extending from $-\infty$ to $\infty$. As a consequence, also sinusoidal signals only have a Fourier transform but no $\mathcal{Z}$-transform, as already mentioned in GKH's answer.

Another important group of sequences with only a Fourier transform but no $\mathcal{Z}$-transform are impulse responses of ideal frequency-selective filters, such as ideal low pass, high pass, or band pass filters. The same holds for the impulse responses of the ideal Hilbert transformer and the ideal differentiator.

In general, two-sided signals with a constant envelope (such as the complex exponential and sinusoids), and two-sided signals that decay only as $1/n$ (such as impulse responses of ideal filters) don't have a $\mathcal{Z}$-transform, but they can have a Fourier transform. The corresponding expressions for the Fourier transform may contain Dirac impulses, or they are discontinuous.


Find Fourier series $f(t) = e^{jx t}$ , $−\pi < t < \pi$ I need to find the Fourier series of the $f(t) = e^{jxt}$ ,  $− \pi < t < \pi$
What will be the first step to solve it?

There is no need for you to find Fourier series just compare with equation $x(t)=\sum_{-\infty}^{+\infty}C_n*e^{j*w_n*t}$ we will know $C_n,w_n.$


Time series of Poisson Process I want to run a numerical simulation of a system (eg first order low pass filter) that is stimulated with a photocurrent-like signal. My preferred language is Python.

I'm failing to understand how to generate a proper time series for the photocurrent with a consistent noise power, especially when the photocurrent becomes very small.

Let's assume:

simulation interval time is dT
the mean photocurrent is I
q is the elementary charge
T is the time series length

In my understanding there are two main ways to model the photocurrent time series.

A) I calculate the average number of electrons and use this number to generate a poisson time series that I convert back into current

N = I * dT / q # mean number of electrons in interval dT
N_ts = np.random.poisson(N,int(T/dT)) # electrons time series
I_ts = N_ts * q / dT # resulting photocurrent time series

The problem I have with this is that if I becomes small enough then average number of electrons per dT is very low and when I draw samples from the distribution it behaves like a quantum current (some samples will have 0,1,2,3,.. electrons). Is this quantization a corrected model or is it a limitation of it? What is the meaning of dT in this case? Is it some sort of averaging? Should the dT be increased to avoid the quantization? If I increase dT, how much should it be increased? Does a change in dT change the noise power of a Poisson time series? How is dT related to the noise bandwidth of the signal? How can I rescale the noise after a dT change to maintain the same input noise power?

B) A second approach could be to avoid the photon counting and directly model a photocurrent time series drawing samples from a normal distribution with $\mu_I = I$ and $ \sigma_I = \sqrt{2qIB} $
However, what is the correct bandwidth to use in this case? Is it 1/(2*dT) because I'm assuming the sampling frequency is two times the noise bandwidth? A problem with modeling using a normal distribution is that if dT=1E-6 and I<1.602E-13 then $\sigma_I=\mu_I$, which means that if I draw samples from such distribution I might end up drawing negative currents. The normal distribution is not a good approximation of a Poissonian distribution for small currents. Should I use a lognormal distribution instead?

Are these approaches fundamentally correct? How do they differ from each other? Is there an equivalence between the two? Does Physics come actually into play with quantization when considering very small currents or is it an artifact of the modeling?

A)

The quantization of the photocurrent is an actual physical phenomenon.

I'll use $\Delta T$ in place of your $dT$. It is not a variable describing the actual physical process but rather the time step of your approximate model. Effectively you are filtering the continuous-time Poisson point process by a continuous-time "moving average" filter with a rectangular impulse response of length $\Delta T$ and sampling the filtered signal uniformly with a sampling period $\Delta T$. The filter acts similarly to a low-pass filter, but is not very effective at band-limiting the process. Therefore you have aliasing and it becomes difficult to interpret your results. If you don't know what aliasing means, have a look at the question: "explanation of Hybrid systems?" and its answers for an introduction on sampling of continuous-time signals.

Sampling the RC-filtered Poisson Process

Non-uniform sampling: The impulse response of a first-order low-pass filter is a right-sided exponential function. It would be easy to convolve the Poisson point process directly by it. Non-uniformly sampling the filtered signal at the event times would give an unambiguous representation of the signal that allows exact reconstruction of the signal from the samples. This would not be sampling of a band-limited signal but would take advantage of the special form of the signal. See math later below.

Uniform sampling: If the cut-off frequency of the first-order low-pass filter is low enough relative to the sampling frequency, then uniform sampling of the filtered signal does not cause much aliasing and you will have quite a clean discrete-time signal (time series) to work with. When time stepping, the signal decays exponentially from the previous sample to the current sample, unless there have been events within the last time step. In that case you'd calculate the contribution of the filtered events and add that into the current signal value. That's the beauty of the right-sided exponential impulse response: it allows to forget when the past events occurred, combining their contribution to a single value that keeps decaying exponentially. The drawback of uniform sampling is that if there are two or more events between successive samples, then their exact times cannot in general be reconstructed from the samples of the filtered signal.

There are two easy ways to generate the event times of a Poisson process: 1) The time between successive events follows an exponential distribution, so you just need a pseudo-random number generator with an exponential distribution, and generate the times of the events in order until you are at the simulation end time. The time of an event is obtained adding the exponentially distributed random number to the time of the previous event. Note that it is only a coincidence that also the low-pass filter impulse response is a one-sided exponential function. 2) Generate a Poisson distributed random number to get the number of events over the simulation length, pick that many event times from an uniform distribution that is as long as the simulation and sort the event times. Or you can generate event times for each sampling interval individually.

In this presentation I will not use units, for my own convenience. The impulse response $h$ of a resistor-capacitor (RC) filter, as function of continuous time $t$, is of form:

$$h(t) = \begin{cases}0&\text{if }t<0,\\\displaystyle a e^{-at}&\text{if }t\ge0,\end{cases}\tag{1}$$

Figure 1. RC-filter impulse response (Eq. 1), with decay parameter $a = 1$.

where $a$ is a positive constant that determines the decay rate. If this is convolved with a realization of a Poisson process with $K$ Dirac impulse events at times $p_k$, $0 < k \le K$, then we get the signal $y$ of which you can find some literature under the keyword filtered Poisson process:

$$y(t) = \sum_k^{K-1}h\left(t-p_k\right)\tag{2}$$

Figure 2. A realization of an RC-filtered Poisson process (Eq. 2), with $a = 1$.

We can set up a sample point at time $t_n$. If we know the signal value  $y(t_n)$ at that point, then we can forget the events at or before that time and still reconstruct $y(t)$ at later times $t \ge t_n$:

$$y(t) = y(t_n)\,e^{-a(t-t_n)} + \sum_{k=0}^{K-1}\begin{cases}0&\text{if }p_k \le t_n,\\h(t-p_k)&\text{if }p_k > t_n,\end{cases}\quad t \ge t_n,\tag{3}$$

or in a less formal short-hand notation:

$$y(t) = y(t_n)\,e^{-a(t-t_n)} + \sum_{k,\,t_n < p_k}h(t-p_k), \quad t \ge t_n.\tag{4}$$

This is a property unique to the one-sided exponential function $h(t)$. The effect of an event at time $p_k$ is already included in $y(p_k)$ with the way Eq. 1 is currently defined. We can make explicit also that because the impulse response is causal, the events happening after time $t$ do not affect $y(t)$:

$$y(t) = y(t_n)\,e^{-a(t-t_n)} + \sum_{k,\,t_n<p_k\le t}h(t-p_k), \quad t \ge t_n.\tag{5}$$

You can do non-uniform sampling at the event times by setting $t_k = p_k$, or uniform sampling with sampling period $\Delta t$ by setting $t_n = t_0 + n\,\Delta t$. In your data structure you'd just store $t_n$ (no need in case of uniform sampling) and $y(t_n)$. Python example:

import numpy as np
import bisect
import matplotlib.pyplot as plt

sim_length = 10  # Simulation length
rate = 100  # Event rate per time unit
a = 1  # Filter decay constant

K = np.random.poisson(rate*sim_length)  # Number of events
N = K + 1  # Number of samples. Always have one sample at simulation start
ty = np.empty(N, dtype = [('time', float), ('value', float)])  # Sample times t and signal values y
ty['time'][0] = 0
ty['value'][0] = 0  # Start signal at value 0
ty['time'][1:] = np.random.uniform(0, sim_length, K)  # Event times
ty[1:].sort(order='time')  # Sort events by time
for k in range(K):  # Sample signal at event times
    ty['value'][1 + k] = ty['value'][k]*np.exp(-a*(ty['time'][1 + k] - ty['time'][k])) + a

def sample_y(ty, t):  # Sample signal at arbitrary time t
    k = bisect.bisect(ty['time'], t) - 1
    return ty['value'][k]*np.exp(-a*(t-ty['time'][k]))

# Create a time series of the RC-filtered Poisson process
t = np.arange(0, sim_length, 1/10000)  # Delta t = 1/10000
y = np.zeros(t.size)
for k in range(t.size):
    y[k] = sample_y(ty, t[k])

# Plot
plt.plot(t, y)
plt.xlabel('t')
plt.ylabel('y(t)')
plt.show()

Figure 3. RC-filtered Poisson process plotted by the above Python script.

Figure 4. RC-filtered Poisson process reparameterized with rate = 1 and a = 10.


What Is the MATLAB `imreconstruct()` Useful For? Could some one tell reasons to use the MATLAB imreconstruct function when processing images? I already have studied the topic and what the function actually does to the image. However, I cannot figure what we would use it for. What output should we expect from calling doing this function?

It could be useful in the field of OCR as a pre processing step.

Think about badly scanned data.
You'd like to convert it into binary image.
So the first step would be applying some kind of thresholding.
Since no thresholding is perfect, There will be some "Holes" / "Gaps" within the text.

Closing those "Holes" / "Gaps" can be done using morphological operations.

See example here: StackOverlfow - Adaptive Thresholding - Implementation of the Minimum Error Thresholding Method.


How to calculate bitrate of PAL? I have to calculate the bitrate of PAL format $\ (576*720) $. I have the following data:

$\ 25 $ frames per second
Interlaced
Interior representation is: $\ 4:2:0 $
$\ 8 $ bpp

I have tried the following formula but I am really not sure at all if it is correct because I have not clearly understood the difference between interlaced and non-interlaced video and how this difference changes the calculation fo the bitrate.

This is what I tried:

$\ (576*720)*(1+1/4+1/4)*8*25 $

This part: $\ (1+1/4+1/4) $ is due to the specific representation $\ 4:2:0$

In the early days of analog television broadcast with CRT display receivers, it was found that instead of sending a whole picture frame information (576 PAL video lines) at once at a frame rate of 25 Hz, it was visually better to send half the lines as even field (0,2,4,...,574) and odd field (1,3,5,..,575) at 50 Hz. This is called interlaced video 576i for PAL and 480i for NTSC. If you send all the lines at once (and at half the rate of 25 Hz) then it's called progressive and denoted as 576p and 480p. These days they are used for 720 or 1080 etc line HD or FHD or 4k formats in the same line.

The reason for better visual characteristic is due to the fact that faster field rate yields smoother motion and yet reduced vertical resolution is compansated by several factors.


How to shift FFT spectrum by half a bin? I am trying to find a way to shift the outputs of the FFT (any algorithm implementation, currently using radix-2) by half a bin with the same bandwidth.  To put it another way, say you have a 512pt FFT sampling data for a bandwidth of 100 kHz. The frequency bins are going to fall ~390 Hz apart (bandwidth) For example 195 to 585 is effectively the band pass range for bin 2 which happens to fall at the center frequency of 390 Hz. Bin 3 is from 585 to 975 with Fc being 790 Hz.  Say I want Fc to be 585 instead of 790. This represents a shift of all the FFT outputs by half a band.

I do not want to double the size of the FFT to get the result. I do not want to multiply the signal in time domain prior to FFT by an offset frequency as I don’t want to deal with extra products within bandwidth of the overall FFT.

Is there a way through modification of the twiddle factors to get the desired effect?

Well, if you're willing to perform multiplications in the time domain, you can cyclically shift the whole spectrum in the frequency domain using:

$$x'[n] = e^{j\pi \frac{f}{F_s/2}n} x[n]$$

Where $x[n]$ is your original input sequence, $F_s$ is your sample rate, and $f$ is the amount of frequency shift you desire.

The DFT of $x'[n]$ will have a spectrum that is cyclically shifted by $f$ Hertz, compared to the DFT of $x[n]$. And there will be no "extra products within the bandwidth of the overall FFT", if by that you meant mixing products, from say multiplying by a cosine to effect a shift.

I suppose you could merge that complex exponential factor into your twiddle factors of a custom FFT implementation, but I would advise against that.  Page 19 of these notes shows how that merge would happen: http://eeweb.poly.edu/iselesni/EL713/zoom/dftprop.pdf if you are dead set on modifying the twiddle factors instead of a separate multiplication.


Gaussian Filter Close to Image Border When applying Gaussian filters close to the borders of an image, algorithms usually zero-pad or replicate/mirror/wrap the borders. This is not good enough for my case, so I wonder if there is something better out there.

As the idea of smoothing is to absorb the information of surrounding points, I wonder if it would make sense to use a different kernel in the border cases. This kernel would only take values from available surrounding points, e.g. if my kernel is something like:

\begin{bmatrix}
    0.0113 & 0.0838 & 0.0113 \\
    0.0838 & 0.6193 & 0.0838 \\
    0.0113 & 0.0838 & 0.0113
\end{bmatrix}

To calculate a pixel on the top margin (but not close to the corner), I would use

\begin{bmatrix}
    0.0938 & 0.6933 & 0.0938 \\
    0.0126 & 0.0938 & 0.0126
\end{bmatrix}

The latter matrix was calculated by cropping out the part of the kernel that would be out of the image (the top row) and renormalizing. Is there something wrong with this approach, or an alternative?

Here's some code to explain my idea better:

kernel = [0.0113 0.0838 0.0113; 0.0838 0.6193 0.0838; 0.0113 0.0838 0.0113]
kernel2 = [0.0938 0.6933 0.0938; 0.0126 0.0938 0.0126];
M = [1 2 3; 4 5 6]; % my image

% application of kernel to point 2 in M would provide the value:
res1 = sum(sum( kernel(2:3,:).*M )); % 0 padding, equals 2.1058
res2 = sum(sum( kernel.*[1 2 3; M] )); % row replication, equals 2.3186

% application of kernel2 to the same point:
res3 = sum(sum( kernel2 .* M )); % equals 2.3568

Note that in using kernel2 I don't extrapolate a row of [0 0 0] or [1 2 3], which would underestimate the point.

Pixels outside the image borders must be extrapolated.
Now, you need to chose the model of your extrapolation.

For instance, if you're working within the Discrete Fourier model a periodic extrapolation is the native choice.
Usually the most "Eye Friendly" solution is replication based on "Nearest Neighbor".

Your choice is as valid as any choice, you just need to figure which model it satisfies and if that model suits you.

Yet it has 1 noticeable caveat, it makes your kernel Spatially Variant, which means, without a caring treatment it will hurt performance.
The right thing to do is isolate the pixels which should be carefully treated because of the boundary condition and work very efficiently on the rest.

Tim Holy talks about it in the video Image Representation and Analysis (Go to 42:30).


Nyquist sampling rate Please i need help in understanding Nyquist sampling rate. What is the implication of sampling at a lower rate than the Nyquist rate, at exactly the Nyquist rate and at a rate higher than the Nyquist rate.

Implications of the Shannon-Nyquist sampling :

a. Sampling rate below Nyquist rate: Aliasing happens: shifted spectra of the sampled signal will overlap; i.e., high frequency signal content will be aliased into lower frequency regions of the spectrum and low frequencies into higher, therefore, irrecoverable misinterpretations happen about the authentic (true) content of the sampled signal spectrum;

b. Sampling rate equals Nyquist rate : Depends on the signal type. If there is no impulse at the Nyquist frequency, or the spectrum is zero at the Nyquist frequency, then this will in principle satisfy Nyquist sampling rate, but requires an ideal anti-aliasing filter to be practically implemented and thus usually avoided and a slightly higher sampling rate is preferred instead.

c. Sampling rate larger than Nyquist rate : Satifies Shannon-Nyquist sampling theorem, allowing perfect and unique recovery of the original continuous-time signal from its samples. However, going too high, unless necessary will increase system cost. So practically it's best to limit the sampling rate to a satisfactory minimum or an achievable maximum.


How to account for phase and shift when composing a new signal based upon 2 or more recordings? I have 2 signals that were independently recorded, and want to produce a new signal that would represent what would be recorded if both signal sources were emitting at the same time together. Specifically my case is that I have a recording of dolphins, and a recording of beluga whales and although both recordings were produced independently I would like a sound file with both biologics in the signal.
This post on a MATLAB forum link recommends to multiply the time ordered values together $x_1(t) \times x_2(t)$ to produce a new signal as a result of 'mixing'. This blog post link describes the addition of noise as an additive process between the signals. A question on this forum relevant question asks about the composition from multiple signals where the addition is prescribed. But these approaches do not account for phase and I would assume that this would be necessary. It also does not account for the fact that the time of the recording for the onset of each signal greatly affects the alignment of the signals where there is irregular activity. There is quite a large amount of variation in the acoustic activity where the start time changes whether these animal sounds overlap or not.
The approach I currently use, is to take the FFT of each $x_i(t)$ separately and then generate signals using iFFT from each FFT output. This produces separate signals that I then add together their values at every time point.

Is there a better way to do this?

In order to create a more realistic environment would I have to segment the different states of the signals (biologic activity and non-active states) in order to simulate the randomness in their activity. Would that be an approach seen elsewhere? In general a more simulation like approach to account for the stochastic changes in the state of the system where phase then is taken into account?

In EE, waveform mixing is used for wave modulation. It is implemented via mixer circuits whose operation principle can be loosely described as multiplication of currents. It does not apply to your scenario of SPL data acquisition and processing: you have no waves that modulate one another.
Your sound sources being independent and also removed at distances larger than a sound wavelength from each other (1.5 m at 1kHz in water), you can safely ignore phase effects and, specifically, acoustic interference. You can model the composite signal SPL as a linear combination of SPL's of your prerecorded signals.


NORM_HAMMING2 vs NORM_HAMMING I'm trying to use ORB features with a cv::BFMatcher. ORB is based on binary descriptors so HAMMING must be the choice, but I'm not sure I understand the difference between NORM_HAMMING2 vs NORM_HAMMING and witch one is the right choice.

I made some test but I couldn't find any reliable  difference.

From BFMatcher constructor documentation:

  NORM_HAMMING should be used with ORB, BRISK and BRIEF, NORM_HAMMING2 should be used with ORB when WTA_K==3 or 4 (see ORB::ORB constructor description)

And ORB constructor documentation:

  WTA_K - [...] Other possible values are 3 and 4. [...] Such output will occupy 2 bits, and therefore it will need a special variant of Hamming distance, denoted as NORM_HAMMING2 (2 bits per bin).

In other words the right choice is the one that matches your ORB configuration (NORM_HAMMING when WTA_K==2, and NORM_HAMMING2 when WTA_K==3 or WTA_K==4), otherwise you'll wind up with incorrect norms and correspondingly poorer/inconsistent matching performance.


How to match a drawn Chinese character to a database of stroke animations I have a database of Chinese character stroke animations. Each entry is basically a set of points which shows how and in which order a character should be drawn.

To make lookup easier, I would like to allow the user to manually draw a character with the mouse, then my software would propose some possible matches (like it's done on the iPhone).

I searched on Google but didn't find anything relevant. I'd like to know what algorithms or research paper I need to look at. I guess it has something to do with shape recognition, except that the order in which the shape (Chinese character) is drawn is also relevant. Any suggestion?

You were looking for scholarly articles, I am surprised you did not find anything. Here's what a quick search through Google Scholar with keywords path stroke recognition and handwriting stroke recognition got me:

Plamondon, R.: Online and off-line handwriting recognition: a comprehensive survey: Looks like a good place to start. An explanation of digital handwriting representation, and handwriting recognition techniques (with and without the trajectory data of the stroke)
Does not deal specifically with Chinese characters, but I think it might be applicable.

Hideo Ogawa, Keiji Taniguchi: Thinning and stroke segmentation for handwritten Chinese character recognition: It looks quite old, but is dealing with stroke segmentation, or to quote "a way of breaking down a character to a set of consecutive partial strokes".
It is the reverse of what you want (you said you already have info about strokes), but maybe understanding the stroke representation might help you get on the right way.

Koschinski, M.: Segmentation and recognition of symbols within handwritten mathematical expressions: Again, not Chinese characters, but mathematical (which is still more "exotic" than regular text). Nice because it proposed a technique based on probabilities.

Cheng-Lin Liu, In-Jung Kim, Jin H. Kim: Model-based stroke extraction and matching for handwritten Chinese character recognition: I think the title says it all. Looks like exactly what you want, and promises promising results.

One more thing to add: The most recent reference in the above bunch is as old as 2000. That indicates that what you want to do is probably a well-researched area, with tried out state-of-the-art techniques. It should not have been hard to look for reference material.
Just to clarify, I did not actually read any of these articles. I just did my usual preliminary research step: go to article databases (google scholar, mendeley, ieeexplore...), type in a few keywords. Open the ones with promising titles, close one with completely uninteresting abstracts. Read the remaining abstracts more carefully, and select the ones that seem relevant. Having said this, not all the articles I proposed might be relevant, but judging from experience, they most probably are.
I hope this helps. Also, if you have done any research in the meantime, please edit it in your question and share with the community. It looks like an interesting topic.