OpenCV Haar Classifier - best type of positive images I'm currently collating images to try and develop a HAAR classifier using the OpenCV traincascade functionality. The objects I want to detect are small and numerous (bees), and I'm curious as to whether it would be better to train the image with lots of small cutouts of the bees as positive images, or to use bigger positive images with the bees present along with a chunk of the background.

You probably want to train on as tightly cropped rectangular region around the bee as possible.

I'm not sure if you already thought of this but you should remember that this classifier can only handle recognizing objects in rectangles whose axes are aligned with image, i.e. no rotated rectangles. This is going to be tricky if the bees can be facing any which way. A human face detector has the advantage that faces are almost always found right side up in the image. If you wanted to recognize upside down faces too you would probably train only on right side up images and then check your input image twice--once with the original input and once with it flipped upside down.

You probably want to do something similar for the bees. All of your training images should be aligned so that the bee has approximately same orientation and scale. Then if you want to detect different orientations you can rotate your input image.


Image Processing - Counting nuclei I am trying to create a program that can count the number of nuclei in such an image:

What I've already done is the following, step by step:

Apply an Alternating Sequential Filter (closing and opening the image with gradually bigger structuring elements)
Apply a distance transform
Apply watershed segmentation using the distance transformed image to detect minima

Which yields the following result (where each color represents a new nucleus counted):

As we can see, there are lots of imperfections, specifically, overcounted nuclei. I'd say that the reason for that problem is the way I imposed minima for the Watershed Transform (using the distance transform), but I really have no other ideas for imposing minima in that case.

As the Distance Transform generates minima based on the roundness of objects, I would like to know a better alternative to round up the nuclei than the Alternating Sequential Filter (looking at the image above, we can infer that most of the "overcounts" came from the less rounded nuclei). I would also like to know better ways to impose minima for the Watershed Transform.

There are numerous articles on how to handle the oversegmentation problem of watershed, but I think you should read Robust Cell Image Segmentation Methods (scientific article from 2004 by Bengtsson et al).

It covers various methods for segmenting cell images and includes real-world examples that show how to handle oversegmentation from watershed on fluorescence microscopy images similar to yours (it also has examples for bright-field images and confocal microscopy images). It uses seeds from the distance transform, similar to your approach, and merges regions with weak borders. The article reads well and the concepts are pretty straight-forward to implement in Matlab.

For an even more current approach, you can read A Decomposition Scheme for 3D Fuzzy Objects Based on Fuzzy Distance Information by Svensson. It uses a similar method as in Bengtsson et al, but works on the fuzzy distance transform which gives a better density representation for the objects used in the article.


blur detection using opencv I'm writing a script to detect blur images using OpenCV by applying Laplacian filter and calculate the std but there is a problem the std for images that contain motion blur is very close to those images which contain any other type of blur. (my propose is to find detect unsuitable images that can't detect the object on it like a dog image that contain motion blur for example)

I used this piece of code cv2.Laplacian(image, cv2.CV_64F, ksize=3 ).std()
 to find the blur ratio and make a threshold for std < 40 is considered a blurry image
 so I want a method to can differentiate between images that contain motion blur images and other kinds of blur

You can use the following variance of Laplacian responses:

cv2.Laplacian(gray_image, cv2.CV_64F).var()

More details at https://www.pyimagesearch.com/2015/09/07/blur-detection-with-opencv/


Do I need equally paced values to do a Fourier Cosine Transform? I am gathering data from a sensor, but the data gathering depends on a code running, so the timing is not precise to make it equally paced.

What I mean is this, the time it takes to do measurement #1 and #2 is not exactly the same between measurement #2 and #3 and so one.

I may have data like this:

value = 20 taken at t=0s
value = 21 taken at t=0.2s
value = 44 taken at t=0.44s
value = 22 taken at t=0.67s
value = 33 taken at t=0.9s
value = 32 taken at t=1.3s

As you can see the time between measurements is not exactly the same and can vary slightly.

The idea is to perform a Fourier Cosine Transform on the discreet values.

First question is:

will this not equally time between the measurements affect the DCT?
what can be done to "normalize" the data to have constant times?

All standard transforms in DSP (including Fourier transforms or Cosine transforms etc.) assume uniformly sampled (time, space or frequency) input data.

When your data is nonuniformly sampled, you cannot directly use the standard algorithms.

You should either convert them into uniform samples and then use the standard algorithms or use specific algorithms that work on nonuniform data directly.

For the conversion, look for non-uniform to uniform interpolators.


Comparing two colored lines I have two colored lines, which can be thought of as functions I1,I2:[0,1]->R. I want to check whether I1 and I2 'topologically' represent the same object. In other words, if I stretch some parts of I1 and squeeze some other parts, I want to get I2.

Obviously they have to have the same colors and 'in the same order'.

Is there a metric that measures this kind of similarity?

If I consider that your colored lines are multivariate signals (in Red, Green, Blue), that might be related through local offsets and scalings (stretch and squeeze), I would suggest the method in Correlation based dynamic time warping of multivariate time series, Expert Systems with Applications, 2012.

It combines Dynamic Time Warping (DTW), a widely used technique for comparison of time series data, and Principal Component Analysis (PCA)  to address the correlation between the variables.

yet, you have to balance, in your metric, the 'shape' part and the 'color' part of the data. The associated quantities do not live in the same space.


How to reduce PAPR with constellations choice? How to tune constellations mapping M ary-psk or x-qam to reduce PAPR ?

Peak to Average Power Ratio (PAPR) is defined as:

$$
\text{PAPR} = \frac{\text{Peak Power}}{\text{Average Power}}
$$

If you're working with a $M$-PSK constellation, then the peak power is $P_s$ and the average power is also $P_s$ assuming each symbol is equally likely. For $M$-PSK, the $\text{PAPR}=1$.

If you working a $M$-QAM constellation, then the peak power is $\text{max}_{i \in \mathbb{C}} |s_i|^2$, where $\mathbb{C}$ is the constellation. Again assuming all the symbols are equally likely, the average power is $\frac{1}{M}\sum_{i \in \mathbb{C}} |s_i|^2$. The specifics should be worked out for a constellation but is hard to write down because in general there are many possible QAM constellations.

One thing that is directly clear is that $\text{PAPR}_{\text{QAM}}>\text{PAPR}_{\text{PSK}}$. In a QAM constellation, you have some symbols with less power, some symbols with more power, and there is an average somewhere in between. The point is that there are some symbols with more power, and this makes the peak power part of PAPR greater than the average power part. This is in contrast to a PSK constellation, where the symbols lie on a circle and no symbol has more or less power than each other.


How to make a 2D gaussian Blur terms? I want to make a 2D gaussian Blur coeffients.
Now I am making 1D gaussian Blur term as following.

double Cmfc_test5Dlg::Guassian(int x, double sigma)
{
        double c = 2.0 * sigma * sigma;
        return exp(-x * x / c) / sqrt(c * 3.141592f);
}
double Cmfc_test5Dlg::GuassianTerms(int kernalSize, double sigma, int *terms)
{
        for (int i = 0; i < kernalSize; ++i) {
                terms[i] = Guassian(i - kernalSize / 2, sigma)*1024;
        }
        return 1;
}

But I want to make 2D instead of 1D. What should I do ?

It's real easy. Multi-dimensional Gaussian functions are separable, that is they can be calculated as the product of the Gaussian function along all axes. For the 2-d case it means:

terms[x, y] = terms[x]*terms[y]

Note that separability also means that you obtain the same result as by using the 2-d filter, by first 1-d filtering along one axis and then 1-d filtering the output along the other axis.


How to determine if an IIR filter is causal or non-causal? It may sound impractical, as only causal system are passively realizable, but what I am facing is a general statement made regarding recursive digital filters: for them to be stable, the poles should  be inside the unit circle. I am not able to see that how without knowing type of system that the  system is causal or non causal one can comment on location of poles for system to  be stable.

I also know that for system to be stable its output should only be dependent on values of input at or before observation time for instant below is causal

But what about below equation

First of all, it's not correct to say "poles should (always) be inside the unit circle for an LTI system to be stable" ; unless it's implied that system is also causal. Otherwise, if the system is noncausal, then its poles should be outside of unit circle for the system being stable.

For IIR systems that are described by LCCDEs causality must be externally imposed on the system (consistent wit the initial conditions), and cannot be derived from the equation alone. For example, the system

$$y[n] + a y[n-1] = b x[n]$$ cannot be concluded to be causal, as the equation can also be written as

$$ y[n] = (-1/a) y[n+1] -(b/a) x[n+1] $$
which is (apparently) indicating a noncausal dependence on the current output on both future input and future output.

Therefore this difference equation can signify both a causal and a noncausal system. Its solution should be derived by assuming a causal system or noncausal system.

This is also understood by the fact that a given transfer function $H(z)$ will have a corresponding LCCDE representation, but will have multiple ROCs. For each ROC the LCCDE should be solved accordingly, yielding a different solution for each assumption. But the LCCDE is the same.

FIR systems do not have poles (except at origin or infinity) and have no issues with stability and they don't posses regions of convergences on the z-plane. Hence their equation will be signifying the causality such as:
$$ y[n] = a x[n+1] + b x[n] + c x[n-1]$$
is a noncausal FIR system and you cannot manipulate the given equation into a causal form.

Or the following system
$$ y[n] = a x[n] + b x[n-1] + c x[n-2]$$
is a causal FIR system and you cannot manipulate the given equation into a noncausal form...

So causality shall be an assumed property for IIR systems.


Remove background noise from recording I have a sample speech recording with background noise. I managed to extract part of the noise signal, and put the noise and the contaminated signal in the same plot.

I have done fft-analysis on both speech and noise, and since their frequency bands with centralized power are close to each other, I cannot simply apply a bandpass filter. Is there another way to remove the noise from the speech signal with knowledge of the noise itself? (No, subtraction does not work)

The signal seems non-stationary, at least less stationary than the noise, and the noise is of relatively lower amplitude. Hence, a class of methods resides in performing one or several non-stationary transformations (time-frequency, time-scale, or  spectrograms and wavelets). They can concentrate useful signals on a little amount of coefficients. Then, thresholding or shrinkage techniques, reducing low-coefficients and preserving higher ones can be interesting.

Free softwares like audacity have build-in functions. For more techniques, you can check for instance:

VOICEBOX: Speech Processing Toolbox for MATLAB, esp. Spectral noise subtraction.

More advanced techniques, termed "source separation", are for instance described by:

FASST: Flexible Audio Source Separation Toolbox


Interlaced to De-interlaced I have an interlaced image and I like to make a deinterlaced image out of it. Currently I just take average of two rows to create the missing row. Is there a better way to do this? What is the state of the art?

I don't have the previous or the next frames, just this one.

If you want to vertically upsample an image, one method possibility better than linear interpolation is to FFT, zero pad, and IFFT; or a similar process using the DCT/IDCT and/or sub-blocks of the image.


Signal recovery from multiple noisy sensors I'm working on detecting signals buried under the sensor noise using a multiple sensor system and using some technique (correlation?) to take what is common to all the sensors and cancel out the random, uncorrelated part of the sensors.
The signal I want to recover is some 30 times smaller than the rms noise. The smaller I could recover the better, though.
can anyone point to any technique for this case?

Your ability to achieve coherent gain is going to depend on a number of factors and your upper bound on gain is 10log10(N)dB, where N is the number of sensors.

The essential idea is that you need to align your signals in time and add them resulting in a superposition.  The signals need to be nearly identical at each sensor output.

In beamforming, there is a notion of signal source located at some distant point and plane waves are incident on the sensors.  One can infer the anticipated time delays from geometry.  One can also assume other wavefront geometries such as spherical wave incidence.

Effects such as polarization and sensor orientation complicate the similarity of the received signal at each sensor.

Some signals have multiple propagation paths and these paths aren’t all going to be free space paths. Seismic signals have pressure and sheer wave propagation.

The assumption that noise is independent at each sensor isn’t a given for every application.

These are all general considerations and you need to offer more detail on your problem if you want more help.  The 10log10(N) upper bound should apply in any case.


Inverse $z$-transform of a transfer function in MATLAB I have designed a Butterworth highpass filter (HPF) of 4th order with cutoff frequency high enough to give a gain of $3$ at high frequencies. I want to find the inverse $z$-transform using MATLAB.

L=4;
OSR = 16;
[b,a] = butter(L,0.259,'high');
b = b/b(1);
fvtool(b,a);
H = tf(b,a,[],'variable','z^-1');

iztrans doesn't seem to support tf or zpk.

The function iztrans only supports symbolic expressions. This should do the trick:

[num,den] = tfdata(H);
syms z
H_sym = poly2sym(cell2mat(num),z)/poly2sym(cell2mat(den),z);
h_inv = iztrans(H_sym);


how can I generate SIMO underwater acoustic channel I'm working in underwater acoustic channel, then I need to generate SIMO channel to use it for simulation in MATLAB,  how can I generate SIMO channel?

thank you

While the literature is small, the Oceans are large.  There is a large variation in physical detail.  Unlike RF, there are no ITU standard models. A serious effort really requires real data, You can spend 20 years refining a model.

You could start by looking at,

  Stojanovic, Milica, and James Preisig. "Underwater acoustic
  communication channels: Propagation models and statistical
  characterization." IEEE Communications Magazine 47.1 (2009): 84-89.

At 100Khz, a ray tracing model should be adequate and you can get ray models at Porter's Ocean Acoustic Library at:

http://oalib.hlsresearch.com/

Ray shooting can require a lot of trial and error to get eigenrays. At short ranges, an iso-velocity range independent model might be adequate.  Most models are Fortran, not Matlab, but you can call them from Matlab. If you look around the web, particularly those of Australian establishments , there is code with tighter MATALB integration.  Most prop code is deterministic.  Ira Dyer is a pioneer is random propagation, so papers that cite him, might be useful, and you can do that using Google Scholar. If you can read Russian, a lot of their literature is very good.

Noise models are local but at 100KHz, noise sources are going to be near.

Projector and forward volume scattering effects are hard to model. Without data, you need to make some guesses.


Use of Long Training Field (LTF) of 802.11 standard In 802.11a/n the preamble part has L-LTF(Legacy-Long Training Field) and HT-LTF(High Throughput-Long Training Field). I want to clarify my understanding of these two fields:

Are the only difference between these 2 fields are the number of subcarriers used and the symbols used?
This field is mainly used for Channel estimation. What method is used to do this estimation (say on our home wifi router), is it LS/MMSE/MLS or any other in 802.11a and n? If possible can you explain how is it determined?

I'm not to familiar with the HT-LTF, but OFDM channel est. is done as follows:

rx'd long training symbol ($LTF_{rx}$) is the tx'd LTF ($LTF_{tx}$) convolved with the channel:

$$LTF_{rx}(n) = \sum_{m=-\infty}^{\infty}h(m)LTF_{tx}(n-m)$$

translate this to frequency domain:

$$LTF_{rx}(w) = H(w) LTF_{tx}(w)$$

and the channel inverse is then:

$$H(w)^{-1} = \frac{LTF_{tx}(w)}{LTF_{rx}(w)}$$


Kalman Filter - Optimal Way to Handle "Derived" Measurements? Ie, if you have as state variables position (p) and velocity (v), and I make low-frequency measurements of p, this also indirectly gives me information about v (since it's the derivative of p). What is the best way to to handle such a relationship?

A) At the update step, should I only say I've measured p, and rely on the filtering process, and my accumulated state-covariance matrix (P), to correct v?

B) Should I create an "extra" prediction step, either after or before my update step for the measurement of p, that uses my measured p and (relatively large) delta-time to make a high-variance prediction of v?

C) In my update/measurement step, should I say I've made a measurement of both p and
v, and then somehow encode information about their interdependence into the measurement co-variance matrix (R)?

For a little more background, here's the specific situation in which I've run into the problem:

I'm working with a system where I want to estimate the position (p) of an object, and I make frequent measurements of acceleration (a) and infrequent, high-noise measurements of p.

I'm currently working with a codebase that does this with an Extended Kalman Filter, where it keeps as state variables p and v. It runs a "prediction" step after every acceleration measurement, in which it uses the measured a and delta-time to integrate and predict new p and v.  It then runs an "update"/"measurement" step for every (infrequent) p measurement.

The problem is this - I get occasional high-error measurements of a, which result in highly-erroneous v.  Obviously, further measurements of a will never correct this, but measurements of p should get rid of this.  And, in fact, this does seem to happen... but VERY slowly.

I was thinking that this may be partially because the only way p affects v in this system is through the covariance matrix P - ie, method A) from above - which seems fairly indirect.  I was wondering if there would be a better way to incorporate our knowledge of this relationship between p and v into the model, so that measurements of p would correct v faster.

Thanks!

In ideal world you'd have the correct model and use it.
In your case, the model isn't perfect.
Yet the steps you're suggesting are based on a knowledge you have about the process - which you should incorporate into your process equation using your dynamic model matrix:

The classic and correct way given F matrix is built correctly according to your knowledge.

"Extra" prediction step will yield nothing, since $ {F}_{ik} = {F}_{ij} {F}_{jk} $ and if you reduce the time frame you should alter $ Q $ and $ R $ accordingly which should get you at the end of the chain of small steps the same $P_{k \mid k - 1}$.

If you don't measure V you'll must "Estimate" it somehow. Yet by definition, if your case falls under Kalman's assumptions using Kalman's filter would yield best results.

All in all, stick with the "Classic".


Digital Camera Spectroscopy? Would it be possible to use a conventional Canon DSLR camera to do some kind of vaguely meaningful spectrum analysis with GNU Octave?

My niece needs to do a biology science project for school so I was thinking we could use her digital camera to periodically take images of a sample like a test tube of broth+yeast or proteins+enyzme or similar and then do some kind of image processing using Octave to do some kind of crude spectroscopy.

We have the computer doing tethered image capture no problem and I have octave processing the images and doing basic fft type stuff. But I'm not terribly familiar with signal processing and definitely not image processing so what would you recommend given the problem?

I am thinking the image should be processing in a homogeneous way. Meaning no one area of the image is more important than another. So some kind of averaging might be done first? Then an FFT?

How can I create an actual visible light spectrum?

For example, let's say we use a phenolphthalein type of reagent such that, as the pH of a sample changes, the absorbance of particular colors changes accordingly and so the camera and image processing ultimately yields some values that can be calibrated to different pH values.

Would it be possible to use a conventional Canon DSLR camera to do some kind of vaguely meaningful spectrum analysis with GNU Octave?

Yes. In fact, you don't need it to be a Canon DSLR, you can even do this with a web cam. The main idea is close to the real thing, you use a diffraction grating to separate the light components and a sensing element (e.g. a web cam) to "read" the colours. Spatial displacement corresponds to wavelength and there is also the dimension of the "strength" of a spectral line. In this way, you do away with the transformation between recorded colour and actual wavelength (more on this later).

For more information, please see this link and this link. If you want really fine control over your diffraction grating, you can use something like lightscribe to create a diffraction grating pattern at your exact specifications. There used to be some alternatives out there to create simple CD ISOs that when burned on a CD surface they would actually create any sort of pattern on them, you just had to supply an image, but at the moment I cannot recall any of them.

If you tried to do it without a grating pattern, you would need full knowledge of the camera data OR you would need to calibrate the camera yourself. Every sensor has a spectral response to light which you need to know in advance to convert some 45,62,186 to a wavelength. The spectral response depends on the image processing AFTER the sensor (and until the digital file) and of course the incident light.

So, to do spectroscopy directly on the captured images, you would need an LED source whose current is stabilised (to produce constant quality light), you would need to measure the "temperature" of the light produced, set the camera to manual, select an exposure setting and then use the sensor's spectral response for these particular settings to obtain the correspondence between RGB value and wavelength. Otherwise, you need to calibrate the camera yourself. This still involves the constant current light source but also light sources whose intensity and wavelength you have characterised very well (or patterns). It is doable, but not exactly straightforward.

  How can I create an actual visible light spectrum?

You can use a diffraction grating as described above. You can also use one of the "old-school" transparent rulers. Their divisions are actually "etched" on the ruler so when you shine a thin strip of light on them at a very shallow angle, they act as diffraction gratings. You can actually measure the wavelength of a laser pointer very accurately in this way.

  For example, let's say we use a phenolphthalein type of reagent such that, as the pH of a sample changes, the absorbance of particular colors changes accordingly and so the camera and image processing ultimately yields some values that can be calibrated to different pH values.

Why? This is already "implemented" and calibrated and standardised. Your local chemist probably stocks uranalysis kits. One of the strips in that kit measures pH. BUT! To get an accurate measurement off of the strip, you still need to make sure that you are photographing it under very well controlled light conditions (i.e. Current control light source, or flash photography but under the same conditions every time). Your image processing part could be matching the actual colour of the pH strip, to the "reference" colour chart that indicates which colour corresponds to which pH value. A very simple program for any platform (i.e. Octave, or Scilab).

Hope this helps.


Dynamic range compression - Need pseudo algorithm for normalizing a signal I have a signal that is NOT an audio signal (but looks like an audio signal - probably you can even "play" it). I would like to amplify the signal in regions where its amplitude is low. In very rare occasions the signal amplitude may be too high so I want to decrease it a bit. At the beginning and at the end of the 'track' I have some noise.

I think what what I want is called "Dynamic range compression". Is this correct? How I do this programmatically? I need some pseudo-code.

from your posted waveform, i am assuming that this is a unipolar signal.  that is

$$ x[n] \ge 0 $$

in audio, it would be the same, except that we would be working on $|x[n]|$ instead.

so first you want a sliding maximum of your signal, where the window length is $L$.

$$ x_1[n] = \max_{0 \le i < L} \Big| x[n-i] \Big|  $$

since your input signal $x[n]$ appears to be unipolar, you can leave off the absolute value operation.

then you want to apply a low-pass filter (LPF) to that sliding max.  you want the gain of the LPF at DC (0 Hz) to be 1 (or 0 dB gain).  a simple first-order LPF is

$$ \begin{align}
 x_2[n] &= (1-p) \cdot x_1[n] + p \cdot x_2[n-1] \\
\\
        &= x_1[n] + p \cdot (x_2[n-1] - x_1[n]) \\
\\
        &= x_2[n-1] + (1-p) \cdot (x_1[n] - x_2[n-1]) \\
\end{align} $$

$p$ is the pole value of the LPF and

$$ 0 < 1-p \ll 1 $$

so

$x_2[n]$ will be an envelope for your input signal $x[n]$ and will be delayed by about half of the max window length plus about 4 times the "time constant" of the LPF:

$$d = \frac{L}{2} - \frac{4}{\log(p)} \quad \quad \text{samples}$$.

so, to normalize, you want to invert (compute the reciprocal of) the envelope $x_2[n]$ and multiply your input signal by that inverted envelope.  the inverted envelope is

$$ x_3[n] = \frac{A}{x_2[n] + \epsilon} $$

$A$ is the normalized amplitude you want (it can be $A=1$ or $A=$ any other positive number that you like).  $\epsilon$ is a tiny number, much smaller than most of your non-zero $x[n]$ that you need to add to the denominator to keep from dividing by zero (which is a bad thing).

but you should line up the signal and the delayed inverted envelope, so your normalized output is also delayed by $d$ samples (defined above):

$$ y[n] = x_3[n] \cdot x[n-d] $$

the most computationally expensive operation is the sliding max.  we were just talking about this sliding max in the music-dsp mailing list.  i will dig up C code for it and post that in a following answer.


Why do we follow Shannon's Theorem in audio sampling? My professor told me that in sampling a signal, is not usual to respect the Shannon theorem because the interpolation formula contains sines and fractions which are difficult to calculate with a microprocessor.

Why is audio sampled at 44.1-48 KHz (Nyquist frequency of audio signals)? Is there a dedicated  hardware for accomplishing this? I can sample sound with my PC, so does this mean that I have this hardware in my PC?

What your professor probably meant is that in most audio converters a small amount of aliasing is permitted (for example if you put 24.1khz into a 48khz converter you will see an output component at 23.9kHz that is not very heavily attenuated). The reason for this is that manufacturers like to use so- called half-band filters in the decimation chain because this saves a factor of 2 in the number of multiplications. Halfband filters impose a response symmetry around FS_out/2 so the attenuation at Nyquist is only 6db.  Since this is above the audio band it is not audible. However if you attempt to use this shortcut at lower sample rates like 8khz, it can be quite audible so it is not a recommended technique for lower sample-rates.

Bob


Simulating analog filter using convolution or converting with fft My task is simple; I want to simulate analog low-pass filtering of an input signal, using Python. Note that the input signal is an array of values, not an analytical function.

My first question is if it simply is possible? If so, where is the line drawn between making a digital (FIR/IIR) filter and simulating an analog filter?

My second question is how I, in python, can make two equivalent simulations, using convolution in one, and converting the input signal to frequency domain in the other. Convolution is used in the link below, and it would help me a lot to see how exactly the same thing would be done by instead using fft (even though I think it's using an FIR filter, not simulating an analog filter).

http://glowingpython.blogspot.nl/2012/02/convolution-with-numpy.html

Thanks a lot in advance!

What's not clear to me is what the fundamental difference (if any) is between simulating an analog filter and making a digital filter.

Either way, these functions will produce "ba" transfer function outputs, but the b and a are totally different.

For a 2nd-order filter, for instance, b = [b0, b1, b2] and a = [a0, a1, a2].  These are the coefficients of the transfer function.

For an analog filter, this represents a transfer function like this:

$$H(s) = \frac{b_2 s^2 + b_1 s + b_0}{a_2 s^2 + a_1 s + a_0}$$

For a digital filter, the transfer function is:

$$H(z)=\frac{b_0+b_1z^{-1}+b_2z^{-2}} {a_0+a_1z^{-1}+a_2z^{-2} }$$

  Note that the input signal is an array of values

Then you are doing digital filtering.  You can simulate analog filtering this way, but it's only going to be a simulation.

All the Python IIR filter functions you're talking about, when outputting a digital filter, design it as an approximation of the analog filter using an analog filter transfer function prototype and then converting it to digital using the bilinear transform and frequency warping.  Yes, they can output analog transfer functions, but I think that's only useful for building actually analog electronic filters, or for plotting the frequency response, impulse response, etc.

The Bessel in particular is a bad approximation as you get near fs/4, because the important property of a Bessel filter is group delay, not amplitude, and phase/group delay is not preserved by the bilinear transform (red vertical line is fs/4):

More plots of digital Bessel filter

You can make the simulation closer to what an analog filter would produce by putting a large margin between the highest frequency in your signal and the sampling frequency, called oversampling.  Assuming your "array of values" is correctly bandlimited, then you can upsample them and then apply the digital filter using lfilter and the ba digital filter coefficients you found earlier.  If not, you need to bandlimit the analog signal and then sample it at a significantly higher rate.  How did you get your sampled signal?  If you are generating it digitally, it is probably not bandlimited, unless you specifically thought of bandlimiting while generating it.

  for example a 1st-order Bessel filter

All 1st-order filters are the same, the "1st-order Bessel" is the same thing as a "1st-order Butterworth" or anything else.


Estimate/Compute parameters of image dataset I have a dataset of images taken by Kinect (RGB and D)

How can I extract the following information only from the images:

Frame rate
Camera height

The dataset containes images of person (s) moving in the scene like here. So I do have the bounding box or/and the skeleton of the person.

How can I estimate them, is there any way to do it?

Height is easy. You know that the camera is at $(0,0,0)$. Just fit a plane to the ground pixels and measure the orthogonal distance from the origin to this plane.

Frame rate is more of an issue. It is correlated with the speed of the user moving the camera. I suppose an approach would be to actually take the Kinect camera, move it like a normal person does and record the frame rate. You could then compute the optical flow in your dataset and record the magnitudes of flow vectors. You could then train a classifier (e.g. random forest) to learn to regress the frame rate given the flow vectors. Assuming that the camera motion is very similar, you could then put your dataset into this classifier and regress the frame rates for the whole trajectory. The average over all sets of frames should give you, kind of, the frame rate.


Modeling channel response as complex vs splitting I and Q If I have a communications channel where a digital QAM signal is transmitted, I can either model the channel's impulse response as a complex $h(t)$, or as 2 separate impulse responses $p(t)$ and $q(t)$ for the I and Q channels, respectively. Let $s(t)$ be the transmitted signal, with $s_I(t)$ and $s_Q(t)$ denoting its I and Q components.

For the complex case, the filtered signal would be $$[s_I*\operatorname{Re}(h)-s_Q*\operatorname{Im}(h)]+j[s_I*\operatorname{Im}(h)+s_Q*\operatorname{Re}(h)]$$

For the 2 separate case, it would be $$[s_I*p]+j[s_Q*q]$$

Clearly, they look rather different. Is either one a "better" model than the other? Better in the sense of, for example, using an adaptive filter to equalize the channel. We can train for a complex response with the cross coupling of the former, or run the channels separately. If not, can we talk about some sort of mathematical equivalence between them? Is there a mathematical mapping?

After giving it more thought, I see my original confusion. We need to be careful about what we mean by the channel's impulse response. Normally for an equalizer for this application, we are not modeling the channel alone, but rather also including the frequency mixing up and down, before and after, transmission through the channel. When we include the mixers in this way, we get the cross coupling as I mentioned for the first case. So if $h(t)$ is the response of the channel alone, then $g(t)=e^{-j\omega_ct}h(t)$ is the response of the entire system from baseband to baseband. The second case with independent responses would be for applications where the channels are not mixed together at any stage and are transmitted and received independently.


FIR design for arbitrary magnitude and phase: Why can I separate real and imaginary parts like this? In this question about the design of a FIR filter with arbitrary magnitude and phase specifications user robert bristow-johnson suggested to split the desired complex frequency response $H$ into its real and imaginary parts $H_R$ and $H_I$, then use firpm or firls to design two linear phase filters for both, setting the ftypeoption to hilbert for the imaginary part. Apparently, the sum of the two resulting FIR filters would yield a filter that meets the specification.

I'm trying to understand why this works. Here's my attempt:
$$H_R = |H|\cos(\arg\{H\})$$
$$H_I = |H|\sin(\arg\{H\})$$
Now we design two FIR filters with firpm or firls, setting the ftypeoption to hilbert for the imaginary part. To my understanding, the frequency responses of the two resulting filters $\tilde h_R$ and $\tilde h_I$ will then look like this:
$$\tilde H_R = \Re\{H\}\ e^{jp(k)}$$
$$\tilde H_I = \Im\{H\}\ e^{jp(k)+j\frac\pi2}$$
Where $p(k)$ is the linear phase term and the second filter has the additional phase shift due to the hilbert option that makes it a type III or IV FIR filter. The sum of the coefficients $\tilde h[k] = \tilde h_R[k] + \tilde h_I[k]$ would then be a FIR with the frequency response
$$\tilde H = \tilde H_R + \tilde H_I$$
because adding the coefficients is basically the same as adding the filter output signals for the same filter input signal.

\begin{align}\tilde H =& \Re\{H\}\cos(p(\omega)) + \Im\{H\}\cos(p(\omega) + \frac\pi2) \\&+ j\left[\Re\{H\}\sin(p(\omega)) + \Im\{H\}\sin(p(\omega) + \frac\pi2)\right]\\
=& \Re\{H\}\cos(p(\omega)) - \Im\{H\}\sin(p(\omega)) \\&+ j\left[\Re\{H\}\sin(p(\omega)) - \Im\{H\}\cos(p(\omega))\right]
\end{align}

This looks nothing like the response $H$ that I actually want. And I don't see how this could work. Is the original approach wrong or am I?

Neither of you is wrong. Your final expression is just the same as $H(\omega)\cdot e^{-j\omega(N-1)/2}$, which is (an approximation of) the desired complex frequency response with an additional linear phase term that depends on the chosen filter length $N$. This means that you need to subtract that linear phase term from the desired phase response, because it will be added in the design process.

Let the desired complex frequency response be given by

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

Now you have to compute a new complex response by subtracting the linear phase term $-\omega (N-1)/2$ from the phase of $(1)$:

$$\tilde{D}(\omega)=|D(\omega)|\cdot e^{j[\phi(\omega)+\omega(N-1)/2]}\tag{2}$$

You have to use the real and imaginary parts of $\tilde{D}(\omega)$ to design the filter.

What is meant by the option 'hilbert', is that you tell the design routine that the desired response (the imaginary part of $\tilde{D}(\omega)$) should actually be multiplied by $j$ and it is an odd function of frequency, because otherwise it will be interpreted as a purely real and even function. This will result in an asymmetric impulse response. In principle this can be done, but I don't know if it can actually be done for arbitrary prescribed responses with the current versions of Matlab's filter design routines.

Note that such an approach is actually not really necessary for the least squares design of FIR filters, because solving a complex least squares problem (for non-linear phase designs) is just as simple as solving a real least squares problem (for linear phase designs). As far as I know it's just not implemented in Matlab's Signal Processing Toolbox. I've written a few Matlab/Octave functions that solve the complex least squares approximation problem for FIR filter design: lslevin.m for complex least squares design of real-valued filters, and cfirls.m for filters with complex coefficients.

For equi-ripple designs, the complex (non-linear phase) approximation problem is indeed much harder than the real (linear phase) problem. For the latter we have the Remez exchange algorithm, which was adapted to the linear phase FIR filter design problem by Parks and McClellan and which is implemented in Matlab (firpm.m). There are algorithms for solving the complex equi-ripple design problem, but they are much less efficient than the Parks-McClellan algorithm. One of them is implemented in Matlab (cfirpm.m).

Note that for the least squares approximation, solving two real-valued approximation problems for the real and imaginary parts gives the same solution as solving one complex least squares approximation problem. This is easily seen as follows. The complex approximation error is given by

$$E(\omega)=D(\omega)-H(\omega)\tag{3}$$

where $D(\omega)$ and $H(\omega)$ are the desired and actual complex frequency responses. The complex least squares solution minimizes the integral over $|E(\omega)|^2$, whereas the individual real-valued least squares solutions minimize the integrals over $E_R^2(\omega)$ and $E_I^2(\omega)$, where $E_R(\omega)$ and $E_I(\omega)$ are the real and imaginary parts of $E(\omega)$, respectively. Since

$$|E(\omega)|^2=E_R^2(\omega)+E_I^2(\omega)\tag{4}$$

minimizing the integrals over $E_R^2(\omega)$ and $E_I^2(\omega)$ is equivalent to minimizing the integral over $|E(\omega)|^2$.

For the complex Chebyshev approximation problem the situation is different. The quantity that is minimized is the maximum of the magnitude of the complex error over the bands of interest $\max_{\omega}|E(\omega)|$, but since

$$\max_{\omega}|E(\omega)|\neq \max_{\omega}|E_R(\omega)|+\max_{\omega}|E_I(\omega)|\tag{5}$$

independent Chebyshev approximation of the real and imaginary parts is not equivalent to complex Chebyshev approximation. If for a given complex approximation problem the maximum error of the optimal solution equals $\delta$, then in the worst case the errors of both real-valued approximations will also be equal to $\delta$. This means that in the worst case, the maximum complex error of the solution obtained from independent optimization of the real and imaginary parts becomes

$$\max_{\omega}|E(\omega)|=\max_{\omega}\sqrt{E_R^2(\omega)+E_I^2(\omega)}\le\sqrt{\max_{\omega}E_R^2(\omega)+\max_{\omega}E_I^2(\omega)}\le\sqrt{2}\delta\tag{6}$$

This is equivalent to a maximum of $3$ dB increase in error compared to the optimal solution.

I'll show an example design to illustrate a few points that I made. I designed a bandpass filter with a desired phase response that is linear but with a smaller group delay than an FIR filter with perfectly linear phase. The filter length is $N=61$, so the delay of a linear phase FIR filter would be $(N-1)/2=30$ samples. I chose the desired delay to be $20$ samples instead. I used a complex least squares and a complex Chebyshev criterion to approximate the given complex frequency response. For each criterion I designed the filter in two different ways. First, directly in the complex domain, and second by using the real and imaginary parts of the desired frequency response to design two linear phase filters, which, after combining them, should approximate the original complex frequency response. The figure below shows the design results. As expected, for the least squares design, both solutions are equivalent (up to numerical accuracy) and only one of them is shown in the plots (in blue). For the Chebyshev criterion both solutions are slightly different, as predicted. The solution from splitting the response into real and imaginary parts (in red) has a larger error by about $1.2$dB (as can be seen at the first sidelobe in the stopband). In the passband all three solutions are very similar. Note that all three filters approximate the desired group delay value of $20$ samples relatively well over a large portion of the passband.


What is the x and y axis showing for this matplotlib.pyplot histogram of a .wav file? The python code below generates a histogram from a .wav file, but what exactly does the histogram show? Is the x axis amplitude binned by sampling rate? Is the y axis a count of how many samples are in each amplitude bin? And how is amplitude calculated?
import numpy as np
import matplotlib.pyplot as plt

Fs, data = wavfile.read('audio file')

plt.hist(data, bins='auto')
plt.show()

It seems we have a confusion with the meaning of the term frequency in different contexts.
In signal processing we tend to use the term frequency as a measure of the rate in the a signal completes a cycle.
In statistic the frequency of an event is number of times an event occurs.
The histogram plots frequencies in the statistic sense. Each bin corresponds to an event defined as "the sample has value in the interval $b_k < x_i < b_{k+1}$", where $b_k$ corresponds to the position of the vertical edges of each bar in your histogram.
So what you are plotting in the histogram is associated with how many times the sample amplitude lied in each of those ranges.


Arbitrary shaped RIR generator I wish to simulate some sound sources convolved with a Room Impulse Response (RIR) and impinging on a spherical microphone array. I could find a couple of generators online (such as this one), though all of these are concerned with rectangular rooms.

Is there a way to generate such RIR for an arbitrary shape or to access a bigger variety of shapes for the room?
If there are simulators for such RIR, are there online implementations? (I could not find any)
Is there a common method to convert a rectangular RIR to an RIR for a different shaped room?

Arbitrary shapes is A LOT more difficult than rectangular rooms. You either need to do image method with visibility checking or some sort of ray, particle or cone tracing algorithm.

There are few commercial programs that do this (Ease, Odeon, Bose Modeller, etc. ) and I doubt you'll find a free one. In addition the much more complicated algorithms, these need a way to enter or import complex room geometries, and extensive data bases on material (absorption, reflection, diffusion) and source polar patterns. These are all required to get a reasonably realistic impulse response.

Most of the simpler online programs make so many simplifying assumptions, that the results tend to be very unrealistic and not useful for many purposes.


Correlation between real and imaginary parts of a Fourier transform of zero mean Gaussian From this previous post, the real and imaginary parts of the Fourier transform of a zero mean Gaussian are uncorrelated (and i.i.d. Gaussians)

This some how seems counter-intuitive. It seems if the real part has a large intensity, then the imaginary part should be small. I read through the proof linked from the previous post, but I don't have a good intuitive feel from this proof.

Can you provide an intuitive explanation of why this is true?

Recall the following properties of the Fourier transform:

If $x(t)$ is an even function, then its Fourier transform $X(\omega)$ is purely real.
If $x(t)$ is an odd function, then its Fourier transform $X(\omega)$ is purely imaginary.

Thus, we can think of the real and imaginary parts of the Fourier transform of a zero-mean Gaussian random process $x(t)$ as the Fourier transforms of two separate inputs $x_e(t)$ and $x_o(t)$: the components of the process that have even and odd symmetry, respectively. We split the process into these components as follows:

$$
x_e(t) = \frac{x(t)+x(-t)}{2}
$$

$$
x_o(t) = \frac{x(t)-x(-t)}{2}
$$

Note that $x(t) = x_e(t) + x_o(t)$. Moving through this step by step,

Your observation was that $\text{Re}\{X(\omega)\}$ and $\text{Im}\{X(\omega)\}$ (the real and imaginary parts of the process's DFT) are uncorrelated.
Using the Fourier transform properties mentioned above, we can deduce that $\text{Re}\{X(\omega)\} = X_e(\omega)$ and $\text{Im}\{X(\omega)\} = X_o(\omega)$; the real and imaginary parts of $X(\omega)$ are none other than the Fourier transforms of $x_e(t)$ and $x_o(t)$, respectively.
Therefore, your observation is equivalent to saying that $X_e(\omega)$ and $X_o(\omega)$ are uncorrelated.
Since the Fourier transform is a one-to-one mapping between the time and frequency domains, I posit that the lack of correlation between $X_e(\omega)$ and $X_o(\omega)$ would imply a lack of correlation between $x_e(t)$ and $x_o(t)$ as well.

What is the correlation between $x_e(t)$ and $x_o(t)$? Simple:

$$
\begin{align}
\mathbb{E}(x_e(t)x_o(t)) &= \mathbb{E}\left(\left(\frac{x(t)+x(-t)}{2}\right)\left(\frac{x(t)-x(-t)}{2}\right)\right) \\
&= \mathbb{E}\left(\frac{1}{4}\left(x^2(t) - x^2(-t)\right)\right) \\
&= \frac{1}{4} \left(\mathbb{E}(x^2(t)) - \mathbb{E}(x^2(-t))\right) \\
&= \frac{1}{4} \left(\sigma^2 - \sigma^2\right) \\
&= 0
\end{align}
$$

As expected, the even and odd components are uncorrelated. So, to summarize, I would say the following:

The real and imaginary components of a Fourier transform correspond to the individual Fourier transforms even and odd components of the input function.
For a zero-mean Gaussian random process, these even and odd components are uncorrelated.
Therefore, their Fourier transforms (the real and imaginary components that you asked about) are also uncorrelated.

Edit: To address your followup:

If $x(t)$ is Gaussian, then its even and odd components $x_e(t)$ and $x_o(t)$ are as well, due to the property that any weighted sum of Gaussian random variables is also Gaussian.
If $x_e(t)$ and $x_o(t)$ are Gaussian random processes, then their Fourier transforms $X_e(\omega)$ and $X_o(\omega$) are as well. This follows from the same property as the previous statement; if you look at the transform, you're computing a weighted sum of a bunch of Gaussian random variables.
If $X_e(\omega)$ and $X_o(\omega)$ are Gaussian, and they are uncorrelated with one another (as described above), then they are also independent. This is a property of the Gaussian distribution.


Why does the impulse response determine the transfer function of a system? why does it describe the transfer function.. How come? Especially for LTI systems.
I thinking about the theory about how come an impulse input can provide information about a complete system.. As it doesn't operate fully of all configuration.

Because every signal can be decomposed into a linear series of scaled impulses that are shifted in time. Thus applying a linear time-invariant system on each impulse in the series and then summing the results again will give the same result as applying it on the whole signal. Therefore, one only needs to know how the system responds to an impulse to be able to calculate its response to an arbitrary signal and all properties of the system.
This DSP guide chapter has some nice illustrations on the subject.
This DSP guide chapter (only in PDF) shows how from the impulse response as written in the form of a difference equation the transfer function can be derived using the z-transform. For FIR filters the $b_i$ coefficients will be zero in equation 33-2.


How to convert my transfer function to the frequency domain I have the following transfer function in the time domain:

$h\left(t\right) = B x\left(t\right) + a x^2\left(t\right) + b x^3\left(t\right)$.

My simulation is in the frequency domain, so I would like to have a transfer function of the form $H\left(\omega\right)$.

I didn't find any suitable Fourier/Laplace identities so I tried substituting as follows $\theta_1\left(t\right) = x\left(t\right)$, $\theta_2\left(t\right) = x^2\left(t\right)$ and $\theta_3\left(t\right) = x^3\left(t\right)$. The resultant frequency domain function

$H\left(\omega\right) = B \Theta_1 \left(\omega\right) + a  \Theta_2 \left(\omega\right) + b  \Theta_3 \left(\omega\right)$,

wasn't of much use to me because I only have $\Theta_1\left(\omega\right)$ and I don't know how  $\Theta_2$ and $\Theta_3$ are related to $\Theta_1$.

I would like my transfer function to depend only on $\Theta_1$, which is the variable that I know.

How can I go about solving this problem?

Since this is not an LTI system, we cannot have $H(\omega) = K\Theta_1(\omega)$, where $K$ is a scalar complex number. But as OP mentioned in the comment he has knowledge of $\Theta_1(\omega)$, we can compute $\Theta_2(\omega)$ and $\Theta_3(\omega)$. As $x^2(t)=x(t)\times x(t)$,
$$
\Theta_2(\omega)= \Theta_1(\omega)*\Theta_1(\omega) =\int_{-\infty}^{+\infty} \Theta_1(\alpha)\Theta_1(\omega-\alpha)d\alpha
$$
where $*$ is the convolution operation. Similarly, $x^3(t)=x(t)\times x(t)\times x(t)$
$$
\Theta_3(\omega)= \Theta_1(\omega)*\Theta_1(\omega)*\Theta_1(\omega)=\Theta_1(\omega)*\Theta_2(\omega)
$$
So $H(\omega)=B\Theta_1(\omega)+a\Theta_2(\omega)+b\Theta_3(\omega)$, where $\Theta_2$ and $\Theta_3$ are computed from $\Theta_1$ above.


How do we know the receiver delay in an OFDM transceiver system? I was simulating a simple OFDM transceiver system without a channel (just the transmitter and the receiver directly connected). I saw that if I use the "Find Delay" block in Simulink between the transmitted and received bits, I find that it keeps changing, because of which it becomes difficult to decide a specific receiver delay to calculate the BER. This however shows a constant value of delay in just a simple BPSK transceiver system without the OFDM. Can anyone please tell me what is actually happening and what are the things to consider while calculating the BER of an OFDM system?

In basically all simulations I've done, you just send known (but, depending on what you want to simulate, seemingly random) data. Then you correlate the received bits with the transmitted ones, and look for a peak.

That's your end-to-end bit delay, which you don't really care about, but it's how you know which bit you compare to which.

If you data has no significant autocorrelation, and assuming your transceiver basically works with a BER much smaller than 0.5, this is pretty reliable.


New to audio processing, why is my data so lopsided when I try to break into frequency ranges? I have a wav file audio data, I broke it up into 1024-length windows (no overlap), and performed fft on each one.

If I visualize this data it actually looks pretty good, but the problem is that the data is really lopsided. It seems like buckets 1-3 get good activity, but they're generally much larger than buckets 4-8, so when I visualize the data I have to have a weird conditional multiplier on the higher frequency buckets so I see some activity.

So then, what is the proper way to break my fft into frequency buckets? A simple explanation would be best. Thank you!

Most of the energy in an audio/speech signal is almost always found in the lower bands (roughly below 1kHz) so the lopsided shape that you're observing  is not surprising.

Let me just add some more information regarding fft frequency bins. These bins can be organized into larger buckets if you need this kind of representation.

For an FFT length of 1024, you should end up with 513 frequency bins (N/2+1). The first bin corresponds to DC (0 Hz) and is usually ignored and N/2+1 is the theoretical Nyquist frequency (also ignored).

The bandwidth of a frequency bin is defined as

BW_bin =  Sampling_rate/FFT length;

Note that while the sampling rate isn't necessary to compute the FFT, it is needed to calculate the bandwidth (frequency resolution).

To get the activity (magnitudes in dB), you can use the following equation (one way to compute the magnitudes)

Mag[i] = 10*log*(sqrt(2*(Real[i]^2+Img[i]^2/fftNorm));//iterate through bins [1-511]

fftNorm depends on the kind of the window function used (https://en.wikipedia.org/wiki/Window_function) and is simply N (fft length) in the case of a rectangular window (no window).
The factor 2 in the equation accounts for the upper (discarded) half of the FFT.

For visualization purposes, you can now easily combine several bins into a group of bins coressponding to a particular frequncy range.


Despeckling Ultrasound images using Diffusion Filtering I am working on Despeckling Ultrasound images. I have been trying out some basic filters and their hybrid combination.. Now I would like to try Diffusion Filtering techniques to despeckle the ultrasound images. I don't know anything about the diffusion filtering techniques (Isotropic diffusion or Perona Malik Anisotropic Diffusion Filter). I would like to learn from the basic. Please suggest some books to learn about Diffusion techniques.. Please keep in mind that I am not an advanced learner.

I did google and found this book (Despeckle Filtering Algorithms and Software for Ultrasound Imaging By Christos P. Loizou, Constantinos S. Pattichis, Costas Pattichis) very useful. But I couldn't able to buy it. If anyone have the PDF of this book (even 1 or 2 chapters), please provide me the link. Also provide some useful links to learn these techniques

Reference: Anisotropic Diﬀusion in Image Processing by Joachim Weickert

No matter where would you learn the theoretical background of the diffusion techniques, I would like to apply the Matlab implementation on the 2nd formal definition of anisotropic diffusion, it is not difficult. I used the ultrasound B-scan image you used before.

I=imread('ultrasound.png');
I=rgb2gray(I);
subplot(1,2,1)
imshow(I);
[rows, cols]=size(I);
diff = double(I); % original image
lambda = 0.25;
niter = 10;
Co = 20;
for i = 1:10  % iterations

  % Construct diffl which is the same as diff but
  % has an extra padding of zeros around it.
  diffl = zeros(rows+2, cols+2);
  diffl(2:rows+1, 2:cols+1) = diff;

  % North, South, East and West differences
  deltaN = diffl(1:rows,2:cols+1)   - diff;
  deltaS = diffl(3:rows+2,2:cols+1) - diff;
  deltaE = diffl(2:rows+1,3:cols+2) - diff;
  deltaW = diffl(2:rows+1,1:cols)   - diff;

  cN = 1./(1 + (deltaN/Co).^2);
  cS = 1./(1 + (deltaS/Co).^2);
  cE = 1./(1 + (deltaE/Co).^2);
  cW = 1./(1 + (deltaW/Co).^2);

  diff = diff + lambda*(cN.*deltaN + cS.*deltaS + cE.*deltaE + cW.*deltaW);

end
subplot(1,2,2),imshow(uint8(diff))

You can see the edge preserving effect of diffusion filter (right) compared with the speckled noisy image (left).


Filter coefficients for colored noise Voss Algorithm If I use a generic filter for generating colored noise like pink,brown,white then how do I modify this statement and how do I know what are the coefficients to be used in AR model for different noise.Signal is white gaussian noise and y gives different noise based on the value of a.

Here, $a$ is a filter coefficient for a simple single-pole recursive low-pass filter. If you look at the documentation for MATLAB's filter() function, you'll see this corresponds to the difference equation $y(n) = x(n) - a y(n-1)$`. Let's analyze the filter in the $z$-domain to get its frequency response (just remembering that $z^{-1}$ corresponds to a single sample delay:

$$y(n) + a y(n-1) = x(n)$$
$$Y(z)*(1 + a z^{-1}) = X(z)$$
$$\frac{Y(z)}{X(z)} = \frac{1}{1 + a z^{-1}} = H(z)$$

That's the $z$-domain transfer function. We can substitute $z = e^{j\omega}$ where $\omega$ is the normalized frequency, where 0 is DC and $\pi$ corresponds to the Nyquist frequency of whatever sampling rate you're at. Now you should realize that if you use the same digital coefficients at different sampling rates, you will get filtered noise that varies depending on the sampling rate. If your sampling rate is $F_s$, then $\omega$ corresponds to the frequency $\dfrac{\omega F_s}{2 \pi}$. Let's find the frequency response:

$$G^2(\omega) = |H(e^{j\omega})|^2 = \left| \frac{1}{1 + a e^{-j\omega}} \right|^2$$

$$G^2(\omega) = \frac{1}{|1 + a(\cos(-\omega) + j \sin(-\omega))|^2}$$

$$G^2(\omega) = \frac{1}{|(1 + a \cos(\omega)) - j \sin(\omega))|^2}$$

$$G^2(\omega) = \frac{1}{(1 + a \cos(\omega))^2 + \sin^2(\omega)}$$

$$G^2(\omega) = \frac{1}{1 + 2 a \cos(\omega) + a^2 \cos^2(\omega) + \sin^2(\omega)}$$

$$G^2(\omega) = \frac{1}{1 + 2 a \cos(\omega) + a^2 \cos^2(\omega) + (1 - \cos^2(\omega))}$$

$$G^2(\omega) = \frac{1}{2 + 2 a \cos(\omega) + (a^2 - 1) \cos^2(\omega)}$$

$$G(\omega) = \frac{1}{\sqrt{2 + 2 a \cos(\omega) + (a^2 - 1) \cos^2(\omega)}}$$

Let's plot the frequency response with $a = -0.9$: WolframAlpha plot. We can also take 20 Log10() to get the frequency response in decibels: WolframAlpha plot .

Now we can find out some things about this filter. For example, its gain at DC is greater than unity! $G(0) = \dfrac{1}{\sqrt{1 + 2a + a^2}}$. When $a = -0.9$, the filter increases the DC gain by a factor of 10, or 20 dB. Often filters are characterized in terms of their -3 dB point, which is where their magnitude decreases by $\sqrt{2}$ from their peak. Let's solve for this point:

$$\frac{1} { \sqrt{2} \sqrt{1 + 2a + a^2} } = \frac{1}{\sqrt{2 + 2 a \cos(\omega) + (a^2 - 1) \cos^2(\omega)}}$$

$$\sqrt{2} \sqrt{1 + 2a + a^2} = \sqrt(2 + 2 a \cos(\omega) + (a^2 - 1) \cos^2(\omega))$$

$$2 + 4a + 2 a^2 = 2 + 2 a \cos(\omega) + (a^2 - 1) \cos^2(\omega)$$

Now treat this as a quadratic in $\cos(\omega)$, let's substitute $x = \cos(\omega)$ and use the quadratic formula:

$$(a^2 - 1) x + 2ax + (-4a - 2 a^2) = 0$$
$$x = \frac{-2a \pm \sqrt( 4a^2 - 4(a^2-1)(-4a - 2 a^2)) }{2 (a^2 - 1)}$$

Now, we want the solution corresponding to positive frequencies, and it turns out that choosing the result with the minus sign gives us that. So, we can take the inverse cosine, and get our -3 dB frequency:

$$\omega = \arccos\left( \frac{-2a - \sqrt{4a^2 - 4(a^2-1)(-4a - 2 a^2)}}{2 (a^2 - 1)} \right)$$

Let's plug it in for $a = -0.9$ and verify:

$w = 0.0958385$

$H(z) = 7.07107 = \dfrac{10}{\sqrt{2}}$

It might be useful to go the other way, to design a filter by choosing a value of a to have some desired -3 dB point:

$$2 + 4a + 2 a^2 = 2 + 2 a \cos(\omega) + (a^2 - 1) \cos^2(\omega)$$

$$2 + 4a + 2 a^2 = 2 + 2 a \cos(\omega) + a^2 cos^2(\omega) - \cos^2(\omega)$$

$$(2 - \cos^2(\omega)) a^2 + (4 - 2 \cos(\omega)) a + \cos^2(w) = 0$$

$$\frac{ -(4 - 2 \cos(\omega)) +/- \sqrt{ (4 - 2 \cos^2(\omega))^2 - 4(2 - \cos^2(\omega))(\cos^2(\omega)) } }{2 (2 - \cos^2(\omega))}$$

And again we can verify by plugging in $\omega = 0.0958385$, and we'll see that we get $a = -0.9$.

So, for example, you could use these formulas to pick values of $a$ at different sampling rates, so that you have the same frequency -3 dB point. You can use the formula for gain at DC, for example, to normalize the filter so its gain at DC is unity. You can take the formula for gain, convert it to decibels, and see that it approaches roughly a -6 dB/oct slope, which is why it's being used to shape noise like this.


Is noise figure dependent on input noise power? I was reading about Noise Figure on Wikipedia (https://en.wikipedia.org/wiki/Noise_figure). There, I saw the following definition:

  The noise factor is thus the ratio of actual output noise to that which would remain if the device itself did not introduce noise, or the ratio of input SNR to output SNR.

where noise figure is $10\log_{10} \text{(Noise Factor)}$.

That arises a question in me.

  If Noise Figure is given for an amplifier as a constant (for example Noise-Figure=$5.5$dB), since amplifier noise is additive (i.e. its power adds up to the input noise power), then the amplifier noise must be dependent on input noise power, which contradicts with my intuition.

As an example, is an amplifier has gain of $1$ (i.e. no amplification, no attenuation) with a NF of $5.5$dB, it preserves the input signal power, in a sense that $$P_{\text{out}}=G\cdot P_{\text{in}}=P_{\text{in}}$$and $$\sigma^2_{\text{noise,out}}{=G\cdot\sigma^2_{\text{noise,in}}+\sigma^2_{\text{amplifier}}\\=\sigma^2_{\text{noise,in}}+\sigma^2_{\text{amplifier}}}$$
hence $$\text{Noise Factor}={\sigma^2_{\text{noise,out}}\over \sigma^2_{\text{noise,in}}}=1+{\sigma^2_{\text{amplifier}}\over \sigma^2_{\text{noise,in}}}=10^{0.55}=3.55$$and $$\sigma^2_{\text{amplifier}}=2.55\sigma^2_{\text{noise,in}}\qquad\qquad(!!)$$

However, we know that Noise Figure is highly used in practice. Where is my mistake in understanding it?

To answer your question, Noise Figure is the amount of noise that would be added referenced to the input of the amplifier assuming no prior noise amplification has been provided in a properly matched system. This equates to a change in the SNR of the output relative to the SNR of the input, and the input noise level with no other gain would be typically referred to as the "thermal noise floor".

So given that, and to your example, if the gain of the amplifier was 1 (0 dB) and the noise figure was 5.5 dB, then in this case the signal out would be equal to the signal in (gain = 1) but the noise out would be 5.5 dB higher than the signal in. (and hence if you follow SNR out vs SNR in you see that this would hold).

To further demonstrate, consider if the amplifier had a gain of 10 dB with the same 5.5 dB Noise Figure. Here the signal out would be 10 dB higher than the signal in, but the noise out would be 15.5 dB higher than the noise in.

Now to further demonstrate noise figure, consider the same case with adding a noise free 10 dB gain amplifier in front of our noisy 5.5 dB NF amplifier: At the input to the noisy amp the signal will be 10 dB higher and the noise will be 10 dB higher relative to the thermal noise floor. Our noisy amp adds +5.5 dB of noise relative to the thermal noise floor, so referenced to the noisy amp input according to this noise figure model, in this case, the total noise will be:

$10Log(10^{10/10}+10^{5.5/10}) = 11.32$ dB above the thermal noise floor.

While the signal will be 10 dB above the input signal level at this point.

(Note this isn't the actual level that would be measured at the input of course but is the "input referenced" noise as measured at the output of the amplifier--as the noise is added in the amplification process).

At the output of our noisy amp, due to its +10 dB gain, the signal will be an additional 10 dB higher, as well as the noise, resulting in a final noise output level that is 21.32 dB above the thermal noise floor. The total change in signal level was +20 dB, while the total change in noise level was +21.32 dB. The change is SNR was therefore only 1.32 dB. This is the Cascaded Noise Figure showing how a prior lower noise amplifier can reduce the overall noise figure of the system.  The actual formula for this is given below, but following the details above helps give the background of how this works. The bottom line is that the noise figure is the additional noise that is added referred to the input of the device relative to the thermal noise floor. Increasing signal and noise prior can swamp out this additional added noise later.

$F_T = F_1 + \frac{F_2-1}{G_1} + \frac{F_3-1}{G_1G_2} ...$

Where

$NF = 10Log_{10}(F)$

and gain in dB = $10Log_{10}(G)$


Combining time domain and frequency domain data for BPM detection in real time I would like to write a BPM detector. Does either of the domains provide more information for this type of analysis?

What exactly is being looked for in the data to detect onset for BPM detection in either domain?

Does tracking energy values in the time domain lead into rhythm detection?

Obviously pattern recognition can be used but how can it be done fast enough in real time using c?

BPM detectors work by finding the strongest period in the onset detection function - a function of time indicating how likely it is that there is an onset at time t. Which brings the topic of onset detection...

Tracking energy/envelope in the time domain is a weak strategy for detecting rhythm/notes, because:

For some instruments, there are playing techniques which create new notes without a significant change of energy (for example moving the finger on a string of a cello with no change in bowing, or changing the fingering on a flute with continuous blowing). This has the effect of creating changes in fundamental frequency without a noticeable change in the time-domain envelope.
A lot of commercial music is so heavily compressed that the signal is almost "pumped" to maximum amplitude all the time.
While this is yet another corner case, it is worth reminding that some instruments or playing techniques have very slow attacks.

So naively looking for peaks in the time-domain (via energy) will not be very robust. You have to look for better "clues", and the frequency domain is the answer. The simplest thing that can be done in the frequency domain is to compute the difference of energy between adjacent STFT frames (a measure known as the "spectral energy flux").

A more elaborate metric, which is actually as cheap in terms of computing costs, consists in checking the deviation between the actual and expected value of the complex amplitude in a STFT bin (reference here) - this deviation would be null in the stable segment of a note (as the amplitudes are just sustaining or slowly decaying; and the phases are just keeping on rolling), and high at an onset.

All these techniques have very modest computational costs and can be implemented in realtime even on relatively cheap microcontrollers - a complete beat-tracker (onset detection + BPM detection + beat tracking) is in the 10 MIPS range.


Frequency response of FIR (MTI delay line canceller) with non constant period Short question

How to plot the magnitude of frequency response of the delay line canceller with non constant period using GNU Octave? Or more directly: how to plot magnitude of frequency response of staggered PRF (pulse repetition frequency) MTI (moving target indicator) processors?

Full question

I try to plot magnitude of frequency response of the MTI delay line canceller (some kind of FIR) which has the following structure using GNU Octave:

With tips in my related post I can plot frequency response for constant period:

close all;
clear all;

T   = 1e-3;
f   = linspace (0, 2 / T, 1000);
w   = 2 .* pi .* f;
z   = exp (-j .* w .* T);
H_z = 1 - z .^ -1;

plot (f, abs (H_z), 'r', "linewidth", 2);
hold on;
grid on;
title ("Magnitude of the frequency response of the delay line canceler.");
xlabel ("Frequency.");
ylabel ("Magnitude.");

But now I am interesting in frequency response for non constant period:

determinate values of period;
random values of period.

I try to do in the following way:

close all;
clear all;

T_1   = 0.001;
T_2   = 0.0012;
f     = linspace (0, 2 / min (T_1, T_2), 1000);
w     = 2 .* pi .* f;
z_1   = exp (-j .* w .* T_1);
z_2   = exp (-j .* w .* T_2);
H_z_1 = 1 - z_1 .^ -1;
H_z_2 = 1 - z_2 .^ -1;
k     = 1 / 2;
H_z   = k .* (H_z_1 + H_z_2);

plot (f, abs (H_z), 'r', "linewidth", 2);
hold on;
grid on;
title ("Magnitude of the frequency response of the delay line canceler for alternate values of period T1 and T2.");
xlabel ("Frequency.");
ylabel ("Magnitude.");

Another words, I mean that: common frequency response is the normalized sum of the frequency response for each value of the period.

Am I right? If yes, is this approach will be correct for random values of period?

Notes

Here is some related papers, when mentioned frequency response of staggered PRF MTI processors:

first (Pg 66 in books numeration or 81 in pdf-document numeration);
second (Pg 32).

If you mean that you want a time-varying delay then the whole system becomes time-varying, which means that there is no such thing as a frequency response in the conventional sense. Only linear time-invariant systems can be described by a frequency response. There are of course ways to describe time-varying systems, but the question is what it exactly is that you expect from such a description of the system. What you've done in your example is subtract two delayed signals from the original signal, which is in fact a time-invariant system. In sum, what you're looking for (a frequency response description of a system with a time-varying delay) does not exist.


Reduce the number of frequencies after FFT I have a long signal (1 day sampled at 500 Hz), which results in many frequencies if I calculate the amplitude spectrum (see e.g. here). However, I am only interested in the "big picture", let's say only frequencies of natural numbers (i.e. 1, 2, ..., 250 Hz).

Should I simply remove all other frequencies? Should I do something like a binning? Or is there even a better solution?

Apply a sliding window function of 500 samples length to the signal, for example a Hamming window. Choose a suitable overlap between windows e.g. 75%. Calculate the spectrum of each window and then average over all windows. Then the frequency bins lie exactly at the natural numbers as you wanted.

But beware: each bin does not only represent one exact frequency, but the energy in a frequency band, in this case  $[f-0.5,f+0.5]$Hz


FFT strong signal in middle frequency bin I've developed a radar that I'm using to determine the distance to remote objects. It uses a custom PCB with an onboard FPGA that performs the DSP algorithms. The data from it is then plotted on a host PC. This appears as a 2D histogram where the y-axis denotes the FFT frequency bins (due to the nature of the radar this is proportional to distance) and the x-axis denotes time. The plot (shown below) gives a very strong signal at a distance halfway between the antenna and max range, which I'm unable to explain.

The actual algorithms performed are: an FIR polyphase decimation filter (downsamples from 40MHz to 2MHz) which produces an output 1024 samples in length. Then I run it through a Kaiser window function with a beta of 6, followed by a 1024-point FFT, the result of which is transmitted to the host PC. For each value of t in the plot (the x-axis) the host PC averages over 30 1024-length sequences (averaged element-wise). Since all inputs to the FFT are real, the output is Hermitian symmetric and so I only plot the first 512 values of each output sequence. The strong signal you see above occurs at bins 257 and 258 (indexed from 0). I've tested the radar in an open space where it shouldn't generate any strong signals. I've simulated all of the FPGA logic, so while I can't be sure it's right (I've only formally verified parts of it), I'd be surprised if it wasn't.

What could be the cause of this? Is there some obvious aspect I'm missing? If any of this is unclear or some part of the information I've omitted is important for answering this question (e.g. the equation relating frequency to distance), please let me know and I'll include it.

Edit: more details on acquired signal

This is an FMCW radar. A frequency synthesizer generates sawtooth ramps from 5.3GHz to 5.9GHz over a duration of 1ms. This signal is simultaneously transmitted and mixed back in with the reflected signal. We then measure the difference frequency to back out the distance.

The FPGA modules are timed such that data is only acquired during the synthesizer's ramp period. First, I enable the ramp and power amplifier and (once enabled) begin acquiring data. The data is processed by the FIR filter and then passed through the kaiser window. Once the last sample passes through the FIR filter the ramp and power amplifier are disabled. The processed data (which were stored in a FIFO) are now run through the FFT and then the resulting output is dispatched in packets to the host PC via USB. I use a header sequence, tail sequence, and duplicated packets to try to avoid data corruption/loss. Once the FFT is finished, the process starts again (ramp and power amp enable, etc.).

The FIR filter should take just longer than 0.5ms to acquire all samples, so it should fall within the frequency ramp period.

If your system is powered by a switching power supply that is not synchronized to the PRF (or multiples of), then you may get reliable spurs over time as seen in your histogram.

In your case, it's been found that was indeed the problem! The 500 KHz frequency fell nicely in the middle of your 1 MHz histogram. Hopefully this will help hunt down similar issues in the future.


Implementing Gaussian random variable by using a uniform random variable I'm trying to write a C++ function that will return Gaussian random values, given their means and variances.

There is a library function rand(), which returns random numbers between 0 and RAND_MAX. RAND_MAX does not have a fixed value, but it is guarantied that it will be at least $2^{15}-1$. Its PDF is uniform.

I'm using Central Limit Theorem to transform this rand()into a Gaussian variable. What I'm exactly doing is a to call rand() for a user specified times, then add up their return values, then shift its mean to the user specified mean.

In the plotting above, I called my Gaussian random generator for $10^7$ times, and plotted frequencies of its return values. As you see, its variance is huge, since it is created by sum of a lot of other random values.

It successfully returns a Gaussian variable with a Gaussian PDF and with the specified mean value. However, the problem is its variance. I'm stuck at this point, because I don't know how to change its variance to the user specified value.

This is my code (incomplete for now; the parameter "Variance" is ignored):

template <class T>
T Random::GetGaussian(T Mean /*= 0*/, T Variance /*= 1*/)
{
    T MeanOfSum = NUM_GAUSSIAN_SUMS / static_cast<T>(2);
    T Rand = 0;
    for (uint64_t i=0; i<NUM_GAUSSIAN_SUMS; i++)
    {
        Rand += static_cast<T>(rand()) / RAND_MAX;
    }
    return Rand - (MeanOfSum - Mean);
}

Assume that NUM_GAUSSIAN_SUMS is 100, and RAND_MAX is 32767.

I want to change variance of the random variable according to the parameter of the function. My question is, how can I change variance of this random variable? How can I do it?

Your initial algorithm creates a random variable that's uniformly distributed between 0 and 1. The variance of that is 1/12. If you sum NUM_GAUSSIAN_SUMS instances of that the variance will be NUM_GAUSSIAN_SUMS/12. In order to get to a target variance, V, you need to multiply the summed random variable with sqrt(V*12/NUM_GAUSSIAN_SUMS).

As a side note, a template will work reasonable well for floats and doubles but there will be significant numerical problems with any fixed point type.


The Proper Way to Do Sinc Downsampling (DFT Downsampling) for Uniformly Sampled Discrete Signals with Finite Number of Samples Given a signal $ \left\{ x [ 0 ], x [ 1 ], ..., x [ N - 1 ] \right\} $ what would be the correct way to downsample it in the frequency domain (Sinc interpolation)?

Interpolation in Frequency (DFT Domain)
The implementation is well known. In MATLAB it will be something like:
if(numSamplesO > numSamples)
    % Upsample
    halfNSamples = numSamples / 2;
    if(mod(numSamples, 2) ~= 0) % Odd number of samples
        vXDftInt = interpFactor * [vXDft(1:ceil(halfNSamples)); zeros(numSamplesO - numSamples, 1, 'like', vXDft); vXDft((ceil(halfNSamples) + 1):numSamples)];
    else % Even number of samples -> Special Case
        vXDftInt = interpFactor * [vXDft(1:halfNSamples); vXDft(halfNSamples + 1) / 2; zeros(numSamplesO - numSamples - 1, 1, 'like', vXDft); vXDft(halfNSamples + 1) / 2; vXDft((halfNSamples + 2):numSamples)];
    end
else
    % Downsample
    halfNSamples = numSamplesO / 2;
    if(mod(numSamples, 2) ~= 0) % Odd number of samples
        vXDftInt = interpFactor * [vXDft(1:ceil(halfNSamples)); vXDft((numSamples - floor(halfNSamples) + 1):numSamples)];
    else % Even number of samples -> Special Case
        vXDftInt = interpFactor * [vXDft(1:halfNSamples); vXDft(halfNSamples + 1) / 2; vXDft((numSamples - halfNSamples + 2):numSamples)];
    end
end

So we take care of 2 cases here:

Upsample
We add zero samples to the center part of the DFT to match the number of samples of the output (numSamplesO).
We take care of the case the input number of samples (numSamples) is even. In that case we split the Nyquist sample ($ X \left[ N / 2 \right] $) into 2 where $ N $ is the input number of samples.
Downsample
We remove samples of the center part of the DFT to match the number of samples of the output (numSamplesO).
We take care of the case the output number of samples (numSamplesO) is even. In that case we split the to be Nyquist sample ($ X \left[ M / 2 \right] $) into 2 where $ M $ is the output number of samples.

The question is, why do we do it this way? Why the interpolation factor interpFactor? Where does the splitting factor of $ 0.5 $ come from?
To answer that we need to remember the DFT is basically the Discrete Fourier Series (DFS).
This means the most important assumptions is the data being periodic both in time and frequency domain.
Now, since the DFT is basically the DFS the natural way to interpolate a signal within its period would be using the Fourier Series.
Before going into details let's define 2 sets of integer numbers which will be used to define the values of indices:
$$
\begin{aligned}
\mathcal{K}_{DFS}^{N} & = \left\{- \left\lceil \frac{N - 1}{2} \right\rceil, - \left\lceil \frac{N - 1}{2} \right\rceil + 1, \ldots, -1, 0, 1, \ldots, \left\lceil \frac{N - 1}{2} \right\rceil - 1, \left\lceil \frac{N - 1}{2} \right\rceil \right\} \\
\mathcal{K}_{DFT}^{N} & = \left\{- \left\lceil \frac{N - 1}{2} \right\rceil, - \left\lceil \frac{N - 1}{2} \right\rceil + 1, \ldots, -1, 0, 1, \ldots, \left\lceil \frac{N - 1}{2} \right\rceil - 1, \left\lfloor \frac{N - 1}{2} \right\rfloor \right\} \\
\end{aligned}
$$
This means, for a signal with maximum bandwidths of $ \frac{1}{2 T} $ sampled by Sampling Theorem for $ t \in \left[ 0, N T \right) $ where $ T $ is the sampling period and $ P = N T $ is the function period:
$$
\begin{aligned}
x \left( t \right) {\Big|}_{t = n T} & = \sum_{k = - \left\lceil \frac{N - 1}{2} \right\rceil}^{\left\lceil \frac{N - 1}{2} \right\rceil} {c}_{k} {e}^{ j 2 \pi \frac{k t}{P} } && \text{By Fourier Series} \\
& = \sum_{k = - \left\lceil \frac{N - 1}{2} \right\rceil}^{\left\lceil \frac{N - 1}{2} \right\rceil} {c}_{k} {e}^{ j 2 \pi \frac{k t}{N T} } && \text{By the period of the function / series} \\
& = \sum_{k = - \left\lceil \frac{N - 1}{2} \right\rceil}^{\left\lceil \frac{N - 1}{2} \right\rceil} {c}_{k} {e}^{ j 2 \pi \frac{k n}{N} } && \text{Setting $ t = n T $} \\
& = \frac{1}{N} \sum_{k = - \left\lceil \frac{N - 1}{2} \right\rceil}^{\left\lfloor \frac{N - 1}{2} \right\rfloor} X \left[ k \right] {e}^{ j 2 \pi \frac{k n}{N} } && \text{The  DFT}
\end{aligned}
$$
The formula above works for the even case $ N = 2 l, \; l \in \mathbb{N} $ and for the odd case $ N = 2 l + 1, \; l \in \mathbb{N} $.
The above defines the connection between the DFT coefficients and the Fourier Series Coefficients:
$$ {c}_{k} = \begin{cases}
\frac{ X \left[ k \right ] }{2 N} & \text{ if } k = \frac{N}{2} \\
\frac{ X \left[ k \right ] }{2 N} & \text{ if } k = -\frac{N}{2} \\
\frac{ X \left[ k \right ] }{N} & \text{ if } k \notin \left\{\frac{N}{2}, -\frac{N}{2} \right\}
\end{cases}, \; k \in \mathcal{K}_{DFS}^{N} $$
But there is also nothing stopping us to use other sampling points for any set $ { \left\{ {t}_{m} \right\}}_{m = 0}^{M - 1} $ where $ \forall m, {t}_{m} \in \left[ 0, N T \right) $. Which gives $ x \left( t \right) = \frac{1}{N} \sum_{k = - \left\lceil \frac{N - 1}{2} \right\rceil}^{\left\lfloor \frac{N - 1}{2} \right\rfloor} X \left[ k \right] {e}^{ j 2 \pi \frac{k t}{N T} } $ for $ t \in \left[ 0, N T \right) $. This will work for complex and real signals.
For real signals, $ x \left( t \right) \in \mathbb{R} $ we can also use the Cosine form of the DFT:
$$
\begin{aligned}
x \left( t \right) & = \sum_{k = - \left\lceil \frac{N - 1}{2} \right\rceil}^{\left\lceil \frac{N - 1}{2} \right\rceil} {c}_{k} {e}^{ j 2 \pi \frac{k t}{N T} } && \text{From the above} \\
& = \sum_{k = - \left\lceil \frac{N - 1}{2} \right\rceil}^{\left\lceil \frac{N - 1}{2} \right\rceil} \left| {c}_{k} \right| \cos \left( 2 \pi \frac{k t}{N T} + \angle {c}_{k} \right) && \text{Fourier series in its Cosine form} \\
& = \sum_{k = - \left\lceil \frac{N - 1}{2} \right\rceil}^{\left\lfloor \frac{N - 1}{2} \right\rfloor} \frac{\left| X \left[ k \right] \right|}{N} \cos \left( 2 \pi \frac{k t}{N T} + \angle X \left[ k \right] \right) && \text{Fourier series in its Cosine form} \\
& = \sum_{k = 0}^{\left\lfloor \frac{N - 1}{2} \right\rfloor} {\alpha}_{k} \frac{\left| X \left[ k \right] \right|}{N} \cos \left( 2 \pi \frac{k t}{N T} + \angle X \left[ k \right] \right) && \text{Using the DFT conjugate symmetry of a real signal}
\end{aligned}
$$
Where $ {\alpha}_{k} = \begin{cases} 1 & \text{ if } k \in \left\{ 0, \frac{N}{2} \right\} \\
2 & \text{ else } \end{cases} $.
So now we need to think over what we saw here and how it relates to the algorithm above.
First we need to pay attention that the main trick here is that the native form of the DFT should be when the index goes $ k \in \mathcal{K}_{DFT}^{N} $. Then it is easier to see the connection to the Discrete Fourier Series (DFS) origins of the DFT.
Remark: In practice, the DFT is defined (And computed) with $ k \in \left\{ 0, 1, \ldots, N - 1 \right\} $.
If we chose the set of the output uniform time grid $ { \left\{ {t}_{m} \right\}}_{m = 0}^{M - 1} $ to be in the form $ {t}_{m} = m {T}_{s} $ where the upsampling rate (We'll take care of downsampling later on) $ q = \frac{M}{N} \geq 1 $ then it is clear what need to be done by lookin at the IDFT to recover a grid:
$$ x \left[ m \right] = \frac{1}{M} \sum_{k = 0}^{M - 1} \tilde{X} \left[ k \right] {e}^{j 2 \pi \frac{k m}{M}} = \frac{1}{M} \sum_{k = - \left\lceil \frac{M - 1}{2} \right\rceil}^{\left\lfloor \frac{M - 1}{2} \right\rfloor} \tilde{X} \left[ k \right] {e}^{j 2 \pi \frac{k m}{M}} $$
Now we need to make this match the interpolation formula from above. Since it is a linear transformation multiplying it by $ q $ will take care of the constant. We can also notice that $ \forall m, \frac{m}{M} = \frac{{t}_{m}}{N T} $ hence by setting:
$$ \tilde{X} \left[ k \right] = \begin{cases} X \left[ k \right] & \text{ if } k \in \mathcal{K}_{DFT}^{N} \setminus \left\{ k \mid k = \frac{N}{2} \right\} \\ \frac{X \left[ k \right]}{2} & \text{ if } k = \frac{N}{2} \\ 0 & \text{ if } k \notin \mathcal{K}_{DFT}^{N} \end{cases} $$
From the $ N $ periodicity of the DFT we can write the final interpolation for a uniform grid of time with interpolation factor of $ q $:
$$ x \left[ m \right] = \frac{q}{M} \sum_{k = 0}^{M - 1} \hat{X} \left[ k \right] {e}^{j 2 \pi \frac{k m}{M}} $$
Where $ \hat{X} \left[ k \right] $ is defined as:
$$ \hat{X} \left[ k \right] = \begin{cases} X \left[ k \right] & \text{ if } k \in \left\{ 0, 1, \ldots, N - 1 \right\} \setminus \left\{ \frac{N}{2} \right\} \\ \frac{X \left[ k \right]}{2} & \text{ if } k = \frac{N}{2} \\ 0 & \text{ if } k \in \left\{ N, N + 1, \ldots, M - 1 \right\} \end{cases} $$
Which exactly what we did in the upsample code above.
What about downsample? Well, we can use the same intuition in the DFT domain as the code shows. This is basically because the interpolation using the Fourier Series coefficients is nothing but multiplication in the frequency domain by the Dirichlet Kernel which is the periodic equivalent of the $ \operatorname{sinc} \left( \cdot \right) $ function. This is also the intuition for the $ \frac{1}{2} $ factor, as we multiply with a rectagle with value 1 at the frequency domain which has jump discontinuity. Indeed Fourier Series converges to the mean value of the jump at discontinuties. Since we go from $ 1 $ to $ 0 $, it means the value at the jump is $ 0.5 $.
So the downsmaplign and upsampling code above just applies Dirichlet Kernel to the data according to sampling frequency of the input, in the upsample case and the output in the downsample case.
Another method to downsample would be upsampling to an integer factor of the output number of samples. Then use decimation (Take every ... sample) to get the samples. The 2 will match for the case the data has no energy in the frequency between the low rate and sampled rate. If it does, they won't match.
I will add MATLAB Code...
Remark: This answer also covers the Upsampling. Please consider opening another question on Upsampling or widen this one.


Why does a longer observation time improve DFT resolution, but repeating a signal does not? As was proven here: https://math.stackexchange.com/questions/228614/why-doesnt-repeating-a-signal-give-rise-to-a-finer-resolution-of-dft-fft repeating a certain sequence does not improve DFT frequency resolution.
However, it is known that oberserving a signal for a longer time will in fact improve the ability to resolve two tones in a signal (because the 'beamwidth' of the window in the frequency domain will become narrower: it is 2/NT with NT the total observation time).
I do not understand this discrepancy: If we assume the signal has two tones f1 and f2, and has a period T. Let's say we observe it over a time N*T such that 1/(NT) becomes small enough to distinguish f1 and f2. We have improved the frequency resolution by observing the signal over a longer time. Why is this not equivalent to simply observing the signal for 1 period, and then paste it together N times? This will give the same signal in the time domain, but according to the proof in the link that I posted, it will not improve the frequency resolution.
The only explanation that I can think of myself is: observing for a longer time period does not improve the frequency resolution further once the observation time becomes larger than T. However, I have written some matlab code where I do an FFT on a pure sine over a certain time. I clearly see the 'beam width' in the frequency domain associated with the windowing start to narrow down as I increase the observation time. Given that a pure sine is also periodic and the observation times that I used are bigger than this period, I must thus conclude that my 'explanation' is incorrect.
Could someone try to show me the correct explanation?

Why is this not equivalent to simply observing the signal for 1 period, and then paste it together N times?

It's only equivalent if certain conditions are met. Let's look at a single sine wave with frequency $f_0$ and period $T_0 = 1/f_0$. Sample rate is $f_s$, sampling period $T_s = 1/f_s$ and DFT length is $N$.
The input to the FFT would be
$$x[n] = sin(2\pi \cdot f_0\cdot T_s \cdot n), 0 \leq n < N$$
Now let's create two signals of length $2N$: one with double the length and one with pasting two buffers together.
$$x_1[n] = sin(2\pi \cdot f_0\cdot T_s \cdot n), 0 \leq n < 2N$$
and
$$x_1[n] = sin(2\pi \cdot f_0\cdot T_s \cdot mod(n,N)), 0 \leq n < 2N$$
The only difference is that modulo function in the second one, but it's important. It basically resets the counter to 0 when it reaches an integer multiple of N. In general this will create a discontinuity in the signal. It forces $x_2[N] = x_2[0] = 0$ even if $x[N]$ would be non zero. The first halves of $x_1$  and $x_2$ are the same but the second halves are not. $x_1$ is still a sine wave but $x_2$ isn't anymore.
The ONLY case where pasting is the same as making the signal longer is if the period of the signal is an integer multiple of the sampling period. In this case, there is no discontinuity and both cases are identical. And sure enough, in this case the DFT will resolve two different frequencies perfectly (as long as both meet the integer multiple condition).
A more philosophical interpretation: Making the signal longer provides more information. Pasting buffers just repeats information you already have, so you are not gaining anything new.


Why the Nyquist frequency is 0.5 of Fs, why not 0.55 or 0.65?, brief explanation This my elaboration of the aliasing issue:
a continuous signal can be represented by factors of :
$e^{(i2{\pi}ft)}$ if we sample this signal then I will get:
$e^{(i2{\pi}fk/N)}$ where $k=0,1,2.., N-1$
this point can be represented by the below factor as well:
$e^{(i2{\pi}f_pt)}$
if we put all together:
$e^{(i2{\pi}ft)} = e^{(i2{\pi}f_pt)}$
$e^{(i2{\pi}(f_p-f)t)} = 1$
this is satisfied only if :
$(f_p-f) 2{\pi} /N = 2{\pi}m,$ this lead to
$f_p-f = N*m$
now here it comes the analysis
$N$ is actually the sampling frequency, let say $N$ = 100
so if a say my signal is a single frequency $f$ = 100 then my aliasing frequency (or the other signal frequency) would be:
$f_p = f + m*N = 100 + m*100 = 100(1+m),$
$f_p(m=-1) = 0,$ issue,
$f_p(m=0) = 100,$ no issue,
$f_p(m=1)= 200,$..issue,
now if my $f = 50$
$f_p = f + m*N = 50 + m*100 = 100(0.5+m),$
$f_p(m=-1) = -50,$ issue
$f_p(m=0) = 50,$ no issue
$f_p(m=1)= 150,$..issue
now if my $f = 25$
$f_p = f + m*N = 25 + m*100 = 100(0.25+m),$
$f_p(m=-1) = -75,$ no issue
$f_p(m=0) = 25,$ no issue
$f_p(m=1)= 125,$.no .issue
now if my $f = 75$
$ff = f + m*N = 75 + m*100 = 100(0.75+m),$
$ff(m=-1) = -25,$ issue
$ff(m=0) = 75,$ no issue
$ff(m=1)= 175,$.no .issue
so then after tabulating all possible frequencies I can say that the best frequencies where $f_p$ doesn't exist is when $f < N / 2$
so What do you think about this reasoning. is there other way to explain the 1/2 term that multiplies $N$ to give us the right signal frequency ( that can be correctly sampled with $N$  samples)?
Appreciate your opinion

From the definition of sampling, and the definition of periodicity follows that any signal is identical to signal at a frequency higher by an integer amount of sampling frequencies.
Since for real signals, negative frequencies aren't distinguishable from positive ones, it directly follows that for every $0.5f_\text{sample}\le f \le f_\text{sample}$, there is a negative frequency $-0.5f_\text{sample}\ge \tilde f=1- f \ge f_\text{sample}$ that looks the same as a positive frequency.
Hence, only frequencies lower than $|f| < 0.5f_\text{sample}$ can be uniquely represented.


Combination of accelation signals from two devices I have $x_A$, $y_A$, and $z_A$ acceleration signals collected from a smartphone and I call it: set $A$ and I also  have $x_B$, $y_B$, and $z_B$ same acceleration signals from smartwatch: let's call this set $B$. I need to combine them together to obtain set $C$, but since I don't have good knowledge in signal processing I am not sure if it is this logically correct to combine them and take the average of them, or any other way should I follow to have one set of signals instead of having two separate sets?

How to select the set if it is possible? The signals are raw data.

May be you can find variance of acceleration along x,y, and z direction in SET A and in set B separately and choose the the device which has least variance.

As this Wikipedia page suggests, you can combine the two by weighting by their respective variances.

So, to get $x_C$ simply:
$$
x_C = \frac{1}{\sigma_1^{-2} + \sigma_2^{-2}} \left(  \frac{x_A}{\sigma_1^2} + \frac{x_B}{\sigma_2^2}\right)
$$

Putting some numbers on this, suppose $\sigma^2_1 = 1$ and $\sigma^2_2 = 0.5$, then
$$
\sigma_3^2 = \frac{1}{\sigma_1^{-2} + \sigma_2^{-2}} = \frac{1}{1 + 2} = \frac{1}{3}
$$
which is a lower variance than either of the individual measurements.

Use the same process to find $y_C$ and $z_C$.


Is it always necessary to use two cameras to find the distance between camera and object? I am working on finding the distance between a camera and an object using image processing.

In many places, it's written that we have to use a stereo pair, but my project allows only one camera.  Is it feasible to find distance to an object using just one camera? If so, what are the possible ways to do this?

There are various ways to measure distance and 3D using a single camera.
These include:

Moving the camera and taking another image - this simulates a 2 camera system.
If you are imaging an object of a known size, then you can know the distance with some simple trigonometry.
You can project a known pattern from a known projector location to reconstruct 3D. This is structured light and e.g. how Kinect works. Obviously the image(s) with the pattern will appear different from the normal image, so multiple exposures might be neccesary if you also want the appearance image.
You can use Depth-from-Defocus by capturing multiple images at multiple focii.
Sometimes an auto-focus lens can provide rough estimates on the depth of the object in focus. If the object is big enough in the frame this may be sufficient for some applications.

Each method has different trade-offs. You must consider which one suits your app.