How to set sampling frequency of an acceleration data set I have a set of accelerometer readings in $X, Y, Z$ axes obtained from an android based smart-phone. The data in $(X, Y, Z)$ was recorded at different time stamps and there is no uniform time period for recording the data.

My data set looks as  timestamp $X, Y, Z$. First, I did a filtering (using a low-pass) on this data and now want to perform $\textrm{FFT}$ on a timewindow of 1 min or may be a window containing 250 samples (not very sure about the window length). I am using this FastFourierTransform class of apache commons in Java.

I am getting $\textrm{FFT}$ magnitude but I was wondering how do I know the corresponding frequency of each $\textrm{FFT}$ mag. I know corresponding frequency would be $(n\cdot F_s)/N$ where $n$ is the bin number or data point index, $F_s$ is sampling frequency and $N$ is number of input data points over the window. Now my question is how do I know $F_s$ for a given set of data say for example input data in an array as $\begin{bmatrix}1&2&3&4&5&6& 7\end{bmatrix}$ over a window of size $N_t$ millisec where $N_t=(\textrm{lastTimeStamp}-\textrm{firstTimeStamp})$

Is it $\displaystyle \frac{7}{(N_t\cdot 10^{-2})} \textrm{Hz}$?

This also immediately follows by doing so if this agrees with nyquist sampling frequency? and this says $F_s\ge2F_{max}$ in the signal and I don't know the maximum frequency over the time window $N_t$.

If you don't have a uniform time period (or nearly so) between samples, then you don't have a sample rate, and the FFT result may well be garbage.

If you don't know ahead of the sampling the  Nyquist frequency for the signal, then the sampling might be already contaminated by aliasing.

If the samples are close enough together relative to the highest frequency in the signal's spectrum, then you may be able to interpolate a new set of evenly spaced samples, which will then have a known sample rate (based on the distance of the new even spacing).


what is continuous wavelet (cwt) ,wavelet packet (wpt) and stockwell (S-T) computational complexity? What is computational complexity of continuous wavelet transform
cwt ?
What is computational complexity of wavelet packet transform wpt ?
What is computational complexity of Stockwell transform transform wpt ?
for the cwt : imagine I have a signal with length of $N$ and I
scale vector is 0:1:s
what will happen if I discretize scale vector to 0:dj:s for better accuracy?
for the WPT : Imagine signal size is $N$ and I am decomposing it in
$L$ levels.

note:

for Continuous case I guess its $O(S\times N^{2})$ and for 0:dj:S, It must be $O(\frac{S}{dj}\times N^{2})$ , Am I right?
for WPT transform it must be $O(2^{L}\times NLogN )$ . Is that right?
For S transfrom I think it is $O(N^{2}Log N)$ . Is it true?
I have calculated the running time of Stockwell and Continuous wavelet transforms. but the time results are strange.
N=3000;            $O(Stockwell)=1.039\times 10^{8}$  it takes 2.5 sec
N=3000, scale=300; $O(cwt)=2.7\times 10^{9}$ it takes 0.58 sec
why this is happening here ? [P.S I have run them 100 times and got the time average]

For a dsicretized signal with N samples. In terms of complexity for the DWT it is O(N) and for WPT its O(Nlog(N)) REFER. For each scale one calculates a convolution and not a single correlation value and the scale values need to be integral values. They can be constrained by using the iterative structure of quadrature mirror LP and HP filters, which is basically a explicit way of coding the scaling function. This is basically because one needs a finite number of scalings of the wavelet to represent the signal (wavelet is overcomplete frame). REFER2. For the wavelet packet tree the extra cost is due to the encoding of the residue/detail coefficients.

EDIT:
You are right the discretization of the scale function in the CWT parameter decides the number of levels of decomposition. Matlab implements the CWT and describes that the DFT can be used in calculating the CWT by zeropadding the wavelet and performing circular convolutions.

Here Δt is the sampling interval (period) which basically discretizes the signal and the wavelet. So i guess you are right about the complexities for both the CWT and the WPT versions! Sorry about the mix up!


Finding radius r of the overlappable sphere(s) in 3D image My current problem:

I have an input 3D binary image (a 3D matrix that has only 0 and 1) that consists of random numbers of sphere with radius r.
We do not know how many spheres are there in the image.
All spheres have the same radius r, but we do not know radius r.
Spheres are all over the place in the image and can overlap each other.
example image is given below.

My requirement:

what is the radius r?

Currently, I simply just flatten the image to get rid of the z axis and perform edge detection and I am trying Hough Transform using: http://rsbweb.nih.gov/ij/plugins/hough-circles.html

However, with Hough Transform, I see that the variables minimum radius, maximum radius, and number of circles have to be specified. I have tried a few attempts below:

Given the correct parameters, Hough Transform can detect the circles just fine. But in the real application, I do not know how many spheres are there, and making the program attempting to guess minimum and maximum radius seems not feasible. Are there other ways to accomplish this?

Cross-link: https://math.stackexchange.com/questions/118815/finding-radius-r-of-the-overlappable-spheres-in-3d-image

A simpler solution and much more computationally efficient when compared to Hough Transform is to use the distance transform:

Find the surface of your spheres (i.e. the pixels that have value 1 and have at least one neighboring 0 pixel);
Compute the distance transform with respect to the spheres surface, but constrain the computation only to pixels that are internal to the spheres. The output will be a distance map;
The radius will be exactly the maximum value in your distance map.

Another advantage of this solution when compared to Hough transform is that it provides a much more precise value for the radius.


Why do we read that Hilbert transform can be used for envelope detection? We can often read that Hilbert transform is useful for envelope detection (e.g.
Hilbert transform to compute signal envelope?)

I have done some tests with various soundfiles, and

x[n] -> absolue value -> 1-pole low pass filter -> envelope

or

x[n] -> Hilbert transform -> absolute value -> 1-pole low pass filter -> envelope

gives the same kind of result, it's not better with a Hilbert transform (top : envelope computed with hilbert, bottom : envelope computed without hilbert) :

Moreover, I know that computing the Hilbert transform is very time-consuming (big FIR filters involved).

So it is really a good method for envelope detection ?

Using a Hilbert transform with a sufficiently narrow-band signal, the low pass filter may not be needed even for continuous envelope estimation.


Gaussian, Rayleigh, and Exponential RVs I am doing some simulations on the performance of OFDM systems in multi-path channels and I am confused about something. It is explained below:

Assume that we have a multi-path channel with a number of taps equal to $L$. Each tap has a power of $\sigma$, hence, the total power in the multi-path channel in time domain is $\sigma\:L$. If the generated channel in time domain of length $L$ is denoted as $h$, then the channel in frequency domain of length $N$ is $H = \rm{FFT}(h,N)$, where $N$ is the number of OFDM subcarriers.

If we have $X \sim N(0,\frac{\sigma^2}{2})$ and $Y \sim N(0,\frac{\sigma^2}{2})$, then $h = \sqrt{X^2 + Y^2} \sim \rm{Rayleigh}(\frac{\sigma}{\sqrt(2)})$.
Additionally, if channel in frequency domain $H\sim \rm{Rayleigh}(\frac{1}{\sqrt{2\lambda}})$, then $|H|^2 \sim e^{\lambda}$? In this line correct? Is $\lambda = \frac{1}{\sigma^2}$?

I am confused about the transition from $h$ to $|H|^2$ and the relation between their means and variances.

I appreciate your help in advance.

So you have a complex impulse response, and the real and imaginary parts of each tap are modeled as i.i.d. zero mean Gaussian variables with variance $\sigma^2/2$ (correct me if I'm wrong). The DFT of the impulse response is a weighted sum of the filter taps:

$$H[k]=\sum_{n=0}^{N-1}h[n]e^{-j2\pi nk/N}\tag{1}$$

where I assume $N\ge L$ (where $L$ is the number of taps), and in Eq. (1) $h[n]=0$ for $n\ge L$. From (1) it follows that $H[k]$ is also a complex Gaussian RV with zero mean and with variance $L\sigma^2$ (the variance of the real and imaginary part is $L\sigma^2/2$, respectively). From this it follows that the magnitude $|H[k]|$ has a Rayleigh distribution:

$$f_{|H|}(x)=\frac{2x}{L\sigma^2}e^{-x^2/(L\sigma^2)}\tag{2}$$

And, accordingly, the squared magnitude $|H[k]|^2$ has an exponential distribution:

$$f_{|H|^2}(x)=\lambda e^{-\lambda x},\quad x\ge 0\tag{3}$$

With $\lambda=1/(L\sigma^2)$, i.e. with mean $1/\lambda=L\sigma^2$, and variance $1/\lambda^2=L^2\sigma^4$.


Derive Bode Diagram from difference equation I have designed a 2nd order low pass butterworth filter for the following:

Sampling Frequency: 10kHz. Cut-off Frequency: 1kHz

After doing a bit of maths I worked out the difference equation:

w=tan(1000*pi / 10000) = 0.3249
H(Z) = 0.06747(z + 1)^2 / z^2 - 1.143z +0.413

Poles / Zeroes:

a1=1.143, a2=-0.413, b0=0.06747, b1=0.13494, b2=0.06747

Difference Equation:

$$y[n] = 0.06747x[n] + 0.13494x[n-1] + 0.06747x[n-2] + 1.143y[n-1] - 0.413y[n-2]$$

Difference Equation in C#

Y[n] = a1 * Y[n - 1] + a2 * Y[n - 2] + b0 * Signal[n] + b1 * Signal[n - 1] + b2 * Signal[n - 2];

The problem I'm having is trying to plot the frequency response using a Bode Diagram (i.e the low pass filter curve).

I'm struggling to understand how the difference equation can be used to plot the bode diagram. How can I do it? I understand that it must be dB magnitude.

To get Y[n] I need to input a signal into the Signal[n] array. I should be able to produce a graph that looks something like this:

You can evaluate the frequency response by evaluating $H(z)$ in the unit circle: $H(e^{j\theta})$, where $\theta$ is the discrete-time frequency.

If you prefer a frequency in Hz, then use $\theta = 2\pi/f_s$ (i.e. the sampling frequency corresponds to $\theta = 2\pi$).

The above expression will give you the frequency response as a complex number, signifying both amplitude and phase response.

If you don't have access to complex number algebra on your platform, use either GNU Octave or Matlab to plot the frequency response (using the function freqz, which takes the a and b coefficients you have already).

Note that Bode diagram refers to an asymptotic approximation of the frequency response, valid only in continuous time systems.


Discrete time signals connected in series Two discrete-time systems are connected in series. Their input-output difference equations are:

  \begin{align}
w[k]&=0.25x[k-2]\\
y[k]&=0.5w[k-1]+0.5w[k-2].
\end{align}

  Determine the overall input-output difference equation relating $y[k]$ to $x[k]$
  Determine the impulse response of the overall system

Julie.

You simply need to put $w[k]$ into $y[k]$:

$$w[k]=0.25x[k-2] \Rightarrow \begin{align}w[k-1]&=0.25x[k-3]\\
 w[k-2]&=0.25x[k-4]
\end{align}$$

\begin{align}
y[k] &= 0.5w[k-1] + 0.5w[k-2]\\
&=0.5\cdot0.25x[k-3] + 0.5\cdot0.25x[k-4]\\
&= 0.125(x[k-3]+x[k-4])
\end{align}

You can directly read the impulse response from it.


Can I use independent component analysis if I have two discrete populations of signal sources mixed together? I understand (mostly) how independent component analysis (ICA) works on a set of signals from one population, but I am failing to make it work if my observations (X matrix) includes signals from two different populations (having different means) and I am wondering if it is an inherent limitation of ICA or if I can resolve this.  My signals are different than the common type being analyzed in that my source vectors are very short (e.g. 3 values long), but I have many (e.g. 1000's) of observations.  Specifically, I am measuring fluorescence in 3 colors where the broad fluorescence signals can "spillover" into other detectors. I have 3 detectors and using 3 different fluorophores on particles. One could think of this as a very poor resolution spectroscopy.  Any fluorescent particle could have an arbitrary amount of any of the 3 different fluorophors.  However, I have a mixed set of particles which tend to have quite distinct concentrations of fluorophores. For example, one set may generally have lots of fluorophore #1 and little fluorophore #2, while the other set has little of #1 and lots of #2.

Basically, I want to deconvolve the spillover effect to estimate the actual amount of each fluorophore on each particle, rather than having a fraction of signal from one fluorophore add to the signal of another.  It seemed like this would be possible for ICA, but after some significant failures (the matrix transform seems to prioritize separating the populations rather than rotating to optimize signal independence), I am wondering if ICA isn't the right solution or if I need to pre-process my data in some other manner to address this.

The graphs show my synthetic data used to demonstrate the problem.  Starting with "true" sources (panel A) consisting of a mixture of 2 populations, I created a "true" mixing (A) matrix and calculated the observation (X) matrix (panel B).  FastICA estimates the S matrix (shown in panel C) and instead of finding my true sources, it seems to me that it rotates the data to minimize the covariance between the 2 populations.

Looking for any suggestions or insight.

Conceptually, for ICA to work, it requires N observations of the SAME phenomenon and it will attempt to decompose these signals (the mixtures) to < N components.

From your description each detector is counting a different energy particle. That is, the "Mixtures" are actually already three different signals. Consider the example where you use the recordings of 4 different microphones to try and distinguish between 2 different (mutually independent) sound sources. In this case each microphone is recording a (different) sum of BOTH sound sources. In your case, because of the phenomenon you are observing, each "microphone" (detector) is NOT recording a (different) sum of BOTH sources, but rather it records each source separately (most of the times).

As you describe, sometimes, the SAME event might be recorded by more than one detectors (spill-over). In this case, if you insist with ICA then it would attempt to decompose THAT sum.

Therefore, perhaps you should consider something simpler, like obtaining the cross correlation between the events recorded by each one of the three different detectors as an indicator of the amount of spill-over.

But in this case, instead of examining 1000 "triplets" of observations, you create three different time series. One time series per detector describing the events (in time) of each detector.

If your detectors and aparatus are fast enough, this signal will look like a series of exponentially decaying pulses. If not, your signal will look like a train of pulses convolved with the impulse response of the system (a smoothed version). In either case, if the same event has been detected by more than one detectors then the output of these detectors will simultaneously rise in level (and increase the cross - correlation or some other form of correlation metric). It will be easier to detect that rise in the first case (fast system).

Therefore, the problem now becomes one where you have to decide which detector to assign that event to for it to be counted correctly. There is a number of ways that you can do this (not using ICA) but you will have to know the spectral response of your detector (Pulse amplitude VS energy or Voltage output VS energy). You would then be able to create a weighted sum of the counted particles by looking up the "weight" for each pulse on your detector's response. Or you could apply a "window" to the pulses. That is, consider a sharp threshold around the peak response of the detector and only count events that generate such (usually high) pulses. Considering a short region around the peak performance of the detector will reject the (lower) pulses created by spill-over.

I hope this helps. I hope to have understood your configuration correctly.


How to derive the "well-known" solution to Unconstrained Array Gain? Can someone point me to a webpage or other resource that shows how to analytically solve the beamformer Unconstrained Array Gain expression in Henry Cox's 1987 IEEE paper "Robust Adaptive Beamforming"?
$$ \max_{\mathbf{w}} \frac{|\mathbf{w}^H\mathbf{d}|^2}{\mathbf{w}^H\mathbf{Q}\mathbf{w}} $$
Cox says:

The well-known solution is $\mathbf{w} = \alpha\mathbf{Q}^{-1}\mathbf{d}$

I'd just like to better understand this by learning how to derive this myself.

You can solve such a problem using the method of Lagrange multipliers. First note that maximizing the expression in your question is equivalent to minimizing the inverse function:
$$\min_{\mathbf{w}}\frac{\mathbf{w}^H\mathbf{Q}\mathbf{w}}{|\mathbf{w}^H\mathbf{d}|^2}\tag{1}$$
Next note that the solution of $(1)$ is invariant to scaling of $\mathbf{w}$, i.e., replacing $\mathbf{w}$ by $c\cdot\mathbf{w}$ in $(1)$ with an arbitrary scalar constant $c$ will not change the value of the function. So we may as well use a scaling such that $\mathbf{w}^H\mathbf{d}=1$ is satisfied. This scaling corresponds to a unity response for the desired signal. With this constraint, problem $(1)$ can be reformulated as
$$\min_{\mathbf{w}}\mathbf{w}^H\mathbf{Q}\mathbf{w}\qquad\textrm{s.t.}\qquad \mathbf{w}^H\mathbf{d}=1\tag{2}$$
We can solve $(2)$ using the method of Lagrange multipliers by minimizing
$$\mathbf{w}^H\mathbf{Q}\mathbf{w}-\lambda(\mathbf{w}^H\mathbf{d}-1)\tag{3}$$
Formally taking the derivative of $(3)$ with respect to $\mathbf{w}^H$ and setting it to zero gives
$$\mathbf{w}=\lambda\mathbf{Q}^{-1}\mathbf{d}\tag{4}$$
The constraint in $(2)$ is satisfied for
$$\lambda=\frac{1}{\mathbf{d}^H\mathbf{Q}^{-1}\mathbf{d}}\tag{5}$$
From $(4)$ and $(5)$ we finally obtain
$$\mathbf{w}=\frac{\mathbf{Q}^{-1}\mathbf{d}}{\mathbf{d}^H\mathbf{Q}^{-1}\mathbf{d}}\tag{6}$$
Note that the scaling in $(6)$ is optional and the general solution is given by $(4)$.


Unit Impulse : Bounded or Unbounded? As the title suggests, I have a confusion here. In a Systems text I am going through, it mentions of the Unit Impulse as an unbounded signal. Yes, the unit impulse "height" is unbounded, but the "strength" of the signal as the way it is mentioned in many texts, is finite, given by the area of the unit impulse. By that line of thinking it is giving me a bounded signal. What am I missing here? How is the unit impulse interpreted/treated in terms of being an input or an output to a system? Thanks!

Bounded or unbounded in that context always refers to amplitude. E.g., BIBO (bounded-input-bounded-output) stability means that a bounded input signal $x(t)$ satisfying $|x(t)|\le A$ always results in a bounded output $y(t)$ satisfying $|y(t)|\le B$ with some finite positive constants $A$ and $B$.

So the Dirac delta impulse is not well-suited for testing stability of a system. E.g., an ideal integrator's response to an impulse at the input is a unit step, which is of course bounded. Yet, the ideal integrator is unstable because its response to a step is unbounded.

In sum, a Dirac delta impulse is unbounded.


How to measure correlation between 1 bit samples and theoretical original waveform? I'm investigating quantization error.
I have an analogue waveform that looks like this and is of theoretically infinite resolution:-

I've sampled it as (8 bit oscilloscope readings & 0b1) to produce a file of single ones or zeros. The samples are partially random as per quantization theory, but I want to check the correlation with the analogue original. Theory suggests that there should be some correlation.
How can I do that? Clearly I can't compare the 1 bit samples with the 8 bit samples they came from.
Put another way, how can I show that my quantization is entirely uncorrelated to the underlying waveform?

Edit & maybe a solution:-
St is the theoretical and perfect(!) signal.
S1 is my sample set as (scope reading) & 0b1.
S2 is another independent sample set as (scope reading) & 0b1.
It occurs to me that I can't possibly compare S1 with St, as quantisation error will be present in both sets. But, if S1 ~ St, then by inference S2 ~ St. Hence I can test the hypothesis that S1 ~ S2 since they're both derived from common signal St. If there is a correlation, then S1 must be correlated to St too.
So I can just use normal correlation analysis (or Hamming distance, Jaccard index) between S1 and S2 to measure any. Yes/no?

I've sampled it at a 1 bit resolution (as 8 bit oscilloscope readings & 0b1)

No. Taking the LSB is NOT 1 bit quantization. A one bit quantizers are not trivial to design. For a zero mean signal that's reasonably symmetrically distributed over positive and negative values, you can just take the sign bit. For signals that are asymmetrically distributed around their mean, designing a decent 1 bit quantizer takes a bit of work.
In your example you should quantize so that you get a "1" of the signal is above the mean and and "0" if it's below the mean. If you do this, you will see that both quantization noise and quantized signal are highly correlated with the original signal. Both will be periodic with the same period.
If you want to study quantization noise, I recommend starting with higher bit levels, i.e. where the amplitude of the signal is much larger than the quantization steps. This makes the choice of quantizer very simple and in most case you can assume that the quantization noise is uniformly distributed. Otherwise, you need to tackle the subject of "optimum quantizer" which requires some non-trivial math.


How to Average Complex Responses (and Justification)? I am developing software that calculates the response of a system by comparing the FFT of input and output signals.  The input and output signals are divided into windows and, for each window, the signals are median-subtracted and multiplied by a Hann function.  The instrument response for that window is then the ratio of the FFTs of the processed data.

I believe the above is standard procedure, although I may be describing it poorly.  My problem comes in how to combine the responses from the multiple windows.

As far as I can see, the correct approach is to average the complex values, across all windows.  The amplitude and phase response are then the amplitude and phase of the average, complex value at each frequency:

av_response = sum_windows(response) / n
av_amplitude = sqrt(real(av_response)**2 + imag(av_response)**2)
av_phase = atan2(imag(av_response), real(av_response))

with implicit loops over frequency bins.

But I have been asked to change this to calculate amplitude and phase in each window first, and then average the amplitudes and phases across all windows:

amplitude = sqrt(real(response)**2 + imag(response)**2)
av_amplitude = sum_windows(amplitude) / n
phase = atan2(imag(response), real(response))
av_phase = sum_windows(phase) / n

I have argued that this is incorrect because averaging angles is "just wrong" - the average of 0 and 360 degrees is 180, for example, but the people I am working with responded by saying "OK, we will only display amplitude".

So my questions are:

Am I correct in thinking that the second approach is generally incorrect
for amplitudes too?
If so, are there any exceptions that may be relevant, and which may explain
why the people I am working with prefer the second method?  For example,
it looks like the two approaches will agree as the noise becomes small, so perhaps this
is an accepted approximation for low noise?
If the second approach is incorrect, are there any convincing, authoritative references that I can use
to show this?
If the second approach is incorrect, are there any good, easy to understand examples that show this for amplitude
(as the average of 0 and 360 degrees does for phase)?
Alternatively, if I am incorrect, what would be a good book for me to educate
myself better?

I have tried to argue that the average of -1 1 1 -1 1 -1 -1 should be zero rather than 1, but that was unconvincing.  And while I think I could, with time, construct an argument based on max likelihood estimation given a particular noise model, it is not the kind of reasoning that the people I am working with will listen to.  So, if I am not wrong, I need either a powerful argument from authority or an "obvious" demonstration.

[I tried to add more tags, but can't find relevant ones and can't define new ones as a new user - sorry]

Transfer function estimation is usually implemented slightly differently than the method you describe.

Your method computes

$$\left\langle \frac{\mathcal{F}[y]}{\mathcal{F}[x]} \right\rangle$$

where $\langle$angle brackets$\rangle$ represent averages taken over data segments, and a windowing function is applied to each data segment before taking the Fourier transform ($\mathcal{F}$).

A more typical implementation will compute the cross spectral density of x and y divided by the power spectral density of x:

$$\frac{\langle \mathcal{F}[y] \cdot \mathcal{F}[x]^* \rangle}{\langle|\mathcal{F}[x]|^2\rangle} = \frac{\langle \mathcal{F}[y] \cdot \mathcal{F}[x]^* \rangle}{\langle\mathcal{F}[x]\cdot\mathcal{F}[x]^*\rangle}$$

Where $\cdot$ represents a pointwise product, and $*$ the complex conjugate.

I believe this is to reduce the effect of data segments where bins of $\mathcal{F}[x]$ are excessively small.

Incoherent estimation

Your employer has suggested that you estimate the transfer function using

$$\frac{|\langle \mathcal{F}[y]|\rangle}{\langle |\mathcal{F}[x]|\rangle}$$

This will work, but has two big disadvantages:

You don't get any phase information.
If your measurements of the input $x$ and output $y$ have any additional noise, then the transfer function estimation will not be correct.

Your method and the method I described circumvent these problems by using coherent averaging.

References

The general idea of using overlapped, averaged segments to compute power spectral densities is known as Welch's method.  I believe the extension to using this to estimate transfer functions is also often known as Welch's method, although I'm not sure if it is mentioned in Welch's paper.  Looking up Welch's paper might be a valuable resource.  A useful monograph on the subject is Bendat and Piersol's book, Random Data: Analysis and Measurement Procedures.

Validation

To validate your software, I suggest applying several test cases, where you generate Gaussian white noise and feed it through a digital filter with a known transfer function.  Feed the inputs and outputs into your transfer function estimation routine and verify that the estimate converges to the known value of the transfer function.


How to sketch the following discrete-time signal? How to sketch the following discrete-time signal:
\begin{gather*}
x[(n-1)^2]
\end{gather*}
for the signal given below:

Any help would be much appreciated.

I give you some hints and then you can solve this homework.

Your $x[n]$ has only $8$ nonzero values. Figuring out what happens to them (in an exhaustive way) is not difficult.
Consider $(n-1)^2$ and change $n$ over $(-\infty, \infty)$. You can see that it is always non-negative. It means that we never can have $(n-1)^2=-1$. Hence the samples of $x[n]$ for $n=-3,-2,-1$ are never used. We can only have samples of $x[n]$ at $n=0,1,2,3$.
Solve $(n-1)^2=0$, $(n-1)^2=1$,$(n-1)^2=2$, $(n-1)^2=3$. Only accept those results that lead to an integer $n$. Then you can see at what $n$ in $x[(n-1)^2]$ we have what value of $x[n]$. For instance, $(n-1)^2=\color{red}1\Rightarrow n=0,2$. It means in $x[(n-1)^2]$ the samples at $n=0$ and $n=2$ are equal to $x[\color{red}1]=1$.
Except at the points mentioned above, the transformed signal would be zero.


Plotting phase delay from first principles giving different result from phasedelay() function. Why? I have a length 10 Finite Impulse response as shown below, and I am trying to find the phase delay of this FIR.

$$h[n]=\begin{cases}0.1,& 0\le n \lt 10\\
0,&\textrm{otherwise}\end{cases}$$

I tried plotting from first principles as well as using the phasedelay() function, however I am getting different results. The results are varying greatly. Any suggestions on what I am forgetting to do in my approach from first principles?

The problem is the definition of the phase. The command angle() computes the phase $\phi(\omega)$ according to

$$H(e^{j\omega})=|H(e^{j\omega})|e^{j\phi(\omega)}\tag{1}$$

where $H(e^{j\omega})$ is the complex frequency response. At frequencies $\omega$ where the frequency response has zeros, the phase jumps by $\pi$. This is shown in the two left figures of the plot below. You can't compute the phase delay from $\phi(\omega)$ because of the jumps. Note that in your computation, the phase delay is correct until the first phase jump due to the first zero of the magnitude response.

The phase delay needs to be computed from a continuous phase $\varphi(\omega)$, which is implicitly defined by

$$H(e^{j\omega})=A(\omega)e^{j\varphi(\omega)}\tag{2}$$

where $A(\omega)$ is a real-valued but bi-polar smooth function. This continuous phase is shown in the top right figure below, and it is obviously linear. From this continuous phase, the phase delay can be computed by

$$\tau_{ph}=-\frac{\varphi(\omega)}{\omega}\tag{3}$$

and for the given linear-phase FIR filter the phase delay is constant:

$$\tau_{ph}=\frac{N-1}{2}\tag{4}$$

where $N$ is the filter length.

Note that the continuous phase $\varphi(\omega)$ cannot be obtained from $\phi(\omega)$ using the Matlab function unwrap(), because the latter just removes jumps by (multiples of) $2\pi$ due to phase ambiguity. Here we have to deal with actual phase jumps due to zeros of $H(e^{j\omega})$.


GCC-PHAT (Generalized cross correlation) always peak at delay=0 on real audio signal I have studied the GCC-PHAT algorithm to estimate TDOA on audio signals at 2 mics.

Here is my MATLAB implementation:

function cc = freq_xcorr_phat(x,y)
  n = length(x)+length(y)-1;
  X = fft(x,n);
  Y = fft(y,n);
  R = X.*conj(Y);
  R = exp(1i*angle(R));
  cc = ifft(R);
endfunction

I use a SIMULATED stereo audio signal (where channel 1 is an real audio recording and channel 2 is channel 1 with a specific sample delay, 7 samples in this case) to test this function. The result GCC-PHAT cc is plotted below, where we can see that it shows the expected result, a peak at -7:

BUT, when I use a REAL stereo audio signal(channel 2 delay 15 samples) to test this function, the result GCC-PHAT plot goes weird. We can see that it has a peak at -15, but the peak at 0 is more stronger:

THE QUESTION IS:

Why the second plot peakd at 0 and stronger than peak at -15? It doesn't make sense to me.

p.s.

The plot is actually the middle part of fftshift(cc).
My question might be relevant to this question.
The real stereo audio signal is recorded from an embedded system(MCU).

I think that's because of the DC offset of the signal. A time domain or frequency domain plot provided by you could have been helpful here for the confirmation of your doubt. Following are suggestions

1.try to filter the signals before FFT depending upon the signal band you are interested in

2.Subtract the signal by the mean of the entire signal (which is a reasonable trick to find DC offset) or by the DC offset directly if you know that.

The signal of your's must have less swing amplitude in comparison to DC offset making the '0' point correlation of greater amplitude than the actual delay.


Looking for project ideas to match my interests I have a Bachelor of Science in Electrical Engineering, concentrating in VLSI, and I am currently in the process of switching fields to go into Physics. I love physics, and mathematics. I like pretty problems that need subtle insight to solve, which is why I'm moving from engineering to physics. Even though the prettiness of a problem may be subjective, I like the "feeling" I get when I'm doing physics, which was totally absent when I was working in VLSI.

That said, before I make the transition complete, I would like to learn about other areas of electrical engineering that may be closer in spirit to physics and applied mathematics. I have a vague idea that signal processing and imaging are areas that may be quite different from VLSI. Are there any others? Which area should I dive into to get a feel for a different kind of electrical engineering? (I have matlab on my computer and I would like to make use of it).

I especially like to use mathematics on real world problems, and the reason I didn't enjoy VLSI too much is because it is full of boolean algebra, which is somewhat boring and non-visual. It would be helpful if, along with your topic suggestion, you could tell me what areas of mathematics it would involve. Thanks in advance.

PS - I don't have enough points to add "projects" and "careers" tag, so it would be nice if someone could do that for me.

I am a former physician, having done half graduation in mechanical engineering, then a master degree in product design. Nowadays, I work for a biomedical/bioelectric-equipment company, doing mostly software design (algorithms and workflow), using Numpy (matlab like) for prototyping.

THAT fulfills my need (incredibly similar to yours) to solve new problems, work with pretty-picture generation via numbercrunching, programming, and putting the head to perform hard intellectual work very often. I can tell no two months from the three years I am working here have been like one another, and the product portfolio keeps increasing.

In my spare-time, I have interest in geospatial stuff, most probably because I like to ride bikes on relatively remote places with maps and GPS, and the very description of geographical stuff allows for a lot of pretty-picturing. Recently, the avalanche of bike-riding GPS data on sites such as Strava, and the current problems CAUSED by traffic-engineering in big cities makes me think these statistical analysis from people displacement's big data can form a body of knowledge for an "anti-traffic-engineering" oriented for actual people instead of private automobiles.

I can "feel" there's a huge future demand, arriving pretty soon, for software tools that transform datasets into useful information. That can happen with medical imaging, geo-statistical imaging, traffic-engineering micro-modelling, and a lot of fields where data collection, transformation, analysis and interpretation is of huge importance, but still needs further development.

I'm not sure I am actually answering your question, but I feel your question is open to some daydreaming, so I hope this helps, and would be glad to answer any comment you may have!


negative value of signals can we say there is no sound produced when the amplitude is 0?
if yes what can we say when the amplitude is negative?

I believe you are talking about soundwaves only, so I'll approach the answer from there.

The instant values a signal takes don't tell much about the sound it may generate. If a signal is zero at a certain instant $t_0$, it doesn't mean that at that moment there is silence. A single sample can't represent any sound because only oscillations can generate sound. I hope you see why a single sample can't possibly represent  an oscillation of any kind.

Take for example a perfect sinewave. The signal $x(t)=\sin(t)$ is $0$ for infinite values of $t$, but the sound you hear is nonetheless uniform and constant. Also, it takes both positive and negative values, but the note you hear is always the same.

If a signal is $0$ for a given period of time though, then for that interval there will be silence indeed. In fact, if a signal is constant for a given period of time (not necessarily $0$), you will hear no sound at all (but you would be wasting power nonsensically).

In a nutshell, the instant values of a signal represent no sound whatsoever. It's the change over time of these values that can represent and generate sounds.


Is it possible to identify the "blank data" in a JPEG without knowing the image history? I am pretty much a begginer in matters of Image signal and to be honest, not that much versed into mathematics. I hope my question will not sound too naive. (sorry in advance)
(Also if it has been answered already, please redirect me to the correct thread)
Well my question is relating to JPEG quality in general.
rephrased question :
Let's imagine a JPEG with a quality under 100%. Let's say 70%.
If I re-encode this image with a quality of 100%, the weight increase a lot (it more than doubles up on my test with ImageMagick).
Obviously some data has been created, but the image can't be of better quality than the original 70% JPEG.
In the process there is what can be called "blank data" added: data that is not improving image.
How can we know if a given JPEG has an amount of such blank data ?

The short answer is that the original encoding of a picture is lost when you re-encode it.  So a tool that takes a 70% encoded picture and re-encodes it to 100% actually does the following:

Decode from 70% to an internal, raw, uncompressed format (similar to .bmp -- basically, it'll be encoded in RGB, probably 8 bits per pixel for each color).
Encode to 100%

Note two things here: that the encoding assumes nothing about the original image, and that changing the encoding of an image (even from "perfect" RGB to "100% quality" JPEG) is going to introduce scrud -- numerical errors that are for all practical purposes impossible to predict and degrade the image.
So your 70% image isn't perfect -- compared to the original it has numerical scrud.  That -- plus a bit more because the conversion isn't perfect -- will be carried forward into the intermediate uncompressed image.
Then, during "compression" to 100% quality any numerical scrud that's in the intermediate uncompressed image will be treated as real image, not as numerical scrud.  So when the intermediate raw image -- complete with numerical scrud -- is converted to 100% image quality jpeg, that numerical scrud will be replicated as faithfully as can be managed.  Note that, because of the number of steps that the JPEG process goes through, even the best possible job you can do, that results in no compression at all, will still introduce its own scrud.
It actually gets worse than that, but not the way that your question addresses.  If you take an image and compress it, then uncompress it, then compress it again, the second compressed image will be lower quality than the first.  I honestly don't know if for certain levels of compression you'd reach some minimum of image quality at which the thing would no longer degrade, but I'm sure that if you were to compress to some slightly different image quality each time (i.e., choose a random number between 70% and 80% each time), that the image would just continue to degrade with compression -- even as the file sizes pretty much stayed the same for each compression level.


Terminology for instantaneous phase of magnitude component of complex signal In my field, we deal with data that are originally complex-valued. Typically, researchers convert their data from real + imaginary to magnitude + phase, and then discard the phase data (i.e., we generally only deal with the magnitude portion of a complex signal).
There are a number of papers on phase synchrony, in which the Hilbert transform is applied to this magnitude-only data to estimate the analytic signal, from which instantaneous phase is extracted. I am very much a novice in signal processing, so I found this a little confusing. The papers refer to the instantaneous phase signal estimated from magnitude-only data just as "phase signal". These folks might not be aware of the complex-valued nature of the raw signal, which is reasonable given how these data are generally handled, but I think there should be a more appropriate term for "phase estimated from magnitude signal using the Hilbert transform", given that there's very little similarity between the real phase signal and this estimated phase signal.
Is this an established concept and, if so, is there a term that would differentiate the two kinds of phase signal?

and then discard the phase data

What's gone is gone. You can't reconstruct the phase of the original data unless you have some additional information.

in which the Hilbert transform is applied to this magnitude-only data to estimate the analytic signal, from which instantaneous phase is extracted.

That's a just mathematical process. You shove data in, and data comes out. Whether the result is useful or not depends on what you want to do with it and why you chose this mathematical process in the first place. If you want to associate any physical meaning with the result you need to justify this with some physical property or law.

To clarify, the Hilbert transform applied to the real component of a complex signal is not a good estimate of the original complex signal?

Correct. In general most signals are NOT analytic unless there is a strong reason why they should. If there is a reason to interpret the results that way, it should be stated in the paper. In order to "estimate" a discarded phase, you need have some additional information or (justified) assumptions about your original signal otherwise you might as well just toss the dice. If there is a additional information, it should be stated AND the phase estimation method should be properly derived from it.


Nonlinear time-invariant frequency doubler Does there exist a causal nonlinear, time-invariant system that maps each input function $k \cdot\cos(fx)$ to $k \cdot\cos(2fx)$ for all choices of $k, f \in \Bbb R$?
If so, can one represent this system via some sort of differential equation, or in discrete terms, as a nonlinear difference equation?

This would be similar to the second Chebyshev polynomial, but that only works for $k=1$.

As already suggested by Robert and Olli, a system that maps $x(t)=k\cos(2\pi f_0t)$ to $y(t)=k\cos(4\pi f_0 t)$ can be formalized as

$$y(t)=|x(t)|_{max}\left(2\left(\frac{x(t)}{|x(t)|_{max}}\right)^2-1\right)\tag{1}$$

which is a time-invariant non-linear system.

However, I doubt that this system works well (i.e., sounds good as a distortion effect) when applied to real-world (non-sinusoidal) input signals.

To get some ideas of digital implementations of distortion effects, take a look at this thesis and the references therein.


Channel weighting of Mono channel in EBU R128 I am developing a loudness meter conforming to the EBU R128 recommendation.
Currently, it works with stereo files but not with mono files.

In the ITU Recommendation, where the metering algorithms are defined, weightings for the individual audio channels are given (table 3, page 6). From that, I assumed that the weighting for a mono channel is 1.0 (0 dB).

As you can see in figure 1 on page 3, after applying k-filter, mean square and channel weighting to each individual channel, the resulting values of each individual channel are summed to one value. Logic dictates that when having two channels of the same audio signal, the value will be higher than when just having one.

As a consequence, a mono file will result in short-term loudness values, that are 3 LU lower in comparison to a stereo file with same content (on both L and R channels). This is clearly wrong, because established meters like "Klangfreund LUFS Meter" and "MeldaProductions MLoudness Analyzer" will output the same short-term loudness with both files.

How am I interpreting the algorithm definition wrong? How to deal with mono signals?

Side note: It seems that FFmpeg once had a very similar issue, but I could not find the solution to it.

Use a channel weight of 1.0 for your mono signal.

The weighting of 1.5 presented in ITU-R BS.1770 addresses the fact that a signal coming from behind is perceived louder by a human than if it is coming from the front (Reason: To better detect an enemy sneaking at you from behind). With a single mono signal this fact doesn't need to be addressed.

The 3 LU difference you mentioned also makes perfect sense from a physical point of view. One speaker in front of you compared to two (equal) speakers only provides half the energy and therefore the SPL will be lower by 3dB.

I'm honored that you use my LUFS Meter to compare to, but I'd recommend to also use the test files from ITU and EBU. The LUFS Meter manual contains the links to these test files.

Are you aware that the sources of the measurement algorithm of the LUFS Meter as well as some development notes are available? (Sadly I'm not allowed to post a third link... They are at docs/developmentNotes/)

Disclaimer: I'm the developer at Klangfreund.


Get histogram of edge lengths/orientations in an image using OpenCV In the process of trying to estimate whether a ROI contains [Japanese] text, I need to get a rough count of edges, hopefully binned by length and/or orientation... I think this would give me strong-enough features to judge if the ROI is likely to contain text or not.

However, I can't seem to figure out an efficient (this eventually has to work on a video stream) way to get such a histogram (or even just a count of edges for certain length/orientation threshold): Hough transforms will tend to focus on certain shapes (lines etc.) while Canny and Sobel will detect any type, but their output is not particularly useful for counting...

I would appreciate any suggestion or hint at some algorithm I can use, to quickly extract a list of edges from an image...

I don't know if OpenCV has some implementation - but i can provide you with a pointer of algorithm for this.

MPEG 7 provides a fairly well defined (and well thought out) edge histogram descriptor [EHD] which essentially serves your purpose.

You can refer to this: Efficient Use of MPEG-7 Edge Histogram Descriptor which i think is an excellent paper on this.

If you want to see some implementation of MPEG-7 you can find it on JOANNEUM RESEARCH MPEG 7 page.

For more information on MPEG-7: List of possible image features for content based image retrieval


What is Multi-band response? What is meaning of multi band response?
I have already tried reading filters texts but no help.
How it is different from simple band pass or band stop

A multi-band filter is, as the name implies, a filter with multiple pass bands and stop bands. Here is an example of a multi-band filter design.


Convolution effects width of the signal? Let's say there are two signal with different frequency:
\begin{align}
X_1(\omega) &= 0\quad\text{for}\quad \lvert \omega\rvert > 1000\pi\\
X_2(\omega) &= 0\quad\text{for}\quad\lvert \omega\rvert > 2000 \pi\\
\text{And}\quad y(t) &= x_1(t) \star x_2(t)
\end{align}

I understand that convolution in time domain means multiplication in frequency domain, but it does not affect the length of the signal. But in this case, what is the maximum frequency of $Y(\omega)$?
If I were to find maximum sampling interval $T_s$ then which maximum frequency do I have to use?

If $y(t)$ is the signal resulting from the convolution of $x_1(t)$ with $x_2(t)$ then it will have the same bandwidth as $x_1(t)$. $x_1(t)$ has (presumably one-sided) frequency support in $[0, 500]$ Hz and $x_2(t)$ has frequency support in $[0, 1000]$ Hz. So, when performing pointwise multiplication in the frequency domain, the higher frequency components in $x_2(t)$ between $500$ Hz and $1000$ Hz will be multiplied by zero (the value of $X_1(f)$ at those frequencies). The Nyquist frequency required to capture $y(t)$ without distortion should be any frequency greater than $1000$ Hz (twice the highest frequency in the signal).


Is there a better Fourier Transform-based algorithm to use in Image recognition I have been trying to use the Fourier transform in recognizing images of the same size 200x400.

I have tried many different ways to do that such as:

Performing the Fourier on the full image and then try to match it with the ones in the database. the accuracy rate was very bad.
Performing the Fourier on the full image and then take only the
sub-matrix contained in centre of the image and then match it with
the ones in the database. the accuracy rate was also very bad.
Performing the Fourier on the full image and then take only the
horizontal and vertical vector that intersect in the centre of the
image, and then match it with the ones in the database. also, the
accuracy rate was very bad.

Is there any better algorithm that might help me solve this problem?

Thank you very much

The Fourier Transform is good for picking up some similarities between images:

Shifts in $x$ and $y$ are just phase changes in the Fourier domain.
Rotations in the $x-y$ plane cause rotations in the Fourier domain.

So, if this is what you are trying to achieve, then it makes sense.

See this paper for a performance analysis of doing this using the following algorithm


Fourier transform of ${\Pi}_{a}(t)\cos(2{\pi}f_{0}t)$ I want to find the following: $\mathfrak{F}[{\Pi}_{a}(t)\cos(2{\pi}f_{0}t)]$

I first did: $\mathfrak{F}[{\Pi}_{a}(t)] = 2a\text{ sinc}(2af)$
then: $\mathfrak{F}[\cos(2{\pi}f_{0}t)] = \frac{1}{2}[\delta(f-f_{0})+\delta(f+f_{0})]$

So: $$\begin{align}\mathfrak{F}[{\Pi}_{a}(t)\cos(2{\pi}f_{0}t)] &= \mathfrak{F}[{\Pi}_{a}(t)]\ast \mathfrak{F}[\cos(2{\pi}f_{0}t)]\\&=2a\text{ sinc}(2af)\ast\frac{1}{2}[\delta(f-f_{0})+\delta(f+f_{0})]\\&=a\text{ sinc}(2af)\ast\delta(f-f_{0})+a\text{ sinc}(2af)\ast\delta(f+f_{0})\\&=a\text{ sinc}(2af)+a\text{ sinc}(2af)=2a\text{ sinc}(2af)\end{align}$$

I would like to know if it is correct and if there is an alternative way to find the result easier.

Your last line is wrong. Note that you have

$$H(f)\star\delta(f-f_0)=H(f-f_0)\tag{1}$$

for any $H(f)$.

So your last line should be

$$\begin{align*}
&a\,\text{sinc}(2af)\star\delta(f-f_{0})+a\text{ sinc}(2af)\star\delta(f+f_{0})=\\&a\,\text{sinc}(2a(f-f_0))+a\,\text{sinc}(2a(f+f_0))\tag{2}
\end{align*}$$

You basically get copies of your original spectrum centered at $f_0$ and $-f_0$.


How to crop a region with smooth edge and without holes? I'm trying to extract a single region from the image in matlab. The light region looks well connected and distinct from surrounding.

However, I applied a hard threshold B = A(A>1e-3) and it produced

The edge of the extracted shape is not smoothing at all and there are lots of isolated pixels.

For an image as simple as this, it might work to just perform opening or closing on the image to smooth out the edges before thresholding. Opening will remove any small, disconnected parts of the object, leading to an undersegmented image (smaller than the actual size). Closing will fill in any small holes in the object, leading to an oversegmented image (larger than the actual object size).

However, I think your image here is clean enough that thresholding should still work. Try a smaller value for the thresold, or use a MATLAB function like imbinarize to threshold it for you.

ex. thresholded_image = imbinarize(image, 'global');


Bandpass Filter window length: filter "in phase" with sinusoidal input Say you have a bandpass filter and you place it on 40 discretely sampled points from a sinusoid with a period of 40.

Why is the bandpass filter (green curve on the image) only "in phase" with the sinusoid if you use a window length (impulse response length) for the filter of 1/4th the period of the sinusoid (here 40/4=10)? So the green filter has the same turning and bending points only if you use the 10 past observations, multiply with the impulse response (of length 10) and add that up..

Am I correct that this is related to how the length of the impulse response affects the phase delay?

I find what you are saying really hard to understand. Let me try to decode it. Let me know if my decoding is correct.

  Say you have a bandpass filter and you place it on 40 discretely sampled points from a sinusoid with a period of 40.

This makes little sense.  What I think you mean is: I have a filter with impulse response $h(n)$ that I am going to apply to a sinusoidal signal of length 40 samples.

  Why is the bandpass filter (green curve on the image) only "in phase" with the sinusoid if you use a window length (impulse response length) for the filter of 1/4th the period of the sinusoid (here 40/4=10)?

Again, I'm finding it hard to fathom what you mean, but I interpret this as: If my impulse response $h(t)$ is of length 10 samples, then the output of my filter is a sinusoid that is in phase with the input.  When I use an impulse response of a different length (9 or 11) then the phase is incorrect.

For a start, your filter output will not generally have precisely the same phase as the filter input.  If your filter is a linear phase FIR filter, then as @MattL says, the filter will hvae a group delay of $(N-1)/2$ samples --- 4.5 samples in the case of a  length 10 filter.

Trying to reproduce your diagram (scilab code below) shows this effect: the green plot is the original, the red plot is the output of a length 9 filter (which has a delay of 4 samples), the black plot is the output of a length 10 filter (delay 4.5), and the blue is the output for a length 11 filter (delay of 5).

CODE BELOW

P = 40;
x = sin(2*%pi*[0:79]/P)

h10 = ones(1,10)/10;
h09 = ones(1,9)/9;
h11 = ones(1,11)/11;

y09 = filter(h09,1,x);
y10 = filter(h10,1,x);
y11 = filter(h11,1,x);

clf;
plot(x,'g');
plot(y10,'k');
plot(y09,'r');
plot(y11,'b');


Hough lines and Convex Hull methods give jagged lines I need to do a perspective transformation on an image section from newspaper, here the girl's image. The section's coordinates are given as rough estimates. I need to find the boundary of the section and do a perspective transformation. I had been referring the following link.

  http://opencv-code.com/tutorials/automatic-perspective-correction-for-quadrilateral-objects/

However, Hough Transform and Convex Hull method on contour finding give jagged lines. This is how the edge detection on the section looks like.

edges = cv2.Canny(gray, 50, 150)

How can I improve the Canny Edge lines to make them straighter so jagged lines are removed? The Convex Hull finds the boundary but not with a number of smaller lines and the boundary is not complete. Here is how the Convex Hull method's boundary looks like.

contours, hierarchy = cv2.findContours(edges, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
    for i, cnt in enumerate(contours):
        if hierarchy[0,i,3] == -1 and cv2.contourArea(cnt)>5000:
            hull = cv2.convexHull(cnt, returnPoints = True)
            print "hull", hull
            for tempodd, itemodd in enumerate(hull):
                if tempodd % 2 == 0:
                    tempeven = tempodd
                    itemeven = itemodd
                    continue
                cv2.line(img, (itemeven[0][0], itemeven[0][1]), (itemodd[0][0], itemodd[0][1]), (255,0,0), 4)
                print "itemodd", itemodd, "itemeven", itemeven

Here since the Canny is getting jagged lines, the convex hull method gives multiple lines too. Point me where I am going wrong. Thank you.

Let's go step by step.

First of all you can remove perspective distortion without camera calibration. e.g.

  Robust Radial Distortion from a Single Image, Faisal Bukhari and Matthew N. Dailey

  Robust Line Based Calibration of Radial Distortion from a Single View, Thorsten Thormahlen, Hellward Broszio, Ingolf Wassermann

For any given image, Hough transform might find multiple lines, many of those would be an indicator of the same edge. But luckily you could easily merge the lines which overlap. To detect the overlap you could well use the angle between the lines and the distance of lines. This way, you obtain a clustering of the lines. In any case, finding a short line on each edge would be sufficient for you to recover the projective quad.
If you are always sure that the orientation of this quad is more or less fronto parallel, then you could as well use directed gradients (X or Y) to emphasize the vertical or horizontal lines, during the edge detection step.

I agree with the other comments on Gaussian smoothing and anti aliasing. I would not recommend the use of convex hulls though, as they are very noise sensitive.


FFT based symbol synchronization for digital demodulation Where can I read more about "slow search" methods for symbol synchronization based on FFT for extracting a clock signal from the modulated signal?

I have read here (8.7 Symbol synchronization)

  One method of synchronization is to extract a harmonic of the symbol
  frequency from the received signal. Then a local symbol clock can be
  synchronized by methods that are very similar to the phase-locked
  loops used to recover the carrier phase. If necessary, a start-up
  procedure, such as one that uses a slow search, can be used for
  initialization. Synchronization is then maintained by locking a
  feedback loop to a clock signal that is extracted from the modulated
  waveform.

Then here I found a diagram which seems to be the kind of algorithm I am looking for:

  The timing tone can be extracted by ... computing DFT at the symbol
  frequency (i.e., a single point of the DFT output is needed for each
  data block)

There must be more literature about this technique? I am not interested in high-performance methods for synchronization, but rather recovering the clock signal by any means without conserving memory or processing time.

I figure the basic technique is to

Pass the signal through a matched filter
"Condition" the signal (and make it purely real) by computing the magnitude (as in the diagram)
Passing the signal through a DFT computation
Searching the DFT bins for the highest power frequency
Use the discovered frequency of the DFT bin to generate a local timing signal
Use the discovered phase of the DFT bin (as in the diagram) to adjust the local timing signal offset

Symbol timing synchronization seems to be a complex topic although once you get some basic principles right, it all makes simple sense. The method you have referred to is known as Digital Filter and Square Timing Recovery$\ ^{[1]}$, also referred to as Oerder and Meyr algorithm.
EDIT:
And later steps in your summary are not correct. There is no search, timing phase is extracted from the first harmonic. Everything else except the DC term is zero due to the signal being bandlimited.

[1]: M. Oerder and H. Meyr, "Digital filter and square timing recovery," in IEEE Transactions on Communications, vol. 36, no. 5, pp. 605-612, May 1988.


Is it numerically more stable to implement filtering as multiplication or convolution? I'm writing a program to filter a 20,000-sample signal with a fifth-order Butterworth filter offline. I've got access to an FFT implementation. There seems to be two alternatives for implementing the filtering:

convolving the signal with the impulse response in the time domain, or
multiplying the signal with the impulse response in the frequency domain and inverse-transforming the result

These methods would be identical in the theoretical FT case. Doing it in real life with the DFT, though, I suppose things are different. Is one of the methods numerically more stable? Are there any other problems that I should be aware of? The number of calculations is not important.

With convolution, you are not going to run into any stability problems, because there is no recursive filtering, so you are not going to accumulate any errors. In other words, the system is all zeros, no poles. I have heard anecdotally, but not verified for myself, that FFT-based convolution has lower error than time-domain convolution, simply because it has O(n log n) arithmetic operations rather than O(n^2).

Typically, as far as I am aware, Butterworth filters are implemented as recursive (IIR) filters, so that is a different topic. IIR filters have poles as well as zeros, so there can be stability issues in practice. Also, for IIR filters, FFT-based methods are not an option, but on the upside, IIR filters tend to be very low order.

As far as stability issues with IIR filters, they tend to have problems at higher orders -- I'll just throw out a number and say roughly 6th order is pushing it. Instead, they are typically implemented as cascaded biquads (2nd order filter sections). For your 5th order filter, write it as a z-domain transfer function (it'll be a 5th degree rational function), and then factor it into its 5 poles and 5 zeros. Collect the complex conjugates, and you'll have two biquads and one first-order filter. In general, stability problems tend to crop up as the poles get closer to the unit circle.

There can also be issues with noise and limit cycles in IIR filters, so there are different filter topologies (i.e. direct form I, direct form II) that have different numeric properties, but I wouldn't overthink this point -- just use double-precision and it'll almost certainly be good enough.


Lowpass : What is wrong with average filtering? The question is : What is wrong with averaging as low pass filter ?

The details :
I want to lowpass filter a signal to downsample it.
The constraints are : I have no RAM available and I work in streaming, therefore I cannot use a weighting window to filter.
Consequently, I have to do either an average on the last X samples, or a classic filter like :

signal_t+1 = signal_t * (1-factor) + new_measurement * factor

with factor a value adapted to the sampling.
What is more, I know there are signals around a frequency just above the downsampling frequency.

Is it better to use the average or to use the 1-order filter described above ?

More details and tests :
I work with real samples, not complex. Sampling frequency = 200 Hz, frequency after downsampling = 4Hz.

Let's compare a 1-order filter, with a cutoff frequency of 1Hz to limitate aliasing, with two averaging filters. The first averaging filter is the average of 50 samples to downsample from 200Hz to 4Hz. The second averaging filter is an average on 66 samples to get as much rejection as with the 1-order filter.

I think that the filter "average on 66 samples" is the best :
The 50-samples average filter has a rejection higher than the 1-order filter, and I really need to limitate aliasing. But the 66-samples filter has a cutoff frequency larger than the 1-order filter and a rejection at least equal after 2Hz. What is more, the 66-samples filter cuts the frequencies around 3Hz and I know there will be signals on 3Hz that will make my signal noisy after the 4Hz downsampling - I work with real samples, not complex.

Then, I would chose the 66-samples averaging filter.

Is it something wrong in the logic above ?
Would you have suggestions ?

Here is the code :

import numpy as np
import matplotlib.pyplot as plt

nfft = 8192
fsampl = 200
dtsampl = 1 / 200
fcutoff = 1
dtcutoff = 1 / (2 * np.pi * fcutoff)
faverage1 = 4
faverage2 = 3

print("generate signal")
raw_signal = np.zeros(nfft)
raw_signal[int(nfft/2)]=1

print("filter signal")
filtered_signal = 0.*np.copy(raw_signal)
for i in range(1, np.size(raw_signal)):
    filtered_signal[i] = filtered_signal[i-1]+ (raw_signal[i-1]-filtered_signal[i-1]) * dtsampl / dtcutoff

averaged_signal_1 = np.copy( raw_signal )
num_average = fsampl / faverage1
for i in range(1, np.size(raw_signal)):
    averaged_signal_1[i] = np.average( raw_signal[int( max( 0, i - num_average ) ):i] )

averaged_signal_2 = np.copy( raw_signal )
num_average = fsampl / faverage2
for i in range(1, np.size(raw_signal)):
    averaged_signal_2[i] = np.average( raw_signal[int( max( 0, i - num_average ) ):i] )

print("get FFT")
raw_signal_fft = np.fft.rfft(raw_signal)
filtered_signal_fft = np.fft.rfft(filtered_signal)
averaged_signal1_fft = np.fft.rfft( averaged_signal_1)
averaged_signal2_fft = np.fft.rfft( averaged_signal_2)

raw_signal_fftlog = 20*np.log10(np.abs(raw_signal_fft))
filtered_signal_fftlog = 20*np.log10(np.abs(filtered_signal_fft))
averaged_signal1_fftlog = 20*np.log10(np.abs(averaged_signal1_fft))
averaged_signal2_fftlog = 20*np.log10(np.abs(averaged_signal2_fft))

print("plot")
freq_axe= np.linspace( 0, fsampl / 2, int( nfft / 2 ) + 1 )
plt.plot(freq_axe, raw_signal_fftlog, '-')
plt.plot(freq_axe, filtered_signal_fftlog, 'r-')
plt.plot(freq_axe, averaged_signal1_fftlog, '-')
plt.plot(freq_axe, averaged_signal2_fftlog, '-')
plt.plot([0.02, 100], [-3, -3], 'k:')
plt.legend(['raw signal (dirac)', '1-order 1Hz filter', 'average on 50', 'average on 66'],  loc="upper right")
plt.ylabel('amplitude (dB)')
plt.xlabel('frequencies (Hz)')
plt.xlim([0,8])
plt.ylim([-30,2])

plt.show()

The question was : What is best between mooving average and 1st order filter when no RAM is available ?

The solution would be : None of them, except for a specific need : rejection of a particular frequency. Otherwise, it is better to use filter of higher orders than 1st order.

Mooving average :

Drawbacks : delay and poor rejection.
Assets : huge rejection at several particular frequencies

1-order filter

Drawbacks : rejection which is poor whatever frequency
Assets : no delay

Higher order filters. Higher order filters include : cascaded 1st order filters (here under an exemple with 5 cascaded 1st order filters), butterworth filters, etc.

Assets : better rejection than 1st order filter and no need of much RAM. Small delay.

View of the result :

Thank you @user883521 and @Florent Ecochard for your help !

Here is the Pyhton code for a 5th order filter (for other orders, cf links here under. for other filters, see signal.lfilter)

import numpy as np
import matplotlib.pyplot as plt

nfft = 8192
fsampl = 200
dtsampl = 1 / 200

fcutoff_1st_order = 1
dtcutoff_1storder = 1 / (2 * np.pi * fcutoff_1st_order)
fcutoff_5thorder =2.635
faverage2 = 3

print("generate signal")
raw_signal = np.zeros(nfft)
raw_signal[int(nfft/2)]=1

print("filter signal")
filtered_signal = 0.*np.copy(raw_signal)
for i in range(1, np.size(raw_signal)):
    filtered_signal[i] = filtered_signal[i-1]+ (raw_signal[i-1]-filtered_signal[i-1]) * dtsampl / dtcutoff_1storder
filtered_signal5= np.copy(raw_signal)
x = np.exp( -2 * np.pi * fcutoff_5thorder / fsampl )
for i in range(5, np.size(raw_signal)):
    filtered_signal5[i] = np.power( 1 - x, 5 ) * raw_signal[i - 1] \
                          + 5 * np.power( x, 1 ) * filtered_signal5[i - 1] + \
                          -10 * np.power( x, 2 ) * filtered_signal5[i - 2] + \
                          +10 * np.power( x, 3 ) * filtered_signal5[i - 3] + \
                          -5 * np.power( x, 4 ) * filtered_signal5[i - 4] + \
                          +np.power( x, 5 ) * filtered_signal5[i-5]
averaged_signal_66 = np.copy( raw_signal )
num_average = fsampl / faverage2
for i in range(1, np.size(raw_signal)):
    averaged_signal_66[i] = np.average( raw_signal[int( max( 0, i - num_average ) ):i] )

print("get FFT")
raw_signal_fft = np.fft.rfft(raw_signal)
filtered_signal_fft = np.fft.rfft(filtered_signal)
filtered_signal5_fft = np.fft.rfft(filtered_signal5)
averaged_signal66_fft = np.fft.rfft( averaged_signal_66 )
raw_signal_fftlog = 20*np.log10(np.abs(raw_signal_fft))
filtered_signal_fftlog = 20*np.log10(np.abs(filtered_signal_fft))
filtered_signal5_fftlog = 20*np.log10(np.abs(filtered_signal5_fft))
averaged_signal66_fftlog = 20 * np.log10( np.abs( averaged_signal66_fft ) )

print("plot")
freq_axe= np.linspace( 0, fsampl / 2, int( nfft / 2 ) + 1 )
plt.plot(freq_axe, raw_signal_fftlog, '-')
plt.plot(freq_axe, filtered_signal_fftlog, '-')
plt.plot(freq_axe, filtered_signal5_fftlog, '-')
plt.plot( freq_axe, averaged_signal66_fftlog, '-' )
plt.plot([0.02, 100], [-3, -3], 'b--')
plt.plot([2, 2],[-50, 0],  'k:')
plt.plot([3, 3],[-50, 0],  'k--')
plt.legend(['raw signal (dirac)', '1-order 1Hz filter', '5-order 1Hz filter', 'mooving average 66 samples',
            '-3dB', 'aliasing freq', 'unwanted signals frequency'],  loc="upper right")
plt.ylabel('amplitude (dB)')
plt.xlabel('frequency (Hz)')
plt.xlim([0,8])
plt.ylim([-30,2])

plt.show()

useful links :
http://www.analog.com/media/en/technical-documentation/dsp-book/dsp_book_Ch14.pdf
http://www.analog.com/media/en/technical-documentation/dsp-book/dsp_book_Ch15.pdf
http://www.analog.com/media/en/technical-documentation/dsp-book/dsp_book_Ch19.pdf


How does one make a Zero Latency, Linear Phase, or Natural Phase filter? Looking at FabFilter's ProQ3 plug-in, they have an option at the bottom of the window to switch between Zero Latency and Linear Phase modes (and I think a Natural Phase option too).

I'm familiar with BiQuad filters, I'm pretty confident I understand them quite well... but I don't know how FabFilter achieve these different latency and phase modes!

If anyone has any info or references about how one can go about implementing these different types of filters, It'd be much appreciated!

To add to hotpaw's answer, Linear Phase filters can only be implemented with all zero structures (FIR Filters) but they can be approximated with IIR filters such as the Lerner Filter (https://ieeexplore.ieee.org/document/1444798) and other equiripple phase approaches. In order to achieve linear phase, the coefficients of the filter must be either symmetric and assymetric as given in the examples further below.

Zero Latency filters usually refer to Zero-phase digital filtering which is strictly a post processing technique by processing the data through the filter in both forward and reverse directions such that the phase in each direction is cancelled. The filter itself has a transfer function that is the square root of the desired since the signal is filtered twice.  For more information on zero-phase digital filtering, see the filtfilt command available in Matlab/Octave and Python scipy.signal. This is not actually a filter with zero processing delay, but given an input and output data set, the two will be perfectly aligned in time as if the filter had no delay if you later compared the data set once processed (hence post-processing). It will take real time to get the first sample out of the filter that can be placed on a plot with t=0 with the input data. This is very useful for analysis in applications where you would otherwise need to compensate for the filter delay, so worth mentioning here.

With regards to Linear Phase Filters, The following are examples that demonstrate "symmetric" and "asymmetric" taps that is a requirement to implement a linear phase filter.

Linear Phase Type 1 FIR (Symmetric, N odd):

[ 1 2 3 2 1 ]

Linear Phase Type 2 FIR (Symmetric, N even):

[ 1 2 3 3 2 1 ]

Linear Phase Type 3 FIR (Asymmetric, N odd):

[ 1 2 0 -2 -1 ]

Linear Phase Type 4 FIR (asymmetric, N even):

[ 1 2 3 -3 -2 -1 ]

Note importantly how the center bin MUST be zero in order to achieve the required asymmetry for the Type 3. The structures above for each type restrict the zeros of the filter to always be on $z = -1$ for Type 2, and $z = \pm 1$ for Type 3 and $z=0$ for Type 4 restricting which applications can be used for each type.

The following plot further demonstrates this and shows the purpose of having the FIR linear phase "Types":


Identify large solid color areas of image This is a simple image classification problem. I want to categorize jigsaw puzzle images to difficulty levels (e.g, easy, medium, hard).

I want to develop an simple algorithm for this classification, which is based on factors like large solid colored areas, repeating patterns on the image (more details). Such classification can be seen in following website.
example 1

I want to automate this without much human interventions. Is this possible and whats the best way to approach it. (not worried about performance, simple is good, does not need huge accuracy)

EDIT

I was able to generate the magnitude and phase plots of the 2d FFT. However, I am not clear how to come up with a classification by looking at magnitude plot. Any ideas.

Perhaps you could classify the difficulty by assuming that higher spatial frequency implies higher difficulty.  Make sure all the images are the same size, do a 2D FFT, and look at the spatial bandwidth in x and y dimensions, would be my suggestion.


Adaptive Gaussian Filter for Image Denoising I am looking for methods to enhance noisy images, where:

some pixels in the image are very noise,
some other pixels do not contain so much noise.

My first thought is to build an adaptive Gaussian filter.
This means that the Gaussian kernel will depend on the (estimated) noise status of the current pixel (large radius of the Gaussian kernel for noisy pixels and a smaller in the converse case).

Could you help me in calculating the kernel values?

It's not clear in your case that Gaussian kernels will help you.

Since you have a noise shape that is closer to salt-and-pepper noise instead of additive white noise, you may want to try other algorithms (that are somehow linked to an $L_1$ best fit instead of a least squares solution):

you can try median filtering, since it is by design more robust to this noise shape,
you can also try to minimize some TV-$L_1$ functional over your image, see Chambolle-Pock's paper on primal dual minimization for a concrete example,
you also can use a bilateral filter. You will have to play with the exponential decay of the pixel similarity function ($\sigma_d$ in Wikipedia page) to get the desired behaviour.


Computing Hartley transform of zero-padded signal from transform of original signal? If I have the discrete Hartley transform of a signal of length 2n, how can I use it to efficiently obtain the discrete Hartley transform of the same signal padded to length 2n+m for some m > 0?

Since the convolution theorem applies to the Hartley transform as well as the DFT, you may be able to use the same technique as can be used for FFT interpolation or spectrum resampling.  Use the Hartley transform of the implied rectangular window as an interpolation kernel, and interpolate to a larger N. Perhaps window the interpolation kernel for a trade-off between accuracy and computational efficiency.

Depending on the size and efficiency of the interpolation needed to meet your requirements, it may be more efficient to do an inverse and forward transform instead.


the dominant frequency is overwritten in every iteration using matlab code am trying to understand how to find a dominant frequency in a signal
i've a csvfile for an analog data wich detects the breathing rate for the human
i need to find the dominant frequency between 0.1hz and 2hz in this data for every 100 sample , so i wrote a matlab code to do so , after i run the code i find it gives me the same value in every iteration , so i think it must overwrite the value in every index,, here is the code

% clear the current variables in workspace
clear all;
close all;
clc;
filename = 'breathing.txt';
try
    T = csvread(filename);
catch
    warning('no such file');
end

% get the requested signals to process
signalx = T(:,15); % its length is 33167*1

%plot(signalx);

 % ----- apply median filter on data ------

windowSize = 5;
signalx = medfilt1(signalx, windowSize);
plot(signalx);

% -------- process the three signals -------

% initialize variables
samplingRate = 100;
movingwindowSize = 10*samplingRate;
incrementalWindow = 1*samplingRate;
% % % % tolerancefreq = 0.05;
dataLength = size(T,1);
dominantVector = zeros(floor(dataLength/movingwindowSize), 1);
count = 1;
% % % %
for i =1 : incrementalWindow : dataLength-movingwindowSize

    % ------ data segmentation -----

    currentwindow_sigx=signalx(i:i+movingwindowSize-1);

     % ----- apply hamming window -----

     hammingDatax=currentwindow_sigx.*hamming(length(currentwindow_sigx));
     plot(hammingDatax);

     %---- remove dc offset from data ----

     data_x = detrend(hammingDatax);

      % ------- apply fft ------------

     fftData_x = fft(data_x);
     plot(abs(fftData_x));

     % ------- find dominant frequency -----

     len = (length(fftData_x)/2)+1;
     datax = fftData_x(1:len);
     plot(abs(datax));

     fftLength = length(fftData_x);
     freq = (0:(samplingRate/fftLength) :(samplingRate/2));
     plot(freq,abs(datax)), grid('on'), title('dominant freq')

        freq =  freq(find(freq>=0.1 & freq<=2));

    [v,k] = max(datax);

      dominantFrequency = freq(k);
      dominantVector(count) = dominantFrequency;
     count = count +1;

end

so when i run the code it gives me a 322*1 vector with all values == .2000
am sure there is something wrong with the foor loop and the index but i can't find where te problem is.
and here is a photo of the workspace and sample of the overwritten values.

You are searching for the maximum on the entire epoched/windowed data-set:

[v,k] = max(datax);

getting the value and index to this maximum. However you should restrict your search on the frequencies of interest:

valid_idx=(freq>=0.1 & freq<=2);
[v,k] = max(datax(valid_idx));

I would leave the find out as you can perform the search on a logical-index variable and thus you do not need to transform the found-index into the original index.

Additionally, I would first demean your data:

data_x = currentwindow_sigx-mean(currentwindow_sigx);

Which would be equal to detrend(currentwindow_sigx,'constant'), without the parameter constant you are detrending and not demeaning. After this I would apply the hamming window. And then calculate the fft.
You should inspect your data afterwards and check if the data looks "physiological", especially in the frequency band of interest. Here you need some experience to judge this. If the data looks okay the maximum search should give you the correct result.

I assume that the fft will be noise, you could apply some smoothing to it before you search for the dominat frequency to avoid outliers.


can't get the right spectrum in scilab I'm trying to write a simple script that should plot the spectrum in scilab, to test it I use a  sinus function with 440hz so that I get my dirac in this position , my problem is that it doesn't work and I don't understand why?
here's the code :

Fs = 8000;
f = 440;
t= 0:1/Fs:1;
y = sin(2*%pi*f*t);

nf = 1024; // number of point in the DFT
Y = fft(y)
f = Fs/2  * linspace(0,1,nf/2+1);
clf();
plot(f,abs(Y(1:nf/2+1)));

and this is what I get :

any idea why I get this ?

Your fft vector Y has the same length as your input signal y, because you just specified nf without letting fft() know what your desired FFT-length is. This is why your peak does not appear at 440 Hz. Your vector Y does not correspond to the frequencies in f. It's just a matter of correct scaling of the x-axis.

EDIT: I do not know scilab, so I don't know how to pass the desired FFT-length to fft(). If you can't do that, you just need to make your time domain signal have the desired FFT-length.


Capture legitimate huge increase/drop in fuel level from noisy measurements I have a fuel sensor installed in a lorry fuel tank. The sensor measures the current fuel level (in percentage) of the tank. However, due to the fact that the fuel is sloshing in the tank, the fuel level fluctuates wildly — sometimes more than 10%.

I decided to use Kalman filter to predict the actual fuel level, using the following model:

F(k+1) = F(k) - C
C = fuel consumed between k and k+1, which is estimated using
    the distance travelled between k and k+1 (measured using GPS),
    and the average fuel efficiency (fixed) of the lorry (e.g. 35 L/100km).

So far, Kalman filter works perfectly when the lorry is moving and consuming its fuel. However, it fails when 2 events happened that are not captured by the model:

refuel (large increase in fuel level in a short span of time)
fuel theft (large decrease in fuel level in a short span of time)

You can see an example below, on the left side of the graph, where a huge drop in fuel level (most probably due to fuel theft) and a huge increase in fuel level (due to refuel) are smoothed out by Kalman filter, thus the information are lost. On the right side of the graph, where the lorry is moving and consuming fuel, the prediction works well.

So my question is, how do I smooth the fluctuations of the fuel level, and at the same time, capture refuel and fuel theft events? Any ideas? Maybe use something else other than Kalman filter?

Presumably, fuel theft and re-fueling require the vehicle to be stopped and its engine turned off (?). Every time that you start the vehicle up, re-set the system to be tracking the fuel level in the tank at start-up. Comparing the "fuel level at (previous) shut down" to the "fuel level at start up" could tell you what is going on.

Sloshing can persist of course, even during brief stops. Adding a low pass filter / integrator to the sensor would improve instantaneous readings.

It would be interesting to obtain the Fourier Transform of the fuel level sensor signal while the vehicle is moving, provided that it is sampled quickly enough. Depending on your specifications for fuel level accuracy, you might find that a simple low pass filter / integrator might be equally effective in "dampening" the sloshing effect on the signal. This is how common "anti-slosh" modules work anyway.

Hope this helps.


How Matched filter is transformed due to frequency modulation? Intro:

Consider a signal source S which emits a PAM modulated signal, bits encoded as +g(t) and -g(t) pulses.

The matched filter for S is f - f is a function of g.

Then S is fed to an FM modulator and yields a modulated signal, S', let's call the corresponding matched filter of the modulated signal f'.

Question:

If f is known how can f' be computed ?

Normally you wouldn't try to do the match filter directly on the received FM signal.  Typically the receiver path would look something like this-

Received Signal -> Noise/Interference Filter -> FM Demod -> Matched Filter


How to check if Fourier components are in phase of out of phase? I have a time series (of 23 timestamps) of which I take the Fourier transform. Now the fourier transform has 23 imaginary values and each has an amplitude and a phase. When I get the phase angle, it given me arctan value of the angle. Now how can I determine if two fourier coefficients are in phase or out of phase.

Depends on the phase reference point in time.

Depending on the reference point, every DFT basis sinusoid can be either in-phase or out-of-phase (by 180 degrees) with every other basis, thus the same with any coefficient weighting of these basis vectors.

If you reference phase to the center of a DFT (using an fftshift), then the angle will tell you whether a particular component has even or odd symmetry relative to the center.


Main differences to take into account between continuous and discrete time signals I started studying DSP and the first things that came out were the differences between continuous and discrete time signals. So, I was wondering if I understood well these concepts before I keep going into more difficult topics.

The first thing I realized is that discrete-time signals replace the continuous parameter $t$ (time) for $n$ (sample) where $n$ is always an integer, leading this to two main differences:

Complex exponential sequences $Ae^{j\omega n}$ with frequencies $\omega + 2\pi k$ where $k$ is an integer are indistinguishable from one another since $e^{j2\pi n} = 1$  always. This is not true for continuous time signals since
 $e^{j2\pi t}$ can possibly be different from $1$.

Concerning periodicity, in the discrete time case, the fundamental period $N$ is always an integer equal to $\frac{k}{F}$ where $k$ is an integer, meaning that $N$ is not necessarily equal to $\frac{1}{F}$.

In the continuous case, this is not true since $T = \frac{1}{F}$.

So, my questions are:

Am I correct about this two differences? If this is not the case, please
explain it to me better.
Are these two the more important to take into account, or am I
missing something (else) big?

In general I'd say that you understand the differences pretty well.

You nailed the first difference- continuous time vs. discrete time (samples).
You were also correct about there being sets of indistinguishable frequencies in the discrete realm.  A corollary of this is that all discrete signals are band-limited, whereas continuous signals do not have to be.  This concept is most often referred to as the Nyquist frequency.
I don't entirely agree with your point about periodicity, though I believe I understand what you're getting at.  If by periodicity you mean having the exact same values repeated, then yes, a discrete period would have to be an integer number of samples.  On the other hand, if you take periodic to mean that the "underlying" signal repeats, then it is possible to have the samples not be the exact same if the period is not an integer multiple of samples.  You could still detect the periodicity with an autocorrelation though.  Anyway, if this point is confusing, don't worry about it- it's not terribly important.

It is, by the way, possible to view discrete signals as a special case of continuous signals.  You can model discrete signals as continuous signals that have been multiplied by an infinite set of evenly spaced delta functions.  If the delta functions are separated by $T_s$ seconds, then in the frequency domain they can be represented by another set of delta functions that are separated by $\frac{1}{T_s}$ Hz.  The multiplication in the time domain is equivalent to convolution in the frequency domain, so copies of the signal are reproduced in the frequency domain at each of the delta functions.

Though the mathematicians in the audience will probably scream at me for saying this, I think that you can likewise look at continuous signals as a special case of discrete signals, where we take the limit of $T_s$ as it approaches zero.


Why does dot product return zero at non-harmonic frequency If I understand correctly, dot product between two sinusoids should return zero if they are orthogonal. Since dot product is at the heart of the DFT, we can only clearly compare the frequencies that are multiples of the frequency from the equation: (sample frequency / number of samples). If our sinusoid is not a multiple of this frequency, we get spectral leakage. This all seems fine, but I can't understand why does the dot product return zeros at some non-harmonic frequencies as seen in the third example. Let's suppose we have the following:
Sampling frequency = 512 Hz
Number of samples = 512
Our fundamental frequency is then: 512 Hz/512 = 1 Hz, so our DFT bins would represent frequencies from 0, 1, 2, ..., N - 1.
Example 1:
sin1 = real sinusoid with frequency 5 Hz
sin2 = real sinusoid with frequency 5 Hz
Dot product with these two real sinusoids returns 256 as expected.
Example 2:
sin1 = real sinusoid with frequency 5 Hz
sin2 = real sinusoid with frequency 6 Hz
Dot product returns zero as expected.
Example 3:
sin1 = real sinusoid with frequency 5 Hz
sin2 = real sinusoid with frequency 5.5 Hz
Dot returns zero, but I don't understand why. Why do we get zero here instead of a non-zero number that would represent spectral leakage (as the second sinusoid is not a multiple of a fundamental frequency)?
Example 4:
sin1 = real sinusoid with frequency 5 Hz
sin2 = real sinusoid with frequency 5.75 Hz
We get non-zero value, so this works as expected as we get spectral leakage.
Is my understanding wrong?

Note that when you compute the DFT, you compute the dot product with a sine term and with a cosine term. In your third example, even though the dot product with the sine term is zero, the dot product with the cosine term isn't, so you do get spectral leakage, as you expected. The dot product with the sine term is zero because you have a whole number of periods plus one half period inside the DFT frame, so the periodized version of the signal is even and, consequently, the imaginary part of its DFT (the dot product with the sine term) vanishes.


How to determine the window size and weights in Weighted Moving Average (WMA), given desired cut-off frequency? I am trying to smooth my discrete-time data points using the method of WMA.

Currently, I am using n as the window size and the weight array, {n/(n(n+1)/2), (n-1)/(n(n+1)/2), ... , 1/(n(n+1)/2)}.

If the y-value of each point is irrelevant, I can just simply randomly choose my size n.

However in my case, I hope to reserve the original values of the data points to the best extent. Thus, I cannot choose a big window that averages everything to flat.

My cut-off freq. is 3Hz and the sampling rate is 50Hz.

How may I choose the size of the window n?

Thanks in advance!

Your normalized window is given by

$$w(n)=\frac{2}{N+1}\frac{N-n}{N},\quad n=0,1,\ldots,N-1$$

The window satisfies

$$\sum_{n=0}^{N-1}w(n)=1$$

which means that the gain of the corresponding moving average filter is 1 at DC.

For determining the cut-off frequency, we need to compute the frequency response of the window:

$$W(e^{j\theta})=\sum_{n=0}^{N-1}w(n)e^{-jn\theta}$$

After some algebra you get

$$W(e^{j\theta})=\frac{2}{N+1}\frac{1-\frac{N+1}{N}e^{-j\theta}+\frac{1}{N}e^{-j(N+1)\theta}}{(1-e^{-j\theta})^2}\tag{1}$$

Now you need to find the value of $N$ for which the magnitude of (1) at the cut-off frequency $\theta_c=2\pi\frac{3Hz}{50Hz}$ becomes $1/\sqrt{2}$ (-3dB). Since $N$ must be integer you cannot get any desired cut-off frequency, but the given cut-off is approximately achieved by $N=9$, for which $|W(e^{j\theta_c})|=0.698$ (-3.13dB).


What does the 420 means in YUV420? I'm not sure if is the correct location for ask this question. If not please info me the correct StackExchange site.
As this wikipedia page show, there is one YUV format of image as YUV420 and YUV411.

I know what is YUV, but what those 420 and 411 means?

The 3-digit number describes the subsampling of the
chroma (U and V) channels.  A detailed explanation is at

http://en.wikipedia.org/wiki/Chroma_subsampling

In particular, YUV420 means that the U and V channels have half the resolution of the luma (Y) channel, in the horizontal and vertical directions. The sampling method usually present in low to medium quality JPEG images is YCrCb420 (high quality JPEGs are usually YCrCb444 (no subsampling).


The small red peak in the CIE standard observer The human eye has 3 types of color cones: red, green and blue. Their spectral sensitivity peaks in short (S, 420–440 nm), middle (M, 530–540 nm), and long (L, 560–580 nm) wavelengths. See diagram below.

Then we have the CIE standard observer. Wikipedia says it can be thought as the spectral sensitivity curves of three linear light detectors yielding the CIE tristimulus values X, Y and Z.

I don't get why the red line in the CIE standard observer color matching function also peaks at 440 nm. To be more specific: What causes the x-line in a CIE standard observer color matching function to peak at around 440nm?

Can someone explain me that?

Probably you've noticed that primarities are $\mathbf{X}$, $\mathbf{Y}$, $\mathbf{Z}$, not $\mathbf{R}$, $\mathbf{G}$, $\mathbf{B}$ (which are corresponding to the color values $R$,$G$,$B$). This is the aftermath of original work conducted by Wright and Guild yielded the Color Matching functions: $r(\lambda)$, $g(\lambda)$, and $b(\lambda)$ having a negative value of the red primarity around $522 \;\texttt{nm}$ (figure below). Everyone will agree that it is non-intuitive and somewhat confusing.

Why values of $r(\lambda)$ are negative? Well after tristimulus color matching experiments, it turned out that not every color can be created with the preset primarities ($700\; \mathtt{nm}$, $564.1\; \mathtt{nm}$, and $435.8\; \mathtt{nm}$). So for example, if you wanted to create the blue-green color of wavelength $\lambda = 490 \; \mathtt{nm}$ it cannot be done only by adding blue and green primarity alone. You also need to subtract the red one. Mathematically speaking, the color equation:

$$\mathbf{C}=r\mathbf{R} + g\mathbf{G} + b\mathbf{B} $$

Must by modified by adding the outer mixture $-r_{out}\mathbf{R}$ to both side of the above equation:

$$\mathbf{C} + r_{out}\mathbf{R} = g\mathbf{G} + b\mathbf{B} $$

That is why CIE introduced new set of so-called virtual primarities: $\mathbf{X}$, $\mathbf{Y}$, $\mathbf{Z}$. In order to obtain only positive or null values of the Color Matching Functions. That's how $2^{\circ}$ observer Standard Color Matching Functions (also called CIE 1931) were created. By doing so they ensured that (amongst other assumptions):

New SCMF's  $x(\lambda)$, $y(\lambda)$, $z(\lambda)$ have values $\ge 0$ (not like $r(\lambda)$).
Are related to previous CMF's ($x$,$y$,$z$) by a appropriate transform.
$\mathbf{X}$, $\mathbf{Y}$, $\mathbf{Z}$ must be independent (kind of orthogonal).
Only color value $Y$ is proportional the the lightness of a colour.

Newly obtained color quantities $X$,$Y$,$Z$ are called Standard Color Values. The relation between CMF's and SCMF's is given by a following linear transform:
$$
\left( \begin{array}{c}
x(\lambda) \\
y(\lambda) \\
z(\lambda) \\
\end{array} \right)
=
\left( \begin{array}{ccc}
2.7689 & 1.7517 & 1.1302 \\
1 & 4.5907 & 0.0601 \\
0 & 0.0565 & 5.5943
\end{array} \right)
\cdot
\left( \begin{array}{c}
r(\lambda) \\
g(\lambda) \\
b(\lambda) \\
\end{array} \right)
$$

After over a three decades people from CIE had realised that $2^{\circ}$ observation angle might not be enough for adequate color assessment. That's why in 1946 they introduced the $10^{\circ}$ observer, also known as CIE 1964 observer. Latter is recommended for visual angles above $4^{\circ}$, whereas CIE 1931 should be used for any angle below that. Of course there is no sudden change in perception of colour at $4^{\circ}$, this number is just arbitrary. Comparison of these two standards is depicted of a following plot, on which you can see the 'bump' mentioned by you for the $x(\lambda)$ standard color value:


why Gaussian noise is usually used? Gaussian noise with different SNR levels are usually used in research works to simulate a realistic environment. How can researchers guarantee that Gaussian noise can simulate the reality of a System?

As pointed out in other answers, the Central Limit Theorem is one reason why Gaussian noise is so important as a model. One other reason that hasn't been explicitly mentioned and that I would like to point out is the fact that Gaussian noise is completely described by second order statistics, which are relatively easy to measure. E.g., if you have a univariate distribution and if it is Gaussian, then you know everything you can ever know about it by knowing (measuring) its mean (first order) and its variance (second order). There are no higher order statistics, which is good because they are much harder to measure reliably. This fact is of course no justification for using a Gaussian model, but it is a very handy property if we think that a Gaussian model is justified.


How can we know the geometry of a blob detected by OpenCV I am performing blob analysis using connectedComponentsWithStats, which returns elementary statistics about the blobs (AABB and area). It is also possible to compute geometric moments.

Anyway, I need to compute my own features and for this I need the detailed geometry, be it a chain of outline pixels, a run-lengh-coded description or similar.

I know that the contours can be obtained by other means (findContours), but I would like to know if the connected components engine can help me achieve my goal or if it is a black box forever.

Instead of connected components you may also use blob detector mentioned here

It provides some information on the shape and geometry of blob

For geometric moments you will have to compute contours.

I am not sure what do you mean by connected components being a black box.


Is it possible to simplify the convolution integral if the functions are non-zero in disjoint areas? Possible Duplicate:
  When convolving two functions that are constants in a region and 0 everywhere else, where does the integration start?

I have a function $f(x,y)$ and $h(x,y)$.  $f(x,y)$ has a value of $\frac{1}{3}$ when $x$ is between $\frac{1}{3}$ and $\frac{5}{9}$, and $y$ is between $0$ and $1$.  The function has a value of $0$ everywhere else.  Meanwhile, $h(x,y)$ has a value of 1 when both $x$ and $y$ are between $-\frac{1}{19}$ and $\frac{1}{19}$.

These two functions are non-zero in completely disjoint regions.  Is there a way to leverage this property to simplify the convolution integral?

Link to the closed-form solution:

wolframalpha

From left to right: $f$, $h$, and $f*h$:


What is a good book reference on the definition of "subharmonic frequency"? Is there a good book reference on the definition of "subharmonic frequency" as integer fractions of the highest available frequency?

"Highest available frequency" sounds fishy but here's a reference from dynamics:

Cunningham, Walter Jack. Introduction to Nonlinear Analysis. 1st ed. New York, USA: McGraw-Hill, 1958, pp. 191–192.

ISBN 10: 0070149178  ISBN 13: 9780070149175


Prefiltering to Avoid Aliasing In processing analog signals using discrete-time systems, it is generally desirable to
minimize the sampling rate. This is because the amount of arithmetic processing required
to implement the system is proportional to the number of samples to be processed.
If the input is not bandlimited or if the Nyquist frequency of the input is too high,
prefiltering maybe necessary. An example of such a situation occurs in processing speech
signals, where often only the low-frequency band up to about 3 to 4 kHz is required for
intelligibility, even though the speech signal may have significant frequency content in
the 4 kHz to 20 kHz range. Also, even if the signal is naturally bandlimited, wideband
additive noise may fill in the higher frequency range, and as a result of sampling, these
noise components would be aliased into the low-frequency band.
Could someone explain please the point below and why it is a problem?
Ie, why is it a problem if noise components are aliased into the low-frequency band?

"Also, even if the signal is naturally bandlimited, wideband
additive noise may fill in the higher frequency range, and as a result of sampling, these
noise components would be aliased into the low-frequency band."

Let's look at a simple example: you have speech signal that has a sinusoidal noise component at 7 kHz.
If you sample this at 8 kHz without pre-filtering, the noise will alias down to 1kHz and contaminate your speech signal in a very important frequency region.
If you pre-filter at 4 kHz, the noise will be gone.


resonance demodulation for bearing faults Can someone explain me how are the different frequency bands useful for detecting bearing faults using resonance demodulation?

This presentation has a pretty good overview of characterizing bearing faults. In essence, once the measured frequency content does not match what is considered a good response from a bearing, it can be flagged as faulty or nearing fault.