2D Gaussian Seperation into 1D Gaussian components I have been reading this http://homepages.inf.ed.ac.uk/rbf/HIPR2/gsmooth.htm and I am having trouble to perform the seperability.

So I kinda did it in paper. Lets say y Gaussian function is G(X,Y), then seperating them will become G(X)G(Y), and then I will need to calculate the 1D component for X and 1D component for Y. Then I can pass over my image twice using the two components each time.

I do have a couple of questions though (one of them is more general):

1 - When I calculate the 2D Gaussian why do I need to normalize it by dividing each pixel of the template with the sum of the template?

2 - Do I need to do the same for the 1D components aswell?

3 - In the link it says:

  Figure 4 shows the 1-D x component kernel that would be used to
  produce the full kernel shown in Figure 3 (after scaling by 273,
  rounding and truncating one row of pixels around the boundary because
  they mostly have the value 0. This reduces the 7x7 matrix to the 5x5
  shown above.). The y component is exactly the same but is oriented
  vertically.

Why scale the 1D component by 273? 273 is the sum of the 2D template not the 1D. And how does he truncate 1 row? The 1D X component is 1 row by itself...

Thanks
I hope my questions make sense.

Normalization is not "required". It only serves to have scale-consistent results, which a not so useful for visualization, but mostly for measurements: if the Gaussian kernel is "sum normalized", the result of the filtering of a constant image is the same constant image. Constant images are "invariant" for smoothing.
You need (or not) to do that for exactly the same reason as above. As it is separable, you can, if you want, normalize first on the 1D sum on the rows, second on the sum onn the columns.
They do not scale the 1D component: the 1D kernel g1d7 $[   0.006  ,  0.061    ,0.242 ,   0.383  ,  0.242  ,  0.061  ,  0.006]$ has almost $1$ average. They do a tensor product of $x$ ($x'. x$) to build a $7\times 7$ kernel, and only keep the center kernel with $273$ normalization.

So,  here is, in pseudo-code with Matlab syntax, a sequence of operations that leads to the $5\times 5$ array g2d5ci from Figure 3.

% 1D 7-tap Gaussian kernel
g1d7 = [0.006,0.061,.242,.383,.242,.061,0.006]';

% 2D 7-tap Gaussian kernel, normalized to 1
g2d7 = g1d7*g1d7';g2d7 =g2d7/sum(g2d7(:));

% 2D 5-tap Gaussian kernel cropped from the central 5x5 part of g2d7, normalized to 1
g2d5c = g2d7(2:6,2:6);g2d5c =g2d5c/sum(g2d5c(:));

% 2D 5-tap Gaussian kernel with integer values
g2d5ci = round(273*g2d5c);

The 1D kernel is not obtained by sampling the Gaussian, but by integrating it over  sub-intervals. For instance, instead of computing $g(x) = \frac{1}{\sqrt{2\pi}\sigma}\exp\left({-\frac{x^2}{2\sigma^2}}\right)$ at $t=3$, you integrate it over the interval $[-3.5,-2.5]$.

The result is in the following figure (pretty close, isn't?):

The Matlab code goes like this:

increment = 0.001;sigma = 1;
time = [-3.5+increment/2:increment:3.5-increment/2]';
data = exp(-0.5*(time/sigma).^2)/(sigma*sqrt(2*pi));
plot(time,data)
data7 = reshape(data,length(data)/7,7);
data7p = sum(data7)';
data7p = data7p/sum(data7p);

figure(1);clf;hold on
plot(time,data)
plot([-3:1:3]',g1d7,'ob')
plot([-3:1:3]',data7p,'xr')
legend('Gaussian kernel','Given values','Recomputed values')
axis tight


Signal encoding identification I've read this signal from a radio remote that opens a gate. This signal is already after AM demodulation, read by using a Saleae logic analyzer:

I don't know much about signal processing, but I would like to emit the same signal with my own transmitter. To do this I would like to know what encoding this signal is using? Is it possible to determine the encoding type only by looking at the signal?

The spikes at the left is the beginning of the signal. I've truncated the rest of the signal because I guess there's no point in showing all of it.

Also could anyone guide me to some good resources about self-clocking signals, usage of low duty cycle clocks in signals, that would help me to decode the message in this signal?

This signal is already after AM demodulation

That means this is only meaningful if the signal actually only contains information in its amplitude; luckily, most simple remotes just use some form of On-Off-Keying or Pulse Position Modulation.

  Is it possible to determine the encoding type only by looking at the signal?

If this really is just a signal that is either on or off, then this looks like on-off keying. But really, your visualization might have hidden a lot of relevant information.

This might be a lot. It might simply be that the receiver uses the clock it has, or derives its clock from the preamble. Then, it could just be on/off, or it could also be manchester encoded.


Corner detector vs feature detector In the literature, I see "corner detector" and "feature detector" are interchangeably used for these methods e.g. Harris, Shi-Tomasi, FAST, etc.
What is the difference between a corner detector and a feature detector? Are corner detectors (Harris, Shi-Tomasi, etc.) one kind of feature detector?

A couple of confusions here.

Features refer to some form of lower dimensional descriptors that explain a (potentially local) region of interest. They are useful in converting the appearance into certain signatures that are easier to handle / more robust than naively using pixels.

What you are trying to ask is probably keypoint detection. This term is used to refer to certain points of interest in the image, that are interesting or worthy to extract features at. Indeed, this is a more general term for corner detection. Most of the time (especially nowadays in the context of deep learning), keypoints need not be some geometrically salient corners but rather points that are critical in describing the image/signal as a whole. The state of the art methods learn this notion of saliency and critical-ity.

In fact, the term feature detection should probably refer to the combination of two steps: first detecting the keypoints and subsequently extracting descriptors on those locations.
Hope this makes sense.


Expressing 2D Gabor transform as 1D spectrum The Gabor transform (Short-Time Fourier Transform with Gaussian window) is represented as $U(t,\omega)$, where $t$ is the time, and $\omega$ is the frequency.  The author of a monograph (Seismic Inverse Q filtering, pg. 190) is asking for the Gabor spectrum to be transformed from 2D to 1D in the following fashion.

Define $\chi = \omega t$
Take the 2D amplitude spectrum $|U(t,\omega)|$ into a 1D amplitude spectrum as $|U(\chi)|$

The issue that I am struggling with is that $\chi$ is the product of $\omega$ and $t$.

The following Matlab code demonstrates that $\chi$ does not continue to increase monotonically.  Note the "sawtooth" appearance of $\chi$.

So how do I take the 2D signal as a 1D signal?  For each $\chi$, there should be exactly one $|U(\chi)|$.

N0 = 51;
N1 = 1604;
fs = 1.0e3;
f = (fs/2) * linspace(0,1,N0);
omega = 2 * pi .* f;
dt = 1 / fs;
T = (N1-1) * dt;
t = linspace(0,T,N1);

chi = [];
cnt = 1;
for i = 1:length(omega)
   for j = 1:length(t)
        chi(cnt) = omega(i) * t(j);
        cnt = cnt + 1;
   end
end

figure;
plot(chi);
xlabel('\chi')
ylabel('Value')

However, the Gabor spectrum is a 2D matrix, similar to this example, where $\omega$ is on one axis, and $t$ is on the other axis:

Equation 11.1 shows the equation for $A$:
$$A(\tau,\omega) = A_0 \exp\left(-\frac{\omega \tau}{2Q}\right)$$

We see that $A$ really only depends on the product of $\tau$ and $\omega$, therefore by defining a new variable $\chi = \omega \tau$ we can rewrite this expression as

$$A(\chi) = A_0 \exp\left(-\frac{\chi}{2Q}\right)$$

Perhaps this helps.


Estimating Real-Time Vehicle Velocity from low-frequency GPS, accelerometer, gyro and magnetometer in Smartphone I am trying to write a system that would be able to estimate the velocity of a traveling vehicle (usually a car) using a smartphone. Most smartphones have the next sensors:

1Hz Low-Frequency GPS
10-100Hz Accelerometer
10-100Hz Gyroscope
10-100Hz Magnetometer

I would also impose a limitation, that device must be in either portrait or landscape position, relative to vehicle's direction of travel.

Now, the GPS will tell me speed and position once per second, but I would like to estimate the velocity in time between GPS data comes in, up to 10-100 times per second, for which I would like to use other sensors. The speed should be accurately represented, even if the car is decelerating.

For now, I've implemented Kalman filter which fuses together GPS, accelerometer and gyroscope and estimates the velocity. But the biggest problem becomes the deceleration, as in all usual references I could find, only the absolute velocity from all axes is taken into consideration.

How to approach this problem? Am I even correct in using Kalman? How to get correct deceleration? What if vehicle is cornering?

Thank you for your help!

How to approach this problem? Am I even correct in using Kalman?

Yes, a Kalman filter (or a derivative) is correct.  Rotations are involved, and 3D rotations are highly nonlinear, so you'll pretty much have to use an extended Kalman or an unscented Kalman.

The basic approach is to model the vehicle motion as position, velocity, and angle driven by the IMU output, with the phone position estimate being corrected by the GPS.

You can correct the phone orientation estimate with the magnetometer, but then you need to take into account the various compass offsets.

  How to get correct deceleration? What if vehicle is cornering?

That will pretty much fall out of the Kalman solution.  As a plus, the phone's orientation with respect to the vehicle will only matter to the extent that you need to know where the vehicle is pointing, and that should be a pretty easy, fixed transform.


Intuitive interpretation of Fourier transform in beamforming When applying Fourier transform in time domain the signal will be carried to frequency domain. And when we apply Fourier transform in spatial domain the signal will be carried to directivity domain.

$$ D(\theta) = \sum\limits_{n = \frac{-N}{2}}^\frac{N}{2} v(n) \exp\left(\frac{j2\pi}{\lambda} nd \sin(\theta)\right)\tag{1}$$

And let's simplify  $ d = \lambda / 2 $ so

$$ D(\theta) = \sum\limits_{n = \frac{-N}{2}}^\frac{N}{2} v(n) \exp(j\pi n \sin(\theta))\tag{2}$$

I can see the mathematically this is $\rm sinc$ if I take $v(n)$ as a rectangular window. I have no problem with the result.

First $\sin(\theta) $ confuses me extremely  how does this $\sin$ performs within Euler's formula ?
Secondly I can understand clear in time domain - frequency relationship: $\exp\left(\frac{j2\pi}{T}t\right) $. From Euler's formula $\cos$ and $\sin$ waves are created every $\frac{t}{T} $ step from $0$ to $ 2\pi$.  I can not understand this when  I look at beamforming function above (formulas (1) and (2)).
And finally why result is in $\theta$ domain how that happened ?

The equations relate to how the directivity of a uniformly spaced linear array of sensors works.

Let's start with an array:

where $\theta$ is our "look" direction (where we want to steer the array), $d$ is the sensor spacing, and, for our purposes, the angles of each segment of the array are all 0 (i.e. $\theta_1 = \theta_2 = \ldots = 0$) so the array is in a straight line along the $x$ axis.

In that case, the time at which the incoming wavefront hits each sensor, $t_j$ is:
$$
t_0 = 0\\
t_1 = d\sin(\theta)\\
t_2 = 2d\sin(\theta)\\
\vdots\\
t_j = jd\sin(\theta)\\
\vdots\\
$$

The idea is that you delay each signal from each sensor by $t_j$, weight them by $v(j)$ and then sum them you'll get an array response like:
$$
a_{\rm tot}(t) = \sum_{j=0}^{N-1} v(j) \delta(t-t_j)
$$
in the time domain.

In the frequency domain, this looks like:
$$
{\cal F}\left[a_{\rm tot}(t)\right] = A_{\rm tot}(f) = \sum_{j=0}^{N-1} v(j) e^{i2\pi f jd\sin(\theta)} \tag{A}
$$

If you now assume that the speed of the incoming waves is $c = f\lambda$ so that $f = \frac{c}{\lambda}$ we get:

$$A_{\rm tot}(\lambda ) = \sum_{j=0}^{N-1} v(j) e^{i2\pi \frac{c}{\lambda} jd\sin(\theta)} \\
= \sum_{j=0}^{N-1} v(j) e^{i\pi c j \sin(\theta)}
$$
if we assume that $\lambda = \frac{d}{2}$.

This last assumption means that now our $A_{\rm tot}$ is not dependent on the frequency variable, and just shows us how the array response changes with bearing angle, $\theta$.

So, really, while the Fourier transform is used in this derivation, equations (1) and (2) in your question are not Fourier transforms themselves.

I don't believe I've directly answered your questions, but feel free to comment and see if I can remove your doubts.

  if beamforming equations are not Fourier transforms why the result of a uniformly weighted beamformer is sinc function in $\theta$ domain just like Fourier transform ? And also does this means equation 11 in the paper is wrong

So equation 11 of that paper is the same as equation (A) in what I wrote above (I just added the tag).

So, no, the equation is not wrong.

The reason there is a $\rm sinc$-like function is simply the maths of summing complex exponentials.  The function is actually not a $\rm sinc$ but a only like it for small values of the argument in the denominator.

As you can see from (4.5) in the link, it actually looks like:
$$
{\rm asinc(\omega)} = \frac{\sin(N\omega/2)}{\sin(\omega/2)} \approx \frac{\sin(N\omega/2)}{\omega/2}
$$
where the approximation only holds if $\omega/2$ is small.

  But in the formula there is $\sin(\theta)$ instead of $\theta$ how $\sin$ affects this equation?

As @Marcus says:

  using $D(\theta)$ and calling it Fourier Transform is wrong; you could use $D(\sin(\theta))$ and get away with it :)

This just means that the $\theta$ axis is scaled slightly differently and makes it more periodic. See the picture below which plots:

$\require{color}\color{black}{\bf black}$ : $\rm sinc(\theta)$
$\require{color}\color{red}{\bf red}$ : $\rm asinc(\theta)$, and
$\require{color}\color{green}{\bf green}$ : $\rm asinc(\sin(\theta))$.

R Code Below

# 30425

M <- 10
N <- 100
theta <- seq(-pi*5,pi*5,pi/N)

sinc <- sin(M*theta/2) / (theta/2)
asinc <- sin(M*theta/2) / sin(theta/2)
asinctheta <- sin(M*sin(theta)/2) / sin(sin(theta)/2)

plot(theta, sinc, col=1, type="l", lwd = 4)
lines(theta, asinc, col=2, lwd = 2)
lines(theta, asinctheta, col=3, lwd = 1)

legend(-10, 10, c("sinc()", "asinc","asinctheta"), col=c(1,2,3), lwd = c(4,2,1))


What is the best to confirm that the output of a matched filter for a linear FM chirp is correct? I have written a matched filter in matlab to compress a linear FM chirp signal and would like to confirm that the results of my range compression. My question is what is the best way to do this? I know that when looking at the plot of the output of the matched filter you should see something like the image below, but was wondering if there is a more scientific way to confirm that the output of the matched filter is correct.

In the absence of noise, the output of the matched filter is exactly the autocorrelation function of the input signal, delayed in time so that the peak autocorrelation is at the chosen sampling time (cf. the first part of this answer on this forum). I don't recall off the top of my head what the autocorrelation function of a linear FM chirp signal is, but the sharp peak that in your bottom trace suggests that the autocorrelation function is sharply peaked, perhaps like a sinc function with very small null-to-null distance.


PSD - Correct Formula for removing Noise Ok, so finally some things are making sense.

I'm trying to scale the results, which I have done by the following equation:

$$
PS_{rms} = \frac{2 \cdot \left |y _{m} \right |^{2}}{
   S_{1}^{2}
} m = 0.....N/2
$$

Where   $y_{m} = (re*re+im*im)$

This produces the following:

(Sorry about the axis - I have not had time to do these)

Prior to scaling the results, I got the following:

Essentially, so far, so good (I hope!)

Now I want to remove the noise around the signal, in theory, so that the main peaks are just showing and remove all the peaks that are not. So in theory, I would just be left between points (60 - 80).

I have been trying to read, and, understand the "Section 10: Averaging and overlapping" in this paper:
http://www.rssd.esa.int/SP/LISAPATHFINDER/docs/Data_Analysis/GH_FFT.pdf

It assumes that the usual remedy is to take the average of $M$ elements and reduce the standard deviation of the averaged result by a factor of $1 / \sqrt{M}$ however, what would be $M$ in this case? How could I therefore reduce the standard deviation, in order to remove noise from the result of scaling the PSD?

This is the only part to which I am now confused about, I essentially need to reduce the noise levels at the start of the result, and, at the end of the result thus giving me only the important peaks within the signal.

What you need is a digital filter. Digital filters are means to shape or "filter" the frequency response of a signal. Specifically, you need something called a band pass filter which would leave the 60-80 frequency mostly alone, and would mostly get rid of everything else. Snoop around this site or Google them. Digital filter are DSP 101, so there are many examples and tutorials out there.


Tikhonov Regularization Alternate Formulation 2.2 Tikhonov Reqularization

  Tikhonov regularization, named for Russian mathematician Andrey Tikhonov, attempts to fix the issue that arises when the least squares method is used with an ill-posed inverse problem by adding an additional constraint on the minimization:
  $$\min_x\{\Vert Ax-b \Vert ^2_2 + \lambda^2\Vert x\Vert^2_2\}$$
  In the new minimization, the additional constraint is the squared norm of $x.$ $\lambda$ is a non-negative constant decided in advance, acting as a weight on the strength of the additional constraint. Some alternate ways to write Tikhonov regularization include
  $$\bbox[lightgreen]{(A^T A + \lambda^2 I)x=A^{T}b}$$
  and
  $$\min \left\Vert\genfrac{[}{]}{0pt}{}{A}{\lambda I}x - \genfrac{[}{]}{0pt}{}{b}{0}\right\Vert$$

Can someone explain the alternate formulation in green. I am not able to convert this formulation to standard one. I get one Atransposedinverse extra. The second formulation can be easily converted to original formulation, though. Thanks.

Just expand the term to be minimized, take the derivative with respect to $x$ and set it equal to zero:

$$\begin{align}||Ax-b||_2^2+\lambda^2||x||_2^2&=(Ax-b)^T(Ax-b)+\lambda^2x^Tx\\&=x^TA^TAx-2b^TAx+b^Tb+\lambda^2x^Tx\tag{1}\end{align}$$

Taking the derivative w.r.t. $x$ and setting it equal to zero gives

$$2A^TAx-2A^Tb+2\lambda^2x=0\tag{2}$$

from which the given equation follows.


What are the spike removal techniques in matlab without using filters output from wavelet filter using daubechies 4 wavelet.

Using wavelet filter this is the output obtained. The problem is elimination of the spikes at 0.023,0.043,0.063 and so on. Please note: no predefined filters present in simulink and dsp toolbox are to be used. Kindly suggest an alternative method to remove the spikes.

based on your results, you can try a median filter, in Matlab - medfilt1; which takes n-length samples of your series x at a time and computes the median, and apply the median value to each value of x. It is a very common non-linear filtering technique.

try something like

result = medfilt1(signal,6)

you and alter the length of evaluation from 6 to whichever number suits the data.


Sampling rate vs sampling time of FFT I am having some trouble understanding the sampling time and frequency used in the computing of the FFT, and I was wondering if someone could make it clear for me :

The Sampling time: the time in which I took my sample for example 5 minutes, its the difference between the time at which i started taking measures and the time when I was done.
The sampling rate: how many samples do I save per second, which means if I took $N$ sample within $13\textrm{ s}$, my Sampling rate will be $f_s=N/13$.

Is this sampling rate the same as the frequency used when calculating the FFT when creating the frequency axis of the FFT?

The sampling time is the time interval between successive samples, also called the sampling interval or the sampling period, and denoted $T$.

The sampling rate is the number of samples per second. It is the reciprocal of the sampling time, i.e. $1/T$, also called the sampling frequency, and denoted $F_s$.

The frequency axis for the FFT is linked to the number $N$ of points in the DFT and the sampling rate $F_s$. It is defined as $f=k\cdot\frac{F_s}{N}$. With $k$ going up to $N$


Edge detection in low contrast images I have some images which have some objects of my interest. The objects have very low contrast w.r.t. background and I would like to detect the edges of these objects in these low-contrast images. Any suggestions?

An example image is shown below (objects are the rectangles in the right side of the image and the edges are their boundaries).

Thanks.

Why don't you initially apply basic gray scale modificiations such as gamma , contrast etc.

Here is the result of performing such a simple processing:

It's now a more recognizable image at least for a human viewer and I believe the software algorithm will also perform better, provided it was not already doing it...


Why do linear phase filters have symmetric impulse responses? It was given as a fact that linear phase filters have symmetric impulse responses, but I don't see why that has to be true. Can somebody please explain or prove this?

Actually, I think I see why.

$$X(j\Omega) = |X(j\Omega)|e^{-j\theta(\Omega)}$$

$|X(j\Omega)|$ is purely real, and therefore if we take the IFT it is even and symmetric.

$\theta(\Omega)= a\Omega$ since the phase is linear, so $e^{-ja\Omega}$ merely shifts the corresponding even and symmetric magnitude in the time domain, so the resulting impulse response is symmetric about $a$.


Is there an a method to fit a wave created from two wave? I need to get the frequency and amplitude for a wave that consists of multiple function.
for example, if I have a sine curve (created from two sine waves), How can I extract the parameters for this complex wave.
I used the curve fitting code provided in this answer, as illustrated below. A simple sinusoidal function is used in the provided code. so I obtain a fitted curve with a single amplitude and frequency.
Data example:
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline

x = np.arange(0, 50, 0.01)
x2 = np.arange(0, 100, 0.02)
x3 = np.arange(0, 150, 0.03)
sin1 = np.sin(x)
sin2 = np.sin(x2)
sin3= np.sin(x3/2)

sin4 = sin1 + sin2+sin3
plt.plot(x, sin4)
plt.show()

import numpy, scipy.optimize

def fit_sin(tt, yy):
    '''Fit sin to the input time sequence, and return fitting parameters "amp", "omega", "phase", "offset", "freq", "period" and "fitfunc"'''
    tt = numpy.array(tt)
    yy = numpy.array(yy)
    ff = numpy.fft.fftfreq(len(tt), (tt[1]-tt[0]))   # assume uniform spacing
    Fyy = abs(numpy.fft.fft(yy))
    guess_freq = abs(ff[numpy.argmax(Fyy[1:])+1])   # excluding the zero frequency "peak", which is related to offset
    guess_amp = numpy.std(yy) * 2.**0.5
    guess_offset = numpy.mean(yy)
    guess = numpy.array([guess_amp, 2.*numpy.pi*guess_freq, 0., guess_offset])

    def sinfunc(t, A, w, p, c):  return A * numpy.sin(w*t + p) + c
    popt, pcov = scipy.optimize.curve_fit(sinfunc, tt, yy, p0=guess)
    A, w, p, c = popt
    f = w/(2.*numpy.pi)
    fitfunc = lambda t: A * numpy.sin(w*t + p) + c
    return {"amp": A, "omega": w, "phase": p, "offset": c, "freq": f, "period": 1./f, "fitfunc": fitfunc, "maxcov": numpy.max(pcov), "rawres": (guess,popt,pcov)}

yy = sin4
tt = x
res = fit_sin(tt, yy)
print(str(i), "Amplitude=%(amp)s, Angular freq.=%(omega)s, phase=%(phase)s, offset=%(offset)s, Max. Cov.=%(maxcov)s" % res )
fit_values=res["fitfunc"](tt)
Frequenc_fit= res['freq']
print(i, Frequenc_fit)
Frequenc_fit=Frequenc_fit
Amp_fit=res['amp']
Omega_fit=res['omega']
Phase_fit=res['phase']
Offset_fit=res['offset']
maxcov_fit=res['maxcov']
plt.plot(tt, yy, "-k", label="y", linewidth=2)
plt.plot(tt,fit_values, "r-", label="y fit curve", linewidth=2)
plt.legend(loc="best")
plt.show()

I got a fitted sine curve with a single frequency and amplitude as follows:

Is there a method to obtain fitted curve matches with the original one? or open software can be used to get the parameters for the best-fitted wave from the original one?
PLEASE, my understanding of FFT is VERY limited! So PLEASE forgive me for my lack of knowledge.

As a precautionary statement, the FFT will give you coefficients of the frequency components closest to the integer multiples of $1/T$ in Hz where $T$ is the length of the sequence in seconds being used for the FFT. If the components are not on these integer boundaries, then you will get several results (called spectral leakage) where the actual components (if information is known such as in the OP’s case of having just two tones) can be estimated as being in between two large adjacent results.
Further when combining the waveform you must pay attention to the phase as well as amplitude of the result in order for it to match the time domain waveform shown. The FFT will provide both a scaled magnitude and phase. The FFT results are the coefficients for complex frequencies (of the form $e^{j\omega t}$ but if it is known that the waveform is real then you can use the first half of the FFT as representing the coefficients of real frequencies in the form of $\cos(\omega t)$, scaled by $N/2$ where $N$ is the number of samples in the sequence used for the FFT.


Optical realization of the Radon Transform I am trying to reproduce this paper concerning the Optical realization of the Radon Transform, especially the simulation (section 3). Shortly, the experiment is just a "Fourier filtering" with a 4F system.

In the paper, they propose to use an optical element defined as $$T(r,\theta)=\frac{e^{2\pi i\theta r \beta }}{r}\text.$$

I tried this and it doesn't give the correct output. I think this is due to the fact that the output should be a circle representing the Radon transform in cartesian coordinates, meaning a polar symmetric element and the optical element obviously isn't, at least in the standard definition of polar coordinates (centered at 0).

Then I tried to change the optical element with respect to the argument mentioned above:

$$T(r,\theta)=\frac{e^{2\pi i r \beta }}{r}\text,$$

getting rid of the $\theta$ in the exponential so my element is now perfectly symmetric (as in Fig. 5 of the paper). The result:

It represents the Radon transform in cartesian coordinates of a Shepp-Logan phantom (240x240 pixels). $\beta$ controls the radius of the circle and is tuned to match the mathematical representation (reference). In image it looks like this:

As you see, the image is reconstructed but there's a problem concerning the amplitude. I do believe the problem can come from 2 parameters.

First, I should maybe use $\theta$ somewhere, but when I generate it using atan2 it doesn't give the correct result.

Second, change the way I define $r=\sqrt{x^2+y^2}$ in order to match the definition of polar representation used for the usual Radon transform (r being the diagonal of the picture taking information twice and centering it).

Has anyone ever encountered a similar problem in image processing, amplitude missing?
Or maybe is someone an expert in Radon transform (computed tomography or such) and knows how to generate the proper polar representation?
Or does someone know what I should try next?

I actually ran into this same paper and had the same problem.

The paper linked below does a good job addressing the issue. Starting on page 13 the discuss the issue with this optical element.

https://mstamenk.github.io/assets/files/OpticalRadon.pdf


Spectral effects of multiplying a signal by itself? This has to do with creating the synchronization circuit for MSK in a DSP as described in "Minimum Shift Keying" by Subbarayan Pasupathy published in IEEE Communications Magazine.

At the beginning of the circuit he shows signal passing through a "squarer".  The explanation reads:

  "Although the MSK signal $s(t)$ has no discrete components which can be
  used for synchronization, it produces string discrete spectral
  components at $2f_+$ and $2f_-$ when passed through the squarer. (The
  squarer, in effect, doubles the modulation index and produces an FSK
  signal with $\delta_f = 1/T$, known as Sunde's FSK.  Sunde's FSK has 50
  percent of its total power in the line components at the two
  transmitter frequencies.)"

$s(t)$ consists of a signal containing components at $\pm 300 \mathrm{Hz}$ around the carrier frequency.

I'm confused by this because when I multiply the signal by itself it doesn't seem to produce a signal with components at twice the offset frequencies.  Is squaring a signal different than multiplying it by itself (sample by sample, not convolution with itself)?

Thank you

Is squaring a signal different than multiplying it by itself

No, it is not different.  Squaring the signal is multiplying it by itself, sample by sample.

In the frequency domain this is equivalent to convolving the original frequency spectrum with itself.  The convolution "spreads" the signal out.  If the signal is at $\pm300$ Hz around a carrier frequency $f_c$, then you should see energy at $\pm(2f_c + 600)$ Hz, $\pm2f_c$ Hz, $\pm(2f_c - 600)$ Hz, $0$ Hz, and $\pm600$ Hz.


Optimal windowing strategy for filtering live signal I have a stream of 250 samples per second arriving at a pipeline stage which looks something like this simplified example.

class PipeLineStage(...):
    historical_buffer = []
    def _process(self, sample):
        historical_buffer.append(sample)
        altered_sample = ...
        return altered_sample

_process gets called indirectly by the pipeline manager whenever a new sample arrived from the device and the returned value it passed to the next stage.

I want the effect of the function to apply a 50hz notch filter.

So my question is: how should I ideally window/buffer the signal of the incoming samples?

Specifically, How should I:

decide how many samples are optimal to filter out a 50hz signal when the sample rate is 250hz? (with minimal delay for live viewing)
manage the windowing of a continuous signal; worth applying an lfilter over overlapping buffers e.g. samples 1-100, then 50-150, then 100-200 etc? or is it just as good to clear the buffer after each application of the lfilter and have separate/discrete windows, e.g. samples 0-100, 100-200 etc?

In other words, my concern it whether applying a butterfilter with lfilter to discrete windows will introduce discontinuities or other artefacts in the filtered signal significant enough to justify the performance cost of a rolling window.

Does the way I'm thinking about this make sense?

If possible, design your filter, and then do what must be done (buffering etc.) to implement that filter properly. When implementing a finite impulse response (FIR) filter such as: (pseudocode:)

out[t] = 0.5*in[t] - 0.6*in[t-1] + 0.7*in[t-2] - 0.8*in[t-3]

where t is an integer time index, you could do it so that you store a few of the past in values in a buffer, like this, stepping t:

t buf[0] buf[1] buf[2]
----------------------
0 0      0      0
1 0      0      in[0]
2 0      in[0]  in[1]
3 in[0]  in[1]  in[2]
4 in[1]  in[2]  in[3]
5 in[2]  in[3]  in[4]
...

That way you always have in[t-1], in[t-2] and in[t-3] available for calculation of out[t]. Perhaps you were worried about the performance cost of moving data in a "rolling" buffer like this. But that extra work can be avoided by using a circular buffer, the length of which is usually rounded up to the nearest 2^n (integer n) so that wrap-around can be done using binary and (2^n-1). The next write position is stored here in an extra state variable i:

t i buf[0] buf[1] buf[2] buf[3]
-------------------------------
0 0 0      0      0      0
1 1 in[0]  0      0      0
2 2 in[0]  in[1]  0      0
3 3 in[0]  in[1]  in[2]  0
4 0 in[0]  in[1]  in[2]  in[3]
5 1 in[4]  in[1]  in[2]  in[3]
6 2 in[4]  in[5]  in[2]  in[3]
7 3 in[4]  in[5]  in[6]  in[3]
8 0 in[4]  in[5]  in[6]  in[7]
...

If you are using single-instruction-multiple-data (SIMD) acceleration, it can be more straightforward to also write the current in[t] into the buffer before calculation of out[t].

If you have an infinite impulse response (IIR) filter (like Butterworth) you can use another circular buffer for the past out values. This would give an IIR implementation structure called direct form I. There are other structures that may or may not have better numerical stability, see http://www.earlevel.com/main/2003/02/28/biquads/.

If you don't care about the phase of the filter's frequency response, a simple two-pole, two-zero notch IIR filter will do the job of eliminating 50 Hz: (pseudocode:)

params:
fs = sampling frequency
f = notch frequency
r = sharpness, 0..1 excluding 1

init:
z1x = cos(2*pi*f/fs)
a0a2 = (1.0 - r)*(1.0 - r)/(2.0*(fabs(z1x) + 1.0)) + r
a1 = -2.0*z1x*a0a2
b1 = 2.0*z1x*r
b2 = -r*r

loop:
out[t] = a0a2*(in[t] + in[t-2]) + a1*in[t-1]
         + b1*out[t-1] + b2*out[t-2]

For such a simple filter you may want to skip implementing an elaborate indexed buffering scheme.


what does %matplotlib inline mean? so i have this lab in coursera course, well i tried first at Linux terminal at the first command %matplotlib inline, it doesn't work :
( bash: fg :%matplotlib no such job ), even if I apply it like matplotlib inline or matplotlib also doesn't work, I know those lines like IN[19] show when you write ipython in the terminal so even when I write ipython then %matplotlib inline it doesnt work, also for the command import ipython it doesn't work, since I am already in iypthon how could I import ipython

%matplotlib inline is a magic command for IPython that allows you to add plots to the browser interface. It's not for the the terminal.
Sounds like you are using the Python interpreter at the command line, which isn't going to be useful for you as a beginner. Install jupyter notebook and use it instead of the command line interpreter to follow along with your course.
You do need to import IPython into your script (like any other module) if you want to use their functions during your script. This is not the same as installing a dependency from pip where it is downloaded and saved to your python local environment.


Can you represent an audio signal in any other domains besides time and frequency? Maybe a stupid question, but if there is a time-domain representation of audio signal, and also frequency-domain, so is there any other domain that signal can be represented in?

Basically it's all about breaking down the information into several 'bits'.  The actual audio signal is a time varying 'value', however often it is useful to consider it in a different form.  As an analogy, consider the number 256:  depending on what you're doing with your numbers, it might be useful to treat the number as 200 + 50 + 6, or 16 + 240, or 16*16, or maybe as 2^8; there are an infinite number of ways to treat the number, and which one you use depends on what you're trying to achieve.

The 'frequency domain representation' is an example of doing the above break-down, but with a signal rather than a number.  In this case, you're representing the original signal as a sum of sinusoids, all with different frequencies, amplitudes, and phases.  If you add them all together, you get your original signal back.  You could alternatively choose to represent it in a different way, such as with wavelets, or any other approach that may or may not even have a name (yet), if it is useful for what you're doing to it.  Maybe you could break it into 3 second chunks, then rearrange the signal in each of those chunks to be monotonically increasing, and remember how you reordered them.  This sounds a bit ridiculous, but there is an approach to manipulating MRI images with this sort of reordering (it doesn't use the chunking part, but that's just an example of how you can do whatever you like).

One advantage of sinusoids is that, like Nathan Day says, it relates to the way our ears interpret the pitch of sounds.  But a more important reason is that sinusoids are complex exponentials, which are eigenfunctions of linear systems; that is to say that linear systems are much much simpler to analyse if you consider the inputs and outputs as sums of sinusoids.  That is the main reason why Fourier analysis is so widespread and important.

Short answer to your question: there are an infinite number of domains you can represent an audio signal in.  For another popular one, see wavelets.


Butterworth filter for signal processing I have a school project to be able to move a wheelchair according to user instruction. I am currently starting my project by using MATLAB and I have acquire several data for the training purpose. I am stuck at applying the butterworth filter to the signal, here is a sample of my training data $\tt go10.wav$. And here is my code:

[x,Fs] = audioread('go10.wav'); % to read the data
plot(x) ; % plot the data before aprplying the filter
n = 3; % the number of order
beginFreq = 800 / (Fs/2);
endFreq = 2000 / (Fs/2);
[b,a] = butter(n, [beginFreq, endFreq], 'bandpass');
y = filter(b, a, x);
figure;
plot(y) ; % plot the data after applying the filter
p = audioplayer(y,Fs);
play(p) ; % play the record after applying the filter

Here the plotting result,

So my question is,

Did I do it right? about the begin and end frequency? and is my n correct?
And what can you tell me about this butterworth filter? I know that it's trying to focus on the human voice, but what exactly it does?

First of all when you work with MATLAB I suggest you to use a time axis t.

[x,Fs] = audioread('go10.wav'); % to read the data
t=(0:length(x)-1)/Fs;
plot(t,x) ; % Now when you plot you can see the time in seconds

The code seems to be correct. A butterworth bandpass filter, like others bandpass filter keep only the frequencies in the interval from beginFreq and endFreq. It is impossible to make a perfect filter, so there are a lot of filters like bessel, cebicev, butterworth,... each one have some particular properties and the choice of the better filter depends on the application. In sound processing butterworth filter is quite good. The human voice goes from 100Hz to 1100Hz, if you want to keep the voice in your signal you should change your values to:

beginFreq = 100/ (Fs/2);
endFreq = 1100/ (Fs/2);

A passband filter, like this butterworth, keep only the frequencies in the range you have chosen, so if you have noises in your signal that can't are human voices this filter will remove them. However this filter keeps all noises that have frequency in the vocal range.


use of complex conjugate of weight vectors in beamforming literature in beamforming literature, people generally quote the following as a starting point:

$y(k) = \sum_{i} w_i^*\cdot x_i(k)$

Here:

i varies from 1 to N where N is the number of sensors

$w_i^*$ is the complex conjugate of the ith weight vector

$x_i(k)$ is the kth sample received at the ith sensor

y(k) is the output of the beamformer for kth sample

The question is why is the complex conjugate of the weight vector $w_i^*$ used(instead of the complex weight vector $w_i$ iteslf). The literature seems to suggest that it simplifies notation but I can't see how.

This comes from the fact that the convolution operator for discrete-time linear systems is typically defined as the inner product between the input signal $x_i(k)$ in your example) and the system's impulse response ($w_i$). The definition of the inner product for complex vectors includes the conjugation:

$$
a \cdot b = \sum a_i b_i^*
$$

For the above definition, note that there doesn't seem to be consistent agreement over which of the arguments to the dot operator gets conjugated, so be careful when reading the literature. For the application of discrete-time FIR filters (which your example is), the typical convention is to conjugate the vector of filter taps, hence the $w_i^*$ term in the sum.

Note that this convention is not universal in the FIR filtering literature; I have seen papers where the convolution sum for a system with complex-valued impulse response does not include the conjugate operator. However, I prefer to use the form that you showed in your question.


What is an exact measure of sparsity? I am currently working on compressed sensing and sparse representation of signals, specificly images.

I am frequently asked "what is sparsity definition?". I answer "if most elements of a signal are zero or close to zero, in some domain like Fourier or Wavelet, then this signal is sparse in that basis." but there is always a problem in this definition, "what does most elements mean? Is it 90 percent? 80 percenct? 92.86 percent?!" Here is where my question arises, is there any exact, i.e. numerical, definition for sparsity?

"Is there any exact, i.e. numerical, definition for sparsity?" And by numerical, I understand both computable, and practically "usable". My take is that: not yet, as least, there is no consensus, yet there are some worthy  contenders. The first option "count only non zero terms" is precise, but inefficient (sensitive to numerical approximation and noise, and very complex to optimize). The second option "most elements of a signal are zero or close to zero" is rather imprecise, either on "most" and "close  to".
So  "an exact measure of sparsity" remains elusive, without more formal aspects.
One recent attempt to define sparsity performed in Hurley and Rickard, 2009 Comparing Measures of Sparsity, IEEE Transactions on Information Theory.
Their idea is to provide a set of axioms that a good sparsity measure ought to fulfill; for instance, a signal $x$ multiplied by a non zero constant,  $\alpha x$, should have the same sparsity. In other terms, a sparsity measure should be $0$-homogeneous. Funnily, the $\ell_1$ proxy in compressive sensing, or in lasso regression is  $1$-homogeneous. This is indeed the case for every norm or quasi-norm $\ell_p$, even if they tend to the (non-robust) count measure  $\ell_0$ as $p\to 0$.
So they detail their six axioms, performed computations, borrowed from wealth analysis:

Robin Hood (take from the rich, give to the poor reduces sparsity),
Scaling (constant multiplication preserves sparsity),
Rising Tide (adding the same non zero account reduces sparsity),
Cloning (duplicating data preserves sparsity),
Bill Gates (One man getting richer increases sparsity),
Babies (adding zero values increases sparsity)

and probe known measures against them, revealing that the Gini index and some norm or quasi-norm ratios could be good candidates (for the latter, some details are provided in Euclid in a Taxicab: Sparse Blind Deconvolution with Smoothed $\ell_1/\ell_2$ Regularization, 2005, IEEE Signal Processing Letters). I sense that this initial work ought to be further developed (stay tune at SPOQ, Smoothed $p$ over $q$ $\ell_p/\ell_q$ quasi-norms/norms ratios). Because for a signal $x$, $0< p\le q$, the norm ratio inequality yields:
$$ 1\le \frac{\ell_p(x)}{\ell_q(x)}\le \ell_0(x)^{1/p-1/q}$$
and tends to $1$ (left-hand side, LHS) when $x$ is sparse, and to the right-hand side (RHS) when not. This work is published in 2020 as SPOQ: smooth ℓp-Over-ℓq Regularization for Sparse Signal Recovery applied to Mass Spectrometry (arxiv preprint, SPOQ published paper, IEEE Transactions on Signal Processing, 2020 and Matlab toolbox).
However, a sound measure of sparsity does not tell you whether the transformed data is sufficiently sparse, or not, for your purpose.
Finally, another concept used in compressive sensing is that of the compressibility of signals, where the re-ordered (descending order) coefficient magnitudes $c_{(k)}$ follow a power law $C_\alpha .(k)^{-\alpha}$, and the bigger the $\alpha$, the sharper the decay.


Fast approximate optical flow / image shift I need to detect how fast a camera is panning (either horizontal/vertical) to give a warning to the operator to slow down.

The entire image is moving as a block, I don't need an actual direction (although H or V would be a bonus) and I only need an approximate magnitude - ie. trigger if more than 'N' pixels shift between frames.

Images are large and generally uniform low contrast scenes, I don't have any obvious highlights to track. I need to do it in realtime (60fps) and without using all of the CPU.

Niave solution is pick an RoI in the center, find edges, calculate similarity between pairs of frames, shift one of frames left/right/up/down by a pixel, repeat - find minima.

I wondered if there was a smarter solution?

This might be a slow terrible solution, but you could do a FFT-based cross-correlation of subsequent frames and then find the peak to identify the offset between frames.  Maybe only do it on a small subset of the image to save processor cycles.

It would not work with rotation or drastic scene changes from one frame to the next, and there are probably better methods.  This is kind of a "I have a hammer so everything looks like a nail" solution.  I guess this is just like your naive solution, except there's no need for edge detection and the FFT makes it a lot faster than explicitly shifting one pixel at a time.

This question is similar, and no one's suggesting anything other than cross-correlation, so maybe it's not so bad: Using MATLAB to calculate offset between successive images


Why oversample data before using raised cosine filter? In the process of creating a baseband signal, the complex symbols are filtered using a raised cosine filter to perform pulse shaping (I'm following some literature as a guide). The pulse shaping is performed with the code:

%Create the filter
over=16;
pulse=rcosine(1,over,'normal',0.35);
[trash, pos]=max(pulse);

%Oversample the symbols
sig=kron(sym,[1 zeros(1,over-1)]);
sig=filter(pulse,1,sig);
sig=sig(pos:end);

So the sampling frequency of the filter is 16, which requires 15 zeros be added between all of the symbols. What is the purpose of this? How does this not introduce "extra data" into the original message?

In digital communications, you can transmit a sequence of numbers using a train of pulses that don't interfere with each other. A given raised cosine pulse $p(t)$ has zero crossings every $T$ seconds around its peak (for a certain $T$). Let's say you want to transmit the number sequence $\lbrace 3, -5 \rbrace$. You can do this with the 2-pulse train $$s(t)=3p(t)-5p(t-T).$$ If you sample $s(t)$ at $t=0$ and $t=T$, you'll get the values 3 and -5.

Your question is related to implementing this idea in Matlab. To do this, you have to solve two different problems: how to create the pulse $p(t)$, and then how to create the train pulse $s(t)$.

Let's start with the pulse $p(t)$. Since you specify Fd=1, then the number of samples in the pulse between the peak and its zero-crossings is given by over. In your case, you can see that the peak occurs at sample 49 and the first zero crossing at sample 65, or 16 samples apart, which is the value of variable over. So, you have defined a pulse $p(t)$ with $T$ equal to 16 sample intervals.

In order to generate the pulse train $s(t)$, you start with the sequence sym. Then you insert 15 zeros between the elements of sym, and finally calculate the convolution of the spaced sequence and the pulse. The result of the convolution will be the addition of pulses $p(t-T)$, scaled with the values of sym.

Let us return to the example of transmitting the sequence $\lbrace 3, -5 \rbrace$. You start with

sym = [3 -5];

from which you generate sig:

sig = [3 0 ... 0 -5 0 ... 0];

The convolution of sig and pulse is the addition of 3*pulse starting at sample 1, and -5*pulse starting at sample 16. The first pulse has a peak equal to 3 at sample number 49, and the second pulse has a peak equal to -5 at sample number 65. The pulses interfere with each other at most sample times but, crucially, the second pulse is zero at sample 50, and the first is zero at sample 66. So, at those precise times, they don't interfere with each other. By sampling the convolution at those times, the receiver can recover the transmitted sequence.

Try this simple code in Matlab:

p=rcosine(1,16,'normal',0.35);
s=[3,-5];
sig=kron(s,[1 zeros(1,15)]);
ss=filter(sig,1,x);
plot(ss);

Using the data cursor, you can verify what I said above. You can also print those values with:

[ss(49) ss(65)]


Image with specific frequencies I have created a frequency domain image filter and would like to test it against some images where it is known that the frequencies only occur within a given band. Trouble is, I don't know how to go about creating such images.

I have though about making an image that consists of intensities varying about the x and y axes as a function of sine waves of frequencies within the band. I don't know if this would be good enough though. Another method (don't know about the practicality of this) might be to define a matrix in the Fourier domain and then do an iDFT. Does anyone have any thoughts on the best way to proceed?

Your second option is the most common: choose a frequency of interest, set the corresponding coefficient in the Fourier domain then take the inverse transform.
You have to pay attention to the correct scaling of the non-zero coefficient and to the particular FFT coefficient arrangement (continuous frequency in the middle or in one corner of the Fourier matrix).


LDPC for channel coding I am working on OFDM over harsh channel, so the channel coding is an essential thing for achieving a reliable communication.
According to my reading, LDPC is almost the best channel coding we can use for channel coding, but I have a question regarding that kind of coding.
Assume we are using LDCP with rate $1/2$ which means that communication systems uses the half of data rate for coding and the other half as data.
If we have the OFDM symbol have $1024$ sub-carriers, with modulation order $M=4$ = which gives us a total bits of symbol $2048$. In that case, will the LPDC will take the whole $2048$ together? I mean it will take $128$ bits from the whole bits for coding and the rest for data OR it will take group by group, which means for example group of $6$ bits, and then code them $3$ for data and $3$ for parity or coding.
Which one is right ?

According to my reading, LDPC is almost the best channel coding we can use for channel coding, but I have a question regarding that kind of coding.

Mentally, something being "the best" should always instantly raise a mental flag for you, saying "under which conditions, according to which measure".
It is right that iterative LDPC decoders can achieve maximum likelihood performance (i.e. being the best possible decoder), but only when

The code is large and
there's an infinite number of iterations.

While 2. is never fulfilled, there's often a number of iterations after which gains are small enough for people to simply stop and limit complexity.
The first condition, however, is often a pretty complicated to fulfill one:

will the LPDC will take the whole 2048 (bits) together?

That is a small LDPC code (128 bits would really be tiny, and I don't think I've seen any sensible OFDM application do that; least-rate IoT modes might be interested in doing that on the uplink, but that doesn't match the OFDM approach).
At a small size like 2048 bits, less complex codes and decoders might be comparable or even better suited for your use case and error model.  (PS: Have an error model before deciding on the code you use! There's an awesome website, http://pretty-good-codes.org/ , which is sadly offline now, which compares many codes with metrics.)
Try to put more bits into a single codeword if you want to harness the powers of LDPC codes. For example, in DVB-T2 (which is an OFDM system), the short codewords are 16000 bits long, the normal ones are 64800 bits long.

I mean it will take bits from the whole bits for coding and the rest for data OR it will take group by group, which means for example group of bits, and then code them for data and for parity or coding.

Neither. You take blocksize · rate (so, in your proposed system, 2048·1/2=1024, but really, use larger blocks and established LDPC codes) information bits and encode them as one. These are not systematic, usually, so you don't get separate redundancy bit. You get a code word that is blocksize bits long, and doesn't contain the original bits in any structured way, usually. (Systematic LDPCs are usually undesirable.)
You will need to use a decoder to get the original bits from the codeword.


How to explain this BER vs SNR behaviour? I wrote a script in MATLAB to simulate an LDPC created using an algorithm of a colleague of mine. The parity check-matrix is regular and has dimensions $1012 \times 1518$ and the code rate is $R=1/3$.

The channel I'm using is a BPSK one with added Gaussian white noise. The mapping of the bits I used were the values $-1$ and $1$. The decoding method is sum-product (or belief propagation).

To perform the simulation, I used as a guide the Example 2.6 in the 36th page of this paper. I'm getting this curve:

This is really different from what I would have expected. The fact that the function is no monotonically decreasing was very shocking. Also, the fact that there is a BER of approximately $1$ for a positive SNR doesn't make any sense to me.

I'm thinking that, maybe, what is wrong is how I defined the a priori LLRs or the SNR. Given that I used $E_b = 1$ to map the bits, then I used some formulas from the section 3.2 of this paper:

$$R\cdot E_b = E_s$$
$$N_0 = 2\sigma^2$$

So, using the formula mentioned in the Example 2.6 of the first paper, we get that

$$r_i = 4y_iR\frac{E_b}{N_0}$$

I believe that this weird behaviour is coming from that formula, maybe a difference of definitions between the two papers that I didn't take into account. I don't really know, what do you think is the reason of the function having that shape?

Extra info:

The SNR axis in the plot corresponds to $R\frac{E_b}{N_0}.$
I've read here a question that could solve my problem, but I'm not sure if it applies directly. Here, I don't have any information about the data bit rate that is mentioned in that question. I don't really know what that quantity they call $R_b$ stands for, either.

In a BPSK system, the SNR at the output of the receiver's matched filter is $$\text{SNR}=\frac{E_b}{\sigma^2_n}.$$ In this formula, $E_b$ is the energy used to transmit one bit, and $\sigma_n^2$ is the power of the noise at the filter output.

Assuming an AWGN channel, the noise process at the filter input has constant PSD; this value is usually denoted $N_0/2$. The Wiener-Khinchine theorem relates the PSD and $\sigma_n^2$: $$\sigma_n^2=\frac{N_0}2,$$ so the SNR may be written as $$\text{SNR}=\frac{2E_b}{N_0}.$$

Note that the convention is to plot the bit error rate vs $E_b/N_0$ (note the missing $2$ in the numerator). The reason is that for quadrature systems, the total noise power in the receiver is $N_0$.

In a coded system, each information bit is transmitted using several coded bits. For example, a Hamming (7,4) code uses 7 coded bits to transmit 4 information bits; the code rate is $R=4/7$. In gneral, $R=k/n$, for a code that transmits $k$ information bits using $n$ coded bits.

To make a fair comparison between different codes, we usually keep $E_b$ constant. This means that the energy of each code bit depends on the code rate: $$E_c=RE_b,$$ or $E_b=E_c/R$.

In simulation, it is often convenient to set $E_c=1$. Then, we have $$\frac{E_b}{N_0}=\frac{E_b}{2\sigma_n^2}=\frac{E_c}{2R\sigma_n^2}=\frac{1} {2R\sigma_n^2}.$$


Calculating the numerical gradient of a volume in one direction at a time I want to calculate the gradient of an image volume in one direction at a time.

Using the built-in function of Matlab gradient() I can get the ∂F/∂x for an F volume (the differences in the x direction) and I can get all the differences along all other directions let say (∂F/∂y and ∂F/∂z) for 3D volume.

FX = gradient(F)
[FX,FY] = gradient(F)
[FX,FY,FZ,...,FN] = gradient(F)

The point is that I need the gradient in just one direction lets say (along the z direction) and I can't do that using this function (to my knowledge and according to the documentation) plus calculating the gradient in all direction is CPU and memory intensive for large volumes (like my case).

I come across the snippet of code  below that calculate the gradient along one direction at a time, but when I compared the results of this with the built-in function of Matlab I didn't get the same results!

function D = gradient3(F,option)
% Example:
% Fx = gradient3(F,'x');

[k,l,m] = size(F);
D  = zeros(size(F),class(F));

switch lower(option)
case 'x'
    % Take forward differences on left and right edges
    D(1,:,:) = (F(2,:,:) - F(1,:,:));
    D(k,:,:) = (F(k,:,:) - F(k-1,:,:));
    % Take centered differences on interior points
    D(2:k-1,:,:) = (F(3:k,:,:)-F(1:k-2,:,:))/2;
case 'y'
    D(:,1,:) = (F(:,2,:) - F(:,1,:));
    D(:,l,:) = (F(:,l,:) - F(:,l-1,:));
    D(:,2:l-1,:) = (F(:,3:l,:)-F(:,1:l-2,:))/2;
case 'z'
    D(:,:,1) = (F(:,:,2) - F(:,:,1));
    D(:,:,m) = (F(:,:,m) - F(:,:,m-1));
    D(:,:,2:m-1) = (F(:,:,3:m)-F(:,:,1:m-2))/2;
otherwise
    disp('Unknown option')
end

Does this code calculate the gradient properly and the one that comes with Matlab calculate it differently? or is there is a better way to do it?

Illustrations of both methods results on the x direction on a volume and displaying a slice from that volume and inspecting some of the same pixels from both images, they both obviously calculate the gradient but the central values aren't the same!

Using Matlab's gradient function gradient():

Using the above code:

The only difference I can see between your gradient3 function and MATLAB's gradient is that the latter returns the horizontal derivative as the first output, and your code returns as "x" derivative the vertical derivative.

Note that MATLAB arrays are stored such that the first index is vertical, and the second index is horizontal. Therefore, the code under case 'x', D(2:k-1,:,:) = (F(3:k,:,:)-F(1:k-2,:,:))/2 computes the vertical difference, which intuitively we'd call y derivative.

Comparing:

F = randn(20,20,20);
[FX,FY,FZ] = gradient(F);
isequal(FX,gradient3(F,'y')) % returns true

By the way, using convn you can speed up the computation of the finite difference derivative by about a factor 2x:

D = convn(F,[0.5,0,-0.5],'same');

This produces a different result at the edges of the image, which might or might not be a problem for you but is typically unimportant.


Note detection in an audio file. I am trying to find the given note in an audio file. It has been specifically told that each file contains only one note. So i directly took fft and found the corresponding frequency and then mapped it to corresponding note.Out of my test cases some are giving the correct answer, but some are not.

The pitch of a musical note is very often not the same thing as the frequency peak in an FFT magnitude computation.

This is due to all the overtones and harmonics within interesting musical sounds, as opposed to single pure (and boring sounding) sinewaves.  Please investigate pitch detection/estimation algorithms that do more than look for a bare FFT magnitude peak.  A little reading on human psychoacoustics, pitch perception, audiology, and the physics of musical instruments might also help.


Mathematical model of equivalent-time sampling, the resulting unevenly spaced periodic signal, and its interpretation I sample a continuous signal $s(T)$ over time. This leads to $s[t]$, which depends on several factors of which some are (pseudo-) periodic. I am interested in the effect of one such periodic factors and want to look at signal changes during this period under differing experimental conditions. I know the start and end of every period $i$ and said period slowly and slightly varies over time.

Raw signal $s[t]$

My first approach was to estimate the size of the periodic changes by looking at the power spectrum, where under certain conditions the largest peak is indeed around the frequency for the mean period which I know a priori.

Power spectrum $S[\omega]$

For each measurement I know the start and the end time of the current period $i$, say $t_{i,start}$ and $t_{i,end}$. If I calculate the "phase" or relative position in the current period $\varphi(t)=\frac{t-t_{i,start}}{t_{i,start}-t_{i,end}}$ I can reorder my data as such and plot it. Because I have many samples and the frequency of the effect is not correlated to my sampling frequency the phase increase in my reordered data averages out to pretty much linear (so the two images below are almost identical).

Reordered raw data for reordered sample number $s[n]$

Data plotted to phase $s[\varphi]$

Looking around with google I found that this is called Equivalent-Time Sampling (ETS) in some fields. However it is difficult to find a mathematical formulation for this operation, so it is difficult to see really what information I keep and lose from my data.

My question now is broad: is there a formal mathematical way to model ETS (maybe as a filter on a transformed domain?) and how can I make valid observations about the size and qualities (for example waveform) for this effect when using ETS?

Above (when calculating $\varphi$) I implicitly asume that the waveform is stretched if a period is longer, but in reality there might just be a longer delay before a certain impulse is processed again. Might it be better to use an absolute time delay since start of current period? (as $\varphi(t)=t-t_{i,start}$)

Update

So I think the data, once in equivalent time, can also be considered as unevenly spaced timeseries, since $\varphi(t)$ is unevenly spaced due to the differing period lengths. An answer here suggests methods to calculate correlation between unevenly spaced timeseries. With those methods it might be possible to make a comparison between measurements without any filtering and without assuming a certain sampling rate (seems artificially high in the reordered data). Another answer (and this related answer) suggests to use a Kalman filter to smooth the unevenly spaced time series.

Since the physical measurement takes a certain time $t_s$ (leading to the original sampling rate), each sample in $s[t]$ can be seen as a convolution of a sampling function (perhaps a boxcar) $f(T)$ with the continuous physical value function $s(T)$ during the measuring time. I suspect this is different in the questions that are linked above, which regard economical data and might be regarded instantaneous values.

If that is correct, then in equivalent time the highest frequency that is sampled should not be higher than in the real / original time. Unless there is some way to deconvolve the measurements in equivalent time, the sampling rate in both timeframes should be equal. So the high number of samples between the start and the end of the average period should be reduced to match the true sampling rate, before comparing two measurements.

I think that a moving average with subsampling can be performed. This should be set such that the resulting sampling rate equals the original one. Since in equivalent time the samples are unevenly spaced this can not be performed by a filter in the frequency domain. The moving average has to be performed in the temporal domain instead.

Scale phase $\varphi$  to the average period $P$ of all $N$ periods so its unit is equivalent to that of $t$

$t'(\varphi) = \frac{\varphi}{N} \sum_{i=1}^{N} t_{i,end} - t_{i,start}$ to give $s[t']$

Moving average of $s[t']$ for every $nt_s$

$s'[nt_s] = \frac{\sum_{t'\geq (n-1)t_s}^{t'\lt nt_s} s[t']}{\sum_{t'\geq (n-1)t_s}^{t'\lt n t_s} \frac{t'}{t'}}$

Im not sure about the notation since every $t_s$ has a different number of samples and the period length is not equal to an integer number $n$ of $t_s$. In practice you could append $s[t'+P] = s[t']$ and then use $t_s' = \frac{P}{ceil(\frac{P}{t_s})}$, leading to a slightly interpolated result, but at least it is truly periodic.

The resulting waveform has no more detail / higher frequency content than the original signal. This signal can be compared between experiments like any other signal, so at least it is a workable solution.

This might be a poorly constructed method compared to other suggested solutions (Kalman filter, non uniform FFT´s gridding and unevenly spaced time series operators), but its probably easier too. Not sure if I lose any information in the process.


Q value for a high pass filter How to calculate the Q value for a high pass filter ?
Ps: I'm a beginner in signal processing, so please put it in simple terms.

You can look at these set of notes provided by Analog Devices as part of their "Basic Linear Design" seminar, or just look for notes on second order filters.

The transfer function of a second order high pass filter is
$$H(s) = \frac{A s^2}{s^2+\frac{w_0}{Q}s+w_0^2}$$ where $w_0$ is the cutoff frequency, $Q$ is the $Q$-factor, and $A$ is the gain.


MATLAB optimization for a subroutine function I am going to find the optimum values for a subroutine function in MATLAB. the subroutine function is not accessible in a mathematical form and we don't know any thing about it. We have ten input parameters and we find the optimum values for them which minimize the output of the subroutine function.

To be clear, you have a Matlab "black box" function which takes a $10\times 1$ input and generates a function value and you want to find the minimizer. Such optimization problems are hard to solve in the sense that we cannot apply any heuristics or tricks because we don't know anything about the function you are trying to optimize. Moreover, you don't have any gradient information either.

I would start by trying a derivative free method like Nelder-Mead (fminsearch in Matlab). Do you have a box constraint on the 10 variables, i.e. do you know lower and upper limits for each coordinate? If so you can do a grid search. Brute force search in 10 dimensions may be hopeless and take a very long time. Instead you can use coordinate descent i.e. change one coordinate at a time to the next grid value and step in a direction that reduces the function value.

You can also try other heuristic global optimization routines such as genetic algorithm, simulated annealing and particle swarm optimization. If you have access to Matlab's global optimization toolbox, you can try their "Optimization app" graphical interface to quickly test some of these out.


Negative values from High Pass filtering I come from Fluid Mechanics domain so pardon me for asking a naive and fundamental question.

I am trying to perform high pass filtering on a set of positive real numbers (~Energy) using FFT. However, the result I obtain has a lot of negative values. I tried to quantify how many of them are negative and it seems $\sim50\%$ of these values are negative.

I tried looking up for an answer and I saw a general approach used in Image processing domain, where people offset these negative values or use some threshold. But in my case, I can't use either of these techniques because I need to look at some statistical quantities after the filtering process.

So I want to ask why do these negative numbers appear for a purely positive input and if so, are there any solutions for it?

A high pass filter will remove a large portion to all of the 0 Hz component or DC offset from a signal.  That's similar to subtracting the mean average, which for a normally distributed signal will make around half of the signal negative (since it was below the original average mean).

If you want the mean to be the same after a high pass filter, you could always add it back.

If you don't care about the actual mean value, you could also add an arbitrarily large positive offset to make sure the sum was always positive.


How to choose a windowing function when performing a Chirp $\mathcal Z$-transform? I've got some data sampled at 32 Hz. The only interesting stuff is happening between about 0-8 Hz and I've got to do real time batch processing and detect events with a delay of no more than a second. So that means I can only perform FFTs with 32 bins or less, to remedy the lack of resolution I'm using a Chirp $\mathcal Z$-transform (czt in Matlab).

I was wondering if there were special considerations to be made when you choose a windowing function on the original time data before when performing a CZT instead of an FFT?

When we add the window function to the definition of the DFT we get

$$ X(k) = \sum_{n=0}^{N-1} w[n]x[n] e^{-j \omega n} $$
with $\omega_k= \frac{2\pi k}{N}, k\in \mathbb{Z} , $ $w[n]$  the window function and $x[n]$ the data. The multiplication of the data $x[n]$ with the window function $w[n]$ in the time domain corresponds to a convolution in the frequency domain.

The CZT and the FFT are basically just different ways to compute this sum.  Each algorithm has its advantages and disadvantages. The FFT is fast, the CZT is more flexible. If you use both algorithms  to compute the same coefficients $X(k)$ you get the same results. So, in this case, the window choice considerations should be the same for CZT and FFT.

So far I can answer only have of the original question: If the CZT is used to compute the same frequency bins as delivered by the FFT, there is no difference with respect to window selection.

Note1: Are you aware of the issues regarding frequency resolution of different DFT algorithms, as discussed e.g. in this thread ? The frequency resolution topic is the same for "CZT vs. FFT" as for "zero padding of FFT" or for "Goertzel vs. FFT".

Note2: This was the basic message of my original post.

If  the CZT is used to compute the DFT on a finer frequency resolution, things are more complicated. In this thread the effect on the DFT is discussed. To make it short, the CZT interpolates the DFT. This interpolation makes window selection more involved.

If (as intended by the OP) the CZT is used to interpolate the DFT at a finer frequency resolution, the underlying maths remain the same, however the different frequency resolution must be kept in mind.

This is illustrated in a very simple example. The first figure shows the absolute values of a rectangular window. The FFT (without zero padding or similar things) computes the DFT right at frequencies where the window is zero, here marked in orange.

With finer frequency resolution this changes, as shown in the next figure.
 This may (as in the example) introduce additional leakage.


Why Should the Sum of Masks / Kernel / Filter in Image Processing be One? I am studying pyramids (expanding and reducing) and I am a newbie.
Why should the sum of the convolution mask be one when we are using gaussian distribution?

According to When should the sum of all elements of a gaussian kernel be zero? and Why is sum of values of an edge detection filter zero? it should be zero not one.

I am confused.  Any precise explanation is really appreciated.

When you're trying to create a blurring kernel (namely a weighted mean of the data) you keep the mask sum as one.
This preserves the mean of the filtered image (think about the DFT, if the sum mask is 1 its DC gain is 1 which means it keeps the sum of the image which means it keeps the mean value of the image.)

When you're trying to apply a high pass filter you want to have gain of zero at the DC coefficient.
This is done by having the sum of the mask as zero (same reason as above.)


Vector of random numbers with fixed sampling rate How to create a vector of random numbers(in Matlab), when each value is [0,1] .

And the sampling rate should be fixed ?

thanks

Sampling rate doesn't change anything. Simply call rand function:

x = rand(1, N);


How to: pixel packs in an image to representative point coordinates (x,y) Let's have an image (gray-scale or even binary) as shown on the following figure at the left hand side, the goal is to generate a list of points i.e., coordinates in the form of (x,y) for each pack of the dark pixels in the image.
What are the proper image processing tools to do this and where are they available?

Updates:
1)
Here you may find some more details about the problem. (Note the variation in the size of packs)

I may suggest having packs detected to compute the convex-hull boundary for each and then find the representative centroid {see this for details}.

2)
Here is the result produced by application of Distance Transform (suggested by "Libor"). Note my annotations on the figure. The method does not work as it was promising!

3)
Erosion eliminates small packs!

from __future__ import division
from scipy import zeros, ndimage as dsp
from pylab import subplot,plot,matshow,show

img = zeros((30,30))
img[10:14,10:14] = 1
img[16:17,16:17] = 1
img[19:23,19] = 1
img[19,19:23] = 1

subplot(221)
matshow(img,0)

subplot(222)
y = dsp.binary_erosion(img,[[1,1],[1,1]])
matshow(y,0)

subplot(223)
y = dsp.binary_erosion(img,[[0,1,0],[1,1,1],[0,1,0]])
matshow(y,0)

subplot(224)
y = dsp.binary_erosion(img,[[1,1,1],[1,1,1],[1,1,1]])
matshow(y,0)

show()

4)
Well here is a Python (i.e., the language of love :) ) implementation of the labeling idea (also proposed by "Jean-Yves" below):

subplot(221)
l,n = dsp.label(img)
sl = dsp.find_objects(l)
for s in sl:
    x = (s[1].start+s[1].stop-1)/2
    y = (s[0].start+s[0].stop-1)/2
    plot(x,y,'wo')

and the result:

Note that although it is done in Python so quick due to Scipy performance, the background procedure in label function should be an exhausting iteration. This may be considered as a trade-off. So for a while I keep being eager to seek more efficient algorithms. And also note that in the given code above I found the center of geometry so simply while for complex or asymmetric shapes this may cause positioning to be biased. That is it is a work in progress ;).

5)
Here is a complex case (a real image) captured from here on which the labelling proposal applied and you see the results. Note that it took only 0.015 s for whole procedure including labelling and finding the objects. Scipy guys, did very well job, I think. Wow! {right-click on image, click on view image for full resolution}

Just a naive suggestion: do you know about component labeling?

The technique is about finding chunks of "touching" pixels and assigning them a label, e.g. an integer number. You can then interrogate each chink separately, looking for the pixel that share the same label.

In MATLAB, here is the function that does it trivially: bwlabel


Cepstrum calculus disambiguation Learning the cepstrum analysis for speech recognition, I have met two different definitions of cepstrum (for discrete signals):

$F^{-1}(ln|F(x[n])|)$ . That is, the cepstrum is the inverse Fourier transform of the logarithm of the magnitude of the Fourier transform of the signal $x[n]$. So, cepstrum coefficients have the dimension of time.
$F(ln|F(x[n])|)$. That is, the cepstrum coefficients are the Fourier transform of the logarithm of the magnitude of the Fourier transform of the signal $x[n]$. Hence, the spectrum of the spectrum. From here, the name of quefrency to indicate the variable of the cepstrum coefficients.

The first equation seems the most used, and intuitively makes more sense to me w.r.t the source separation for speech signal processing. However, I met the second definition quite few times, causing me a great confusion.

How the conceptually differ these two definitions? And which is the most correct one? Always in regards to speech signal processing.

A common technique for computing the inverse fft is to invert the imaginary part of the input array, perform a forward fft and then invert the imaginary part of the output array. In the case of the Cepstrum, the input to the ifft is real-valued (due to the absolute value function) and also symmetric. In this case there is no imaginary part on the input to invert, and the symmetry implies that the ifft will produce no imaginary components on the output that would need to be inverted.  Hence in this case the fft and the ifft are fully equivalent.


Random signal power spectrum I have a signal $X(t)=\sum_{n=-\infty}^{\infty} Z_n \delta(t-n\tau)$, $Z_n$ is a random variable with equal possibility of $\pm 1$ and I know the power spectrum of this signal is $\frac {1}{\tau}$ from the solution.

I have try to calculate the power spectrum by the auto-correlation and then calulated by its Fourier transform. As following

\begin{align}
R_X (t+\tau, t) &=E\left[X(t+\tau) X(t)\right]\\
&=E\left[\sum_{n=-\infty}^{\infty}\sum_{m=-\infty}^{\infty} Z_n Z_m \delta(t+\tau-n \tau)\delta(t-m\tau)\right]\\
&=\sum_{n=-\infty}^{\infty}\sum_{m=-\infty}^{\infty} \delta(t+\tau-n \tau)\delta(t-m\tau)\\
\end{align}

Then, I have no idea how to calulate its Fourier transfrom to get the solution as simple as $\frac {1}{\tau}$. I don't know which part is gone wrong, appreciate any help.

You have to simplify your expression in order to obtain a Dirac Comb (a.k.a. impulse train).

$$R_X(t+\tau,t)=\sum_{n=-\infty}^{\infty}\sum_{m=-\infty}^{\infty} \delta(t+\tau-n \tau)\delta(t-m\tau) = \sum_{n=-\infty}^{\infty}\sum_{m=-\infty}^{\infty} \delta(t-(n+1)\tau)\delta(t-m\tau)$$

Both $\sum_{m=-\infty}^{\infty} \delta(t-m\tau)$ and $\sum_{n=-\infty}^{\infty} \delta(t-(n+1)\tau)$ are actually the same signal: a Dirac Comb with $\tau$ period. The difference would be the that the ''n'' signal would be time shifted by n+1, while the ''m'' signal is shifted by m. As they are defined in $(-\infty,\infty)$, they are defined in exactly the same instants.

Applying delta's well known property:
$$\delta(t-t_0)x(t)=x(t_0)$$

we can obtain that, as we intuitively know because both signals are the same,
$$R_X(t+\tau,t)=\sum_{m=-\infty}^{\infty} \delta(t-m\tau)$$

From here, we just have to obtain the Fourier Transform of the Dirac Comb, which can be easily found, e.g. in Wikipedia, obtaining that all the coefficients are $\frac{1}{\tau}$.


Equalizer coefficients and channel coefficients Blind channel equalization methods equalize the channel without using the source data and without knowing the impulse response of the channel. Consider the channel to be a single input single output FIR system. The popular methods for blind channelequlization are Constant Modulus Algorithm and Least Mean Squares. Certainly, the equalizer coefficients are the not the estimated channel impulse response. It is not clear to me that once the equalizer weights are estimated, how does one get the estimates of the channel coefficients/ impulse response? I am missing a fundamental step somewhere, can somebody please help how we use the equalizer weights to estimate the channel?

With a blind equalization technique like the constant modulus algorithm (which is often implemented using a least mean squares (LMS) filter as you indicated), you aren't directly estimating the channel impulse response itself. Instead, the signal model is like this:

The receiver observes the following signal:
$$
x[k] = s[k] * c[k] + n[k]
$$

where:

$s[k]$ is the underlying desired signal,
$c[k]$ is the channel's impulse response, and
$n[k]$ is a noise process introduced by the channel (typically assumed to be additive white Gaussian noise)

A blind equalizer operates by taking advantage of some constraint on the properties of the desired signal $s[k]$. For instance, if it is constant envelope, like phase-shift keying, then the constraint could be:
$$
|s[k]| = 1
$$
The CMA-LMS blind equalizer operates in two steps. First, it uses the current length-$L$ vector of LMS filter taps (where $L$ is a design parameter of the filter) $\mathbf{w}[k]$ to calculate its output:
$$
y[k] = \mathbf{w}[k]^H\mathbf{x}[k]
$$

where $\mathbf{x}[k]$ is a vector of the last $L$ input samples to the filter. Next, it uses the calculated output to update its set of taps for the next iteration:
$$
\mathbf{w}[k+1] = \mathbf{w}[k] - \mu \mathbf{x}[k](|y[k]|^2 - 1)y^*[k]
$$

When compared to the typical LMS filter structure, one can identify the error signal $e[k]$ as:

$$
e[k] = (|y[k]|^2 - 1)y^*[k]
$$

You can see here how the error captures how far the filter output signal $y[k]$ deviates from the expected constant envelope of the desired signal $s[k]$.

Looking at this structure, you can see that at no point is the equalizer trying to directly estimate the channel's impulse response. Instead, it's more "given the knowledge that I have of the expected signal's structure, how do I need to adjust my filter taps to enhance that structure in the filter output?" So it's more of an inverse problem; the tap vector will tend to more closely resemble the inverse of the channel impulse response. Blind equalization techniques like this can do a surprisingly good job at doing so (even more surprisingly, techniques like CMA can even work for signal types that aren't constant envelope, like QAM).


Video: Detecting a digital timer and when it ticks Input: A video with a digital timer, possibly several hours long, such as this video. The timer might use any font or color, and might be located anywhere on the video, but it probably won't fill the entire screen.

Output: A list of time indices when the timer starts and stops, it might start and stop multiple times. For example:

Timer 1: Starts off paused, starts at 0h20m32s, stops at 2h30m04s, ends up being paused.
Timer 2: Starts off ticking, stops at 0h04m00s, starts at 0h06m38s, stops at 0h32m01s, end up being paused.

Possible frames from different input video files:

How might this be achieved? In particular, what algorithms could be used to categorize the data? How would could you detect a timer that might be everywhere, of any color or font size?

To OCR or not?

If the timer is about $h$ pixels height, you could segment each frame into lines of height $h/2$, run an OCR on each line, and scan the result with a regexp for anything that looks like a timer, eg. \d+[\.:,;-]\d+[\.:,;-]. This gives you $x=time_{real}$ vs. $y=time_{video}$ data pairs. Sure, most of them are wrong, but imposing the restriction $y=1*x+n$ (unless the video speed was changed) eliminates most of them. Just group all data points according to whether they lie on a line of slope m=1.

However, this works only if there are at least some good data points. If we cannot make any assumptions about the color, font, or size (or even script!), OCR doesn't perform well.

Instead, making use of some intrinsic properties of all timers works much better and is agnostic to the exact font or script used. Let us make the following basic assumptions about a timer:

overlayed on a static background, ie. there's at least one bounding box that contains only the timer and no other changing objects
it remains on a fixed position for the entire duration of the video
a ticking timer is changing periodically
it contains a seconds digit with a frequency of 0.1 Hz
it is not animated, ie. the seconds digit remains static and changes only every second

We can use these assumptions to recognize the timer automatically.

Downsampling

Scanning every frame, pixel, and every color channel is a lot of work, takes time, and also completely unnecessary. Let's say we've got a 2 hour video, 720p60, yuv420p. That's 2*3600 seconds * 60fps * 1280*720 pixels/frame * 3 channels/pixel =~ 1200 billion data values to examine.

Grayscale. The color channels are down-sampled already anyway.
1-0.1 fps is usually enough, the timer's frequency is only 0.1Hz.
A resolution of 320x180 pixels is more than enough, unless the timer had been really small to begin with, that still leaves us with ~20-60 pixels to analyze.

This reduces the number of data points by about 3*600*4*4~=30000, for a total of about 40 million data points.

Locating the timer(s)

This is the hard part.

We employ several methods to check each pixel whether it might belong to the timer, then we combine them to larger areas. The video is split into intervals of length dt, and for each interval,n frames get analyzed. 60 frames per interval works pretty nicely.

Check for changes. At the very least, the timer must be changing constantly. We could use the standard deviation of each pixel, but I found that randomizing the sample order and taking the average difference between adjacent pixels gives better results sometimes. The standard deviation tends to be lower.
Simple check for 10-seconds periodicity. Group pixels into those that are pixels 10 seconds apart of each other, take the average difference for each group, then average all groups. Basically the auto-correlation for $\tau=10s$ This is simple, fast, and tends to work well.
Compute the full auto-correlation (via FFT for speed). Sounds better, but you need to find a way to analyze it. What I found works is checking at 10 seconds, 20 seconds etc. whether there is a peak (higher value than at nearby samples), and checking that there is no higher value between two adjacent peaks.
Check step-behavior. Take samples every 0.2 or 0.1 seconds and check when the pixel changes. As each digit remains on screen 1 second, it should do so every 5 or 10 frames. Some digits may look alike, but the difference between two changes should be an integer multiple of 5 or 10.

I found that for the most part, combing the first two methods yields quick, fast, and accurate results.

Now we've got a mask of pixels that possibly belong to the timer. Apply a close operation on them (with a threshold of ~0.5), label connected areas, compute the bounding boxes, eliminate very small or very large ones, and areas that haven't got many pixels that were classified as belonging to the timer.

After doing this for all intervals, group bounding boxes together that are close to each other, eliminate groups with only a few members, and take only those groups with the most members. These are most likely the bounding boxes for the seconds digit of the timer. Average each group and return the result.

(Left: Mask, Right: Image with Mask)

(Top: Raw image, 2nd: checked for changes, 3rd: simple check for 10s periodicity, 4th: auto correlation diagram analyzed, 5th: all masks combinded, Bottom: Bounding box computed)

Note that the close step is not shown above.

Scanning for activity (start/stop)

Now that we found a timer, or rather the seconds digit, we scan the video to analyze when it is ticking and when not. If it is ticking, the digit will keep changing.

extract n frames at an interval of dt , crop to the timer rectangular area
for each pixel, take the squared difference between adjacent frames, and average over all pixels
threshold to 0 (stopped) or 1 (ticking), a good value is about 1.0-1.5
apply a gaussian smooth of about 2-4 samples against rapid changes, and threshold again to 0.5
take difference between adjacent frames again to get the changes between ticking and paused

Also, the timer may not be visible at some parts of the video, add the following step(s) before applying the gaussian smooth

for each timer region found during the scanning process, crop all scanned frames to the timer region
group these images by similarity (squared difference), and select only the ~10-15 groups with the highest frequency, these are most likely the actual digits
average all images of each group to get an idealized sample for each digit
if the timer does is not rectangular (eg. italic text), save the original mask for the timer region and compare only these pixels

Now, after taking the difference between adjacent frames, set all values to 0 (stopped) where the image does not match any of these digits (again, squared difference).

(left: all images from intervals classified as containing the timer extracted from the timer area; right: images grouped by sum of squared difference)

Considerations

There are a few things to consider when applying these methods.

Acquiring frames. Sounds simple, but you don't want to dump them all uncompressed to your HDD. I'm using ffmpeg with fast seeking -ss time -i input.mp4, which seeks to the nearest key frame and decodes frames until it reaches the requested time. Except it won't work some times, especially for bad video files that are common when you are saving a network stream. After transcoding the video file (with ffmpeg), most files that didn't work before were recognized perfectly. If you scale with ffmpeg, encode to 320x180, set a low key frame interval and set -movflags faststart, reading frames from the resulting file is accurate and fast. And at 320x180, it doesn't use much space either. Avoid downsampling the frame rate, as that might destroy your time stamps.
Choosing the right threshold values, time intervals and other parameters. This is what took by far most of my time, especially since testing takes some time.
In particular, it helps to choose a time step such that each digit appears. Choosing 5 seconds for example is a bad choice, as it will give you only 0's and 5's. Choosing a time step lower than 1 second won't work well when trying to analyze how much each pixel changes, as the consecutive digits will be the same.
Related to this, each interval needs to contain enough frames. If the period is 10s and you choose 60 frames for the interval length and a time step of 1.2 seconds, the length after which the digit repeats is $lcm(1.2s,10s)=60s$, ie. 50 frames out of 60. And that is not enough. I found that, in order to eliminate unwanted noise, the number of frames should be least 3-6 times as high as the period length in frames. So a time step of, say 3s, works quite well, as lcm(3,10)=30, which is 10 frames, so you get 6 periods.

Does it work, though?

I made a simple ruby script that implements the method described above.

It's probably full of bugs and the code is a mess, but it seems to work. I tested it on 200 videos of different length (5 minutes to 6 hours, most around 30 minutes). They all contained the same kind of timer (color, font), but at possibly different sizes and positions. After transcoding those video files where ffmpeg didn't manage to seek accurately, the script got the correct result for pretty much all of them, it got only 4 wrong (detected stopped when it had been ticking).

It's possible I might have fine-tuned the parameters for that specific testing set, but at least we can say it does work rather well if you can make some general assumptions about the timer. Also, I tested a few other videos with different font sizes/colors, and they worked as well. I also tested it on a video with two timers (one that kept ticking the entire video), and the script recognized them as well.

As far is accuracy is concerned, it usually determines the starting/ending point at $\pm 1...3$ seconds, (because of some smoothing and averaging to eliminate noise), which is more than enough for me. I guess you could get an even better estimate if, after determining there is a starting/ending point at a certain frame interval, you check exactly when the digit starts/stops changing.

Sample output:

Transcoding video (110.89 min) file before analysis...

ffmpeg -y -loglevel quiet -i "dump_000000038.mp4" -c:v libx264 -an -vf scale=320:180 -preset ultrafast -sws_flags bilinear -g 20 -movflags faststart -crf 12 -f mp4 ./transcode.mp4 < /dev/null
Transcoding took 423.545512367s, 16x real-time.

Analyzing <transcode.mp4>, duration 1:50:53.684...

Scanning for timers...

Trying timestep 11.0, method simple...
Analyzing 0:00:00.000-0:11:00.000...
 - got 0 connected area(s)
Found 0 candidates.

Analyzing 0:11:00.000-0:22:00.000...
 - got 1 connected area(s)
Found 1 candidates.

Analyzing 0:22:00.000-0:33:00.000...
 - got 1 connected area(s)
Found 1 candidates.

Analyzing 0:33:00.000-0:44:00.000...
 - got 0 connected area(s)
Found 0 candidates.

Analyzing 0:44:00.000-0:55:00.000...
 - got 1 connected area(s)
Found 1 candidates.

Analyzing 0:55:00.000-1:06:00.000...
 - got 1 connected area(s)
Found 1 candidates.

Analyzing 1:06:00.000-1:17:00.000...
 - got 1 connected area(s)
Found 1 candidates.

Analyzing 1:17:00.000-1:28:00.000...
 - got 1 connected area(s)
Found 1 candidates.

Analyzing 1:28:00.000-1:39:00.000...
 - got 1 connected area(s)
Found 1 candidates.

Analyzing 1:39:00.000-1:50:00.000...
 - got 0 connected area(s)
Found 0 candidates.

Extracting timer digit images...
Grouping timer digit images...

Timer(s) found:
Timer 1: (176,150), size 5x10

Analyzing timer(s) for event(s)...

Trying coarse timestep 11.0, fine 1.75...
Scanning timers for activity...
Timer 0: Found 4 event(s).
Checking at the beginning of the video at a higher resolution for the timer 0...
Checking at the end of the video at higher a resolution for the timer 0...
Refining estimate for timer 0, event 0...
Refining estimate for timer 0, event 1...
Refining estimate for timer 0, event 2...
Refining estimate for timer 0, event 3...
Reading frame 50/50...

This is what they're doing:
Timer 1:
  - starts off paused.
  - starts at 0:04:53.000
  - stops at 0:37:32.000
  - starts at 0:38:57.000
  - stops at 1:45:26.000
  - ends up paused.

Scanning took 25s, 16265 fps.


2D FFT how is it done? I have a texture that I want to perform the 2DFFT on and I am trying to understand the "Row/Column" or "Column/Row" idea but I am unsure if I have understood it correctly. Hoping some one can explain it a bit better from my current understanding because it seems my current understanding of it requires a lot of FFTs which does not seem right to me.
If I have a 256 by 256 texture and I perform a 2DFFT on the image. Lets say I chose columns/row order so I do columns first, so does that mean I first perform:
256 1DFFT's... one for each column "image.x", containing 256 samples:
1DFFT for x = 0 -> Samples: image[0,0], image[0,1] ... image[0,255]
1DFFT for x = 1 -> Samples: image[1,0], image[1,1] ... image[1,255]
...
1DFFT for x = 255 -> Samples: image[255,0], image[255,1] ... image[255,255]

Then I have to do another 256 1DFFT's along the rows "image.y", containing 256 samples:
1DFFT for y = 0 -> Samples: image[0,0], image[1,0] ... image[255,0]
1DFFT for y = 1 -> Samples: image[0,1], image[1,1] ... image[255,1]
...
1DFFT for y = 255 -> Samples: image[0,255], image[1,255] ... image[255,255]

So in total I do 256 + 256 = 512 1DFFTs ? Is that correct?

Yes, you are right. It is a lot of FFTs, but remember it is a lot of samples as well.
If you have a FFT algorithm with complexity $O(n \cdot log(n))$, then for the 2D scenario you will have an algorithm with complexity $O(M \cdot (N \cdot log(N)) + N\cdot (M \cdot log(M)))$
but if you express this in terms of $M \cdot N$, what you get is.
$O(M \cdot N \cdot (log(N) + log(M)))$ and finally with the log sum identity you get $O(M \cdot N \cdot log( M \cdot N ))$
The complexity in terms of number of samples is the same


equations to calculate coefficients of recursive IIR Gaussian blur filter I am reading the implementation of IIR Gaussian blur filter from this link.

In the source code, the coefficients are generated by the following code

/* Calcualte Gaussian Blur Filter Coefficiens
 *  alpha -> smooting gradient depends on sigma
 *  k = ((1-e^-alpha)^2)/(1+2*alpha*e^-alpha - e^-2alpha)
 *  a0 = k; a1 = k*(alpha-1)*e^-alpha; a2 = k*(alpha+1)*e^-alpha;
 *  a3 = -k*e^(-2*alpha)
 *  b1 = -2*e^-alpha; b2 = e^(-2*alpha)
 */
 void calGaussianCoeff( float sigma,  float *a0, float *a1, float *a2,
   float *a3, float *b1, float *b2, float *cprev, float *cnext)
{
  float alpha, lamma,  k;
  // defensive check
  if (sigma < 0.5f)
     sigma = 0.5f;

  alpha = (float) exp((0.726)*(0.726)) / sigma;
  lamma = (float)exp(-alpha);
  *b2 = (float)exp(-2*alpha);
  k = (1-lamma)*(1-lamma)/(1+2*alpha*lamma- (*b2));
  *a0 = k;
  *a1 = k*(alpha-1)*lamma;
  *a2 = k*(alpha+1)*lamma;
  *a3 = -k* (*b2);
  *b1 = -2*lamma;
  *cprev = (*a0 + *a1)/(1+ *b1 + *b2);
  *cnext = (*a2 + *a3)/(1+ *b1 + *b2);
}

Does anyone know where the equation is from?

@Peter K.'s answer is right to the point. You can also have most proposed implementations side-by-side, with comparisons,  in Pascal Getreuer, A Survey of Gaussian Convolution Algorithms, Image Processing On Line, 2013, with alternative source code.


Why 0.5 in value of heaviside/unitstep? When i want to see heaviside or unit step function in matlab, as shown in attached figure. i get first value of 0.5,why 0.5? why not 1,as unit step is thought to be ?

Well i don't have the Symbolic Math Toolbox.  but i think i know why, and it's not a really good reason.  usually we electrical engineers mean, when we say "the unit step" function, we mean what math and physics people mean by the "Heaviside step function but, there is a convention issue with what this function is with zero for its argument.

the root definition of the Heaviside unit step is

$$ H(x) \triangleq \int_{-\infty}^x { \delta(t)} \, \mathrm{d}t $$

where $\delta(t)$ is the Dirac unit impulse function.  But this does not define $H(x)$ at $x=0$.  It's infinitesimally thin, but it could be anything finite at zero.

in continuous time, it is usually most convenient to define the Heaviside step function in terms of the sign function which is:

$$\operatorname{sgn}(x) \triangleq \begin{cases} -1, & x < 0 \\ 0, & x = 0 \\ 1, & x > 0 \end{cases}$$

then the Heaviside comes out as

$$\begin{align}
 H(x) &= \tfrac12 \big(1 + \operatorname{sgn}(x) \big) \\
\\
&= \begin{cases} 0, & x < 0 \\ \tfrac12, & x = 0 \\ 1, & x > 0 \end{cases} \\
\end{align}$$

because this, $H(x) + H(-x) = 1$, is always true, they are complementary.

the problem is extending these definitions and desired behavior to discrete $x$.  then we like to define the discrete impulse as the Kronecker delta which DSPer's express more simply as

$$ \delta[n] \triangleq \begin{cases} 0, & n \ne 0 \\ 1, & n = 0 \end{cases}  $$

and the discrete unit step as

$$ H[n] \triangleq \sum_\limits{m=-\infty}^{n} \delta[m] $$

but that makes this discrete unit step come out to be:

$$ H[n] = \begin{cases} 0, & n < 0 \\ 1, & n \ge 0 \end{cases}  $$

that is not quite the same as

$$\begin{align}
 H[n] &= \tfrac12 \big(1 + \operatorname{sgn}(n) \big) \\
\\
 &= \begin{cases} 0, & n < 0 \\ \tfrac12, & n = 0 \\ 1, & n > 0 \end{cases} \\
\end{align}$$

now here is the stupid problem:  MATLAB is not a continuous-time tool.  it is a discrete-time tool.

so what they are doing with the discrete Heaviside unit step function appears to me to be another one of MATLAB's mistaken conventions of definitions.  they fucked up the index origin convention which, among other problems fucked up the fft() convention putting the DC term into X(1) rather than have the index represent frequency.  they fucked up the order of polynomial coefficients in polyval() and polyfit().  and they fucked up the definition of the discrete Heaviside unit step function.

EDIT: actually, they defined the Heaviside function correctly, because it has a continuous argument.  but they (and we) need to define a discrete function for the unit step that does not have $\tfrac12$ in it.