What Information does the phase of a cross power spectra give me? I obtain the cross power spectra by the following steps:

compute the FFT of signal A and B
Multiply A with the conjugate of B store it in C (cross power spectra)
Now looking at the Phase of this cross power spectra C, what all information can I deduce about the two signals A and B?

If the two signal show a significant correlation at some frequency $f$ then the phase measured in the cross spectrum gives you the phase lag of the signals at that frequency.

So, take $$ A = \sin(2\pi f t ) $$ and $$ B = \sin(2\pi f t - \phi) $$ then the cross spectral phase at frequency $f$ will give you $\phi$.


Amplify Sound From Known Direction I have microphone array which records sound. There is one main sound which direction I know and some noises in the background. Now I want to selectively amplify sound from known direction. This looks like reverse problem of antenna signal beamforming. I tried to find some information on internet on how to do it but only information I found told that I need to shift sound signals from microphones. But Wikipedia says it works only for narrow band signals while my signal will be more complicated than single sinus.

I tried to search internet for beamforming but all I found is about EM waves and antennas. Are there some known methods to solve this problems?

I just looked at the Wikipedia article on beamforming and i can understand your confusion. It ‘s trying to explain too much at once. There are some nomenclature issues as well, MUSIC is a direction finding technique. It doesn’t form beams in the literal sense.

Unfortunately a lot of articles as well as books are written by people who understand the material for other people who understand the material.

You should also be aware that beamforming isn’t always or even practical solution for many audio situations. Shotgun microphones or parabolic reflectors are often good and potentially cheaper solutions.  People who study bird songs tend to use these.

In addition to your known direction, you need known locations of your microphones.

The simplest way to beamform is to place all your microphones on a plane that is perpendicular to the direction of the source you want to amplify and then sum the outputs of all the microphones. This is based on the assumption that the source is far enough away so that you can neglect the curvature of the wavefront. This assumption depends on the assumption that the source can be considered a point emitter.

Another way to put this is that the acoustic delay time from the source to each microphone is the same.  The Wikipedia article used the term constructive interference to say the same thing.

If locating the microphones on the coherent wavefront is not possible, one introduces delays in each microphone so that each summed microphone has the same acoustic distance travel time from the source when coherently summed.  This introduction of delays can be very general and permit the array to coherently sum the sounds from many chosen directions. Since you say you have recordings, this is your task.

Each microphone should have a clear line-of-sight to your source. Microphones typically have a directional response as well so it helps to point each microphone at the source as well. One cannot compensate for the individual directivity of each microphone by introducing delays. Actually a lot of the material written about arrays makes the unrealistic assumption that each microphone has an omnidirectional response.

There are numerous effects that complicate conventional delay-and-sum beamforming as i have described it.  The source sounds can arrive at the array from more than one direction which is called multi-path. The signal can cause destructive interference with itself in this circumstance.  Inside of rooms there is a lot of potential multipath. The source can’t be considered a simple point source.  These are only a few complications.

You can have another loud source that impinges on the array, so even if not coherently summed, still interfere with the source of interest.  The way to treat this is to steer a second beam towards that source and judiciously subtract its sound from the sound  of interest.  There is an entire literature on optimal beamforming that treats this situation.

Summarizing, beamforming may or may not be worth it. It takes specialized hardware like synchronized multichannel analog-to-digital converters.

It isn’t cheap compared with a shotgun microphone. Beamforming is often done because nothing else is possible.


Discrete-time sampling of filtered white noise I am trying to understand how I can relate a discrete-time random process to a continuous-time random process sampled at discrete times.

Suppose I have a noise source $N_\tau(t)$ which is derived from unit-amplitude additive white Gaussian noise N(t) that is fed through a single-pole low pass filter $H(s) = \frac{1}{\tau s +1}$.

I understand that the power spectral density of $N_\tau(t)$ is $S_{N_\tau}(\omega) = |H(\omega)|^2 = \frac{1}{(\omega\tau)^2  +1}$. So far so good.

Now I want to sample $N_\tau(t)$ at regular intervals $T$ to obtain a signal $n_\tau[k] = N_\tau(kT)$ -- how would I figure out the standard deviation and spectral density of the discrete samples $n_\tau[k]$?

Hmm. For standard deviation, I see https://dsp.stackexchange.com/a/8632/829 which states for uniform power spectral density $N_0/2$, the standard deviation is

$$\sigma^2  = \int_{-\infty}^\infty \frac{N_0}{2}|H(f)|^2\,\mathrm df$$

which in my case would yield

$$\begin{align}
\sigma^2  &= \int_{-\infty}^\infty \frac{N_0}{2}|H(f)|^2\,df \\
          &= \int_{-\infty}^\infty \frac{N_0}{2}|H(\omega/2\pi)|^2\,d(\omega/2\pi) \\
          &= \frac{1}{2\pi}\int_{-\infty}^\infty \frac{N_0}{2}|H(\omega)|^2\,d\omega \\
          &= \frac{1}{2\pi}\int_{-\infty}^\infty \frac{N_0}{2}\frac{1}{(\omega\tau)^2+1}\,d\omega \\
          &= \frac{1}{2\pi\tau}\int_{-\infty}^\infty \frac{N_0}{2}\frac{1}{(\omega\tau)^2+1}\,d(\omega\tau) \\
          &= \frac{1}{2\pi\tau}\int_{-\infty}^\infty \frac{N_0}{2}\frac{1}{u^2+1}\, du \\
          &= \frac{N_0}{4\pi\tau}\left[\tan^{-1} u\right]_{-\infty}^{\infty} \\
          &= \frac{N_0}{4\pi\tau}\left[\pi/2 - (-\pi/2)\right] \\
          &= \frac{N_0}{4\tau},
\end{align}
$$

and with unit WGN $N_0/2 = 1$ so $\sigma = \sqrt{1/2\tau}$.

Not sure how to compute the PSD of a sampled random process, however.


How many psychoacoustical critical bands are there? I understand that the human ear works as a band-pass filter and I have doubt regarding the number of such bandpass filters in the filter bank.

I came across these two possible explanations of how it happens:

There are no fixed bands and the critical bandwidth of a sound depends on which audio frequency the ear is exposed to. In such a case, the audio frequency could be treated as the center frequency of a bandpass filter and another frequency occurring in the same band gets masked out.
There are 24 fixed bandpass filters with their corresponding fixed bandwidths independent of the audio frequency the ear is exposed to. In such a case too, frequencies occurring in the same band sound like a single frequency and the difference goes undetected.

Which one of them is correct? Do the critical bands exist only between the 24 fixed values or are they continuous values?

Neither. The basilar membrane in your ear performs a frequency-to-place transformation which is then picked up by an array of thousands of hair cells (cilia). Therefore there are thousands of heavily overlapped bands. The processing that follows is very complex and not completely understood. The shape of the effective band pass filter at each location on the membrane can be measured in animals, and the masking curves can be measured by listening tests. Masking is measured by increasing the amplitude of a test tone relative to the masking tone until it is audible, so the situation is a bit different that what you assumed in your question (where you assumed that if the masker and the test tone were in the same band, you only hear one tone). The actual result is that as the test tone approaches the masker frequency, it must have a higher amplitude before it can be detected.

Bob


Does the 'FFTW' have a cost of k*N*log_2(N)? @Michael C. Grant wrote that "the cost of FFTW isn't an easy formula based on log_2 anymore."

But the Wikipedia article says that the FFTW "can compute transforms of real and complex-valued arrays of arbitrary size and dimension in O(n log n) time."

So can it compute 1 forward FFT with the size 1920*1080? in kNlog(N)/log(2) =

5*1920*1080*log(1920*1080)/log(2)=
217559066 real operations?

Or should I perphaps take a different value for k?

You are trying to describe a complex phenomenon with a simple formula. That gets you only a rough approximation or an estimate but depending how much accuracy you need, you will have to dive deeper into the specifics. A few points

Any actual operation count will depend heavily on the hardware: Instruction set, pipeline, memory architecture, cache strategies, hardware accelerators, etc. In many cases the choice of the "best" algorithm is actually hardware dependent. Your best shot here is benchmarking and profiling on the target hardware.
The good old N/2*log2(N) formula basically describes the "number of butterfly operation" in a standard Cooley Tuckey FFT. So a 32-point FFT has 16*5 = 80 butterfly operations. There are 5 stages of 16 butterflies each. This is typically (but not always) the bulk of the work so it's a decent estimate for the overall complexity.
One butterfly has one complex multiply and two complex adds (or on
add and one subtract if you are picky about). The number of operations required will depend widely on the Hardware. An Analog Devices Sharc SIMD processor can do this with about 3 operations, an general purpose processor might take up to 20 operations.
There is a lot of other stuff that isn't part of the butterflies: bit-reversing, looping & branching, coefficient calculation storage and retrieval. These can factor heavily into your MIPS count as well and are even more hardware dependent.
The first two stages (for decimation in time) use only simple coefficients like +1, -1, +j, -j and don't need a full complex multiply. Most algorithms take advantage of this and so the first 2*log2(N) butterflies are a little cheaper. By how much depends again on the hardware.
Non power of two FFTs get more complicated. Let's look at a 96 point FFT. First you need 3 FFTs of 32 points each, so that is 3*5*16 = 240 butterflies. Then you need to assemble the result which takes 2 complex multiplies and 2 complex adds per sample, so it's another 96*2=192 extra complex multiplies and 192 complex adds. Total would be 432 mpy and 672 adds. A 128 FFT would have 448 mpy and 896 adds, so it may be better just to zero pad and use the larger power-of-two FFT as it greatly simplifies coefficient calculation and storage.
The more prime factors you have in the FFT length the more complicated and less efficient the algorithm gets.
In your specific case, I'd zero pad the 1920 to 2048 (unless you need true circular convolution, which I don't believe you do here). It's highly unlikely that any prime factorization (128*3*5) would do better. The 1080 is awkward, since it's just a tad over 1024. I'd probably zero pad to 1280 and do a 256*5 or zero pad to 1152 and do 128*9. You can also try the 8*27*5 and see if that's better.


Real and Imaginary Frequency Responses of a single complex pole I'm filtering a real signal with a single complex pole with a complex coefficient (a_re [real part] and a_im [imaginary part]), I also have a gain coefficient but I'm gonna leave it out for the sake of simplicity. The filter has two outputs, one for the real and another for the imaginary part. Here's the filter equation in time for both outputs

Real part output:
$$y_1[n] = x_1[0] + a_{R}  y_1[-1] - a_{I} y_2[-1]$$
Imaginary part output:
$$y_2[n] = a_{I}  y_1[-1] + a_{R} y_2[-1]$$

Where $n$ is the sample number, $x_1[0]$ is the real input, $a_{R}$ & $a_{I}$ are, respectively, the real and imaginary part of the coefficient, $y_1[-1]$ is the previous output of "real part output" and $y_2[-1]$ is the previous output of "imaginary part output".

So, I need to plot the frequency response based on the Z transform. I've already succeeded in dealing with the real version of all this (a real pole), and I'm struggling with the complex version. Anyway, the equation for the real pole is:

$$y[n] = x[0] + a y[-1]$$

Where $a$ is the coefficient of the real pole.

Now, the Z-transform of the real pole is:

$$\tag{1}H(z) = \dfrac{1}{ 1 - a z^{-1}}$$

Where, again, $a$ is the coefficient of the real pole.

In my code I am solving it and dealing  with the z-transform in this way, the input of the function is the $a$ coefficient and the $\omega$ variable (which is the angular frequency in radians per sample). I'm only interested in getting the magnitude output, so I don't care about the phase response of the filter. Anyway, here is how I get it:

$a \rightarrow$ coefficient

$w \rightarrow$ angular frequency

$f_{R} = \cos(\omega)$; $\;\;\;\;f_{I} = \sin(\omega)$

$\mathbf{R} = 1 - a  f_{R}$; $\;\;\mathbf{I} = -a f_{I}$

$\mathbf{Mag} = \sqrt{\mathbf{R}^{2} + \mathbf{I}^{2}}$

$H = \frac{1}{\mathbf{Mag}}$

The whole thing would be:

$$\tag{2}|H(z)| = \dfrac{1}{\sqrt{(1 - a \cos \omega)^2 + (-a \sin \omega)^2)}}$$

Now my trouble is adapting it to the complex version of this z transform, which is supposed to be the same formula as the real pole's z transform in Eq. (1), but I don't know how to adapt the formula/code for the complex pole. Bearing in mind that I want to have the frequency response of the two outputs (the real and imaginary part of the complex filter).

Hopefully, I'd like some help on reaching the frequency response from the coefficient as a function of angular frequency as I've written above for the real pole version of this in Eq.(2)!

Any thoughts, considerations, hints, help is highly welcome.

Thanks a lot

If I understood you correctly you want a formula for the magnitude of the frequency response of a filter with a single complex pole. As you've already noticed, the system function is

$$H(z)=\frac{1}{1-az^{-1}}\tag{1}$$

which is valid regardless of whether $a$ is real-valued or complex-valued. Let's assume that $a=a_r+ja_i$ is a complex-valued number. If we also assume that $|a|<1$ (i.e. the system is causal and stable), then we get the systems frequency response by replacing $z$ by $e^{j\omega}$ in (1):

$$H(e^{j\omega})=\frac{1}{1-ae^{-j\omega}}\tag{2}$$

Now you just have to compute the magnitude of (2). Using

$$|1-u|^2=1-2\text{Re}\{u\}+|u|^2$$

for any complex $u$, we get with $u=ae^{-j\omega}$ from (2)

$$|H(e^{j\omega})|=\frac{1}{\sqrt{1-2\text{Re}\{ae^{-j\omega}\}+|a|^2}}=
\frac{1}{\sqrt{1-2(a_r\cos\omega+a_i\sin\omega)+|a|^2}}\tag{3}$$

EDIT:

This is how you can derive the transfer functions relating the real-valued input to the real part and the imaginary part of the output, respectively. Let $y(n)=y_r(n)+jy_i(n)$. Then you have the following two difference equations for the two real-valued output signals:

$$y_r(n)=x(n)+a_ry_r(n-1)-a_iy_i(n-1)\\
y_i(n)=a_ry_i(n-1)+a_iy_r(n-1)$$

Transforming these equations to the $\mathcal{Z}$-domain gives

$$Y_r(z)=X(z)+a_rz^{-1}Y_r(z)-a_iz^{-1}Y_i(z)\\
Y_i(z)=a_rz^{-1}Y_i(z)+a_iz^{-1}Y_r(z)$$

Now you have two equations for the two unknowns $Y_r(z)$ and $Y_i(z)$. After some algebra you arrive at

$$Y_r(z)=H_r(z)X(z)=\frac{1-a_rz^{-1}}{1-2a_rz^{-1}+(a_r^2+a_i^2)z^{-2}}X(z)\\
Y_i(z)=H_i(z)X(z)=\frac{a_iz^{-1}}{1-2a_rz^{-1}+(a_r^2+a_i^2)z^{-2}}X(z)\tag{4}$$

Note that even though

$$H(e^{j\omega})=H_r(e^{j\omega})+jH_i(e^{j\omega})$$

$H_r(e^{j\omega})$ and $H_i(e^{j\omega})$ are not the real and imaginary parts of the total frequency response $H(e^{j\omega})$ of the complex-valued system because they are both themselves complex-valued. You can obtain the magnitudes of the frequency responses by computing the magnitudes of $H_r(e^{j\omega})$ and $H_i(e^{j\omega})$ as given in (4), i.e. by replacing $z$ by $e^{j\omega}$ and computing the magnitudes of the two complex functions.

Another option is to express the frequency responses $H_r(e^{j\omega})$ and $H_i(e^{j\omega})$ using the total frequency response $H(e^{j\omega})$ (2):

$$H_r(e^{j\omega})=\frac12[H(e^{j\omega})+H^*(e^{-j\omega})]\\
H_i(e^{j\omega})=\frac{1}{2j}[H(e^{j\omega})-H^*(e^{-j\omega})]\tag{5}$$

From (5), the magnitudes can be written as

$$|H_r(e^{j\omega})|=\frac{1}{2}\sqrt{|H(e^{j\omega})|^2+2\text{Re}\{H(e^{j\omega})H(e^{-j\omega})\}+|H(e^{-j\omega})|^2\}}\\
|H_i(e^{j\omega})|=\frac{1}{2}\sqrt{|H(e^{j\omega})|^2-2\text{Re}\{H(e^{j\omega})H(e^{-j\omega})\}+|H(e^{-j\omega})|^2\}}$$

where $|H(e^{j\omega})|$ is given by (3).

EDIT 2:

Here's some Matlab/Octave code using equation (4) to plot the different transfer functions. $|H_r(e^{j\omega})|$ and $|H_i(e^{j\omega})|$ are of course symmetrical because they correspond to real-valued systems, whereas $|H(e^{j\omega})|$ is not due to the complex pole.

w = linspace(0,2*pi,300);          % frequency vector
ar = .4; ai = -.3; a = ar+j*ai;    % complex pole
u = exp(-j*w);                     % u = z^(-1)
Hr = (1-ar*u)./(1-2*ar*u+(ar^2+ai^2)*u.^2);
Hi = (ai*u)./(1-2*ar*u+(ar^2+ai^2)*u.^2);
H = 1./(1-a*u);
plot(w/pi,abs(Hr),w/pi,abs(Hi),w/pi,abs(H))
legend('Hr', 'Hi', 'H');

EDIT 3:

If you need to split the numerators and denominators of (4) into real and imaginary parts, here's what you get:

$$H_r(e^{j\omega})=\frac{1-a_r\cos\omega+ja_r\sin\omega}{\left[1-2a_r\cos\omega+(a_r^2+a_i^2)\cos(2\omega)\right]+j\left[2a_r\sin\omega-(a_r^2+a_i^2)\sin(2\omega)\right]}\\
H_i(e^{j\omega})=\frac{a_i\cos\omega-ja_i\sin\omega}{\left[1-2a_r\cos\omega+(a_r^2+a_i^2)\cos(2\omega)\right]+j\left[2a_r\sin\omega-(a_r^2+a_i^2)\sin(2\omega)\right]}$$


Poles of an analog Gaussian filter I am trying to calculate the poles of an analog Gaussian filter. Its characteristic function, $e^{-\log(2)x^2}$, can be expanded into MacLaurin series:

$$2^{-x^2} = e^{-\log(2)x^2} = \lim_{N \to \infty} \sum^N_{k=0}{\frac{(\log(2))^k}{k!}x^{2k}}$$.

This would go to the denominator of the transfer function, so the roots of the sum are the poles of the filter. The more terms, the better the approximation. So far, so good, classic theory.

But if I calculate the poles for, say $N$=13, I get some strange display (blue dots):

Just by looking at the graph, I thought to rotate them 90 degrees and keep only the left side (red points), while comparing them to the Bessel poles of the same order (green points). The circle is for reference, only, while the Bessel poles have been scaled down by $\sqrt{N}$, for comparison. Interesting enough, if I am to calculate the poles for the normalized transfer function (not the squared one), the series terms would have to be divided by an extra $2^k$, and the resulting transfer function would have to be compared with $\sqrt{|H(s)^2|}$, rather than $|H(s)|$.

I am aware that Bessel isn't Gaussian, or vice-versa, but they are two filters who deal with time-domain, rather than frequency, and their responses are quite similar, hence the comparison.

But now I am puzzled: the blue poles are symmetric on the $X$ axis, and if I were to keep the left hand side, I would end up having two extra, purely imaginary poles for odd orders (as is the case here). Rotating them by 90 degrees seems not only to solve this problem, but also to give the correct transfer function magnitude:

The blue trace is with the original poles, left-hand side poles. One note: keeping the extra, purely imaginary poles would result in a resonance peak (not surprisingly), so I removed those and added an extra, single, real pole on $X$ axis, of magnitude $-2|\max\left({s_k}\right)|$, resulting in the current plot. The red trace is with the rotated poles (also left-hand side), while the green one is the Mathematica expression $\sqrt{e^{-\log(2)*x^2}}$ with a 0.1 offset, for better comparison. $N$=13, as above.

All these seem confusing to me: when calculating the poles from the mathematical expression (the series expansion), do we keep them as they are, or do we rotate them 90 degrees (simple multiplication with $j$), and why?

First of all, computing the poles of an (ideal) Gaussian filter is an impossible task, because its transfer function is not a rational function, and there are simply no poles. This is in contrast to the classic filter approximations, like Butterworth, Chebyshev, and Cauer filters, all of which have rational transfer functions, and, consequently, all of them have poles.

The frequency response of an ideal Gaussian filter is

$$H(\omega)=e^{-a\omega^2}\tag{1}$$

The corresponding transfer function is a function of $s=j\omega$:

$$G(s)=e^{-a(s/j)^2}=e^{as^2}\tag{2}$$

which is of course not rational, and which doesn't have any finite poles.

What you did is compute a rational approximation by developing $1/G(s)$ into a power series, and using only a finite number of terms. However, since you used $\omega$ as a complex variable instead of $s$, the poles are turned by $90$ degrees. This was already mentioned in the comments. The correct way of doing this is as follows:

$$1/G(s)=e^{-as^2}=\sum_{k=0}^{\infty}\frac{(-a)^k}{k!}s^{2k}\tag{3}$$

Using only a finite number of terms gives the following rational approximation of $G(s)$:

$$G(s)\approx\frac{1}{\sum_{k=0}^{n}\frac{(-a)^k}{k!}s^{2k}}\tag{2}$$

the poles of which are of course the roots of the denominator polynomial.


Barcode localization in difficult environments I'm creating a program for localizing types of barcodes.
(I know plenty exist, that's not the question :) )
I'm focusing for the time being on 1D barcodes.

The difficulty I have is that the barcodes I need to be able to locate could be situated anywhere within an image with other information (text/images/...), inequal lighting conditions (-> gradients), and barcodes could be scaled and orientated in any way.

[Edit] Unlike most barcode-related cases however, my application does not need to be run on a smartphone, so can allow myself some processing bandwidth

I have been researching as much as possible in books and publications, but think it's time to ask the experts...

What is the most appropriate approach for reliably localising 1D barcodes in difficult environments?

Currently, I'm applying otsu thresholding, basic edge detection and hough transforms.
Problems i currently face:
- otsu thresholding is global and doesn't handle local gradients on large images well. are there localized algorithms that exist?
- hough locates the lines of the barcodes, but also other items (in particular text), and it's difficult to filter only the barcodes (thresholding the min number of pixels in the line throws out a lot, but also throws out the smaller barcodes)

Could anyone point me in the right direction?

I'm wondering about Fourier and Gabor filters, but they seem more complex, so want opinions before i start.

PS: I am doing everything in Octave and Java

Other than commercial barcode reading algorithms (many of which fail to read challenging codes), I would like to direct you to this paper which is one of the best academic works in that field in my opinion:

Gallo, Manduchi: Reading Challenging Barcodes with Cameras.

Here is a more recent version:

Gallo, Manduchi: Reading 1D Barcodes with Mobile Phones Using Deformable Templates

Judging by the performance they present it can decode really deformed and low quality barcodes.

Following the localization stage, you will require a robust and fast decoding algorithm. Here is one very strong algorithm from Simsekli and Birdal:

Simsekli, Birdal: A Unified Probabilistic Framework for Robust Decoding of Linear Barcodes.


Why DTFT coefficients are periodic and why continuous Fourier transform coefficients are not periodic? As I understand, when the input signal is discrete in time and we want to find the coefficients of Fourier transform, DTFT is used and the coefficients in frequency domain are periodic, but I can't understand why? If DTFT coefficients are periodic, why the coefficients of a continuous signal in frequency domain are not periodic?

The explanation by @Maximilian Matthé is a standard and formal approach for this question. But I think it is not intuitive and easy to understand the reason. In the following, I will try to explain from another aspect.

First of all, periodicity means infinity. For a signal in the temporal/spatial domain, the periodicity is somewhat easy for realization, such as extending it periodically and infinitely. But for the frequency domain, it is quite different. Different position means different component with different changing speed. Therefore, if periodicity is expected in the frequency domain, infinitely rapid changing is expected in the temporal/spatial domain. In the real world of our life, most signals changes smoothly, thus periodicity is impossible in the frequency domain. While for discrete-time signals, the model is based on the $\delta(t)$ function, for which the changing is infinitely sharp and rapid, thus periodicity becomes possible in this case.


How to find orthonormal basis function in the following digital communication problem? I'm trying to express the following set of signals on orthonormal basis functions. These signals are to be transmitted using a 4-ary modulation scheme.How do i write ψ1(t)and ψ2(t) for the same. Also how to find orthonormal function set for given set of signals.

Let's call your $4$ signals $f_0$, $f_1$, $f_2$ & $f_3$.
Just by eyeballing you could have a basis of two function where

b0 = f0, for 0 < t < Ts/2, 0 otherwise and
b1 = f0, for Ts/2 < t < Ts, 0 otherwise and

So you basically have a one first half ramp up in the first half and one second half ramo down. The signals can be expressed simply as

f0 =  b0 + b1
f1 = -b0 - b1
f2 =  b0 - b1
f3 = -b0 + b1

So it's clearly a base. It's also orthogonal since $b_0  b_1$ is zero for all $t$. It's probably not normal but could be easily made so by proper scaling.


Question about decimation/downsampling I am using an on-chip SAR adc where for some reason I can only get accurate signal reconstruction and filtering when I sample at almost 100 times my desired sampling frequency. Specifically, the minimum working sampling rate right now is 25KHz but the desired sampling rate is 500Hz for a signal with a BW of 0.5-100Hz, therefore I am still slightly oversampling this signal but the 25KHz is clearly overkill.

I want to downsample this signal to that desired sampling frequency to aid in the design and implementation of digital filters in that passband of .5-100Hz. I know a key thing to look out for is aliasing and because of that you should use a digital anti-aliasing filter but is it pretty much the same concept as an analog anti-aliasing filter where I simply lowpass filter the signal bandwidth that I want to reconstruct after sampling? Also are there any other pitfalls/design tips that I should know to be able to implement downsampling in a DSP? Thanks in advance.

The key point is that you remove any signal outside of the final bandwidth before downsampling to avoid aliasing. If you want to do this analog or digitally it doesn't really make a difference technically. From an engineering perspective, digital filtering is typically more economical. As for pitfalls/design tips - if you do use digital filtering keep in mind that you can save CPU cycles by not computing filtered output samples which will ultimately be thrown away by downsampling.

So to conclude, design a 100Hz LPF FIR at the 25kHz sample rate, run the filter on every 50th sample and then throw all but those samples away.


Relation between signal processing and control systems engineering? Control systems engineering and Digital signal processing are both important courses/subjects of electrical engineering But how these two subjects/courses are related to each other??

Also please kindly let me know,what are some recommended resources (books, tutorials, lectures, etc.) on control systems engineering, and how to begin working with it on a technical level?

As we have answer in below link but that answer is about dsp resources, I am looking for resources about control systems engineering  example answer of refrence request

There is a lot of overlap but some differences in emphasis. Control Engineering is also older than DSP.  If you have a traditional EE education, you don’t really make much of a distinction.

State variables are the more typical perspective in Controls. The first edition of Oppenheim and Schafer 1975, had a chapter on state variables, but they dropped it over the years.  You need to understand state variables to do Kalman Filtering which is an area of overlap. Linear Estimation and Linear Controls are duals of each-other.

I would also say that hybrid continuous/discrete time systems are more common in Controls but there are many examples for DSP as well.

DSP is almost always done on uniform sampling. State Variables can work with nonuniform sampling as well.

I’ve never heard of anti causal Control System but forward backward filtering in time is common in DSP. Controls are inherently causal. The one sided Laplace transform is more common in controls.

Stability in feed back loops is important in both areas.  An advanced control systems class will cover topics like Lyaponov stability.  You typically don’t see that covered in DSP but there are DSP papers that use that technique.

Control Theory shows up in mechanical engineering. DSP shows up in finance. There is a-lot of both in robotics which also uses computer vision.

In RADAR, waveforms and filtering are more DSP at the front end, but the tracking systems at the back end are more Controls like.

If I had to use a single word to describe each.

Controls: feed back

Signal Processing: sensing

or maybe using a phrase

Controls: in-the-present

DSP: in-the-groove


Moving Average Filter's Magnitude Response I am working on filtering of data using the Moving average filter in MATLAB.
My data is recorded temperature values on the Y-axis and Time(in seconds) on the X-axis
(I'm getting the magnitude response curve of the Moving average filter as attached)

The shape of the curve seems fair for the moving average filter, according to the theory but,
I want to ask, why am I getting all the magnitude value for MA filter positive? (is this correct?)
Implementing other filters like Butterworth etc. on the same data gives the magnitude values to be negative (image attached)

Here is my MATLAB Code for MA (if it helps)
%T.Temp is table column containing Temperature values in double format.

MovAvg = dsp.MovingAverage(5)

movAvgOut = MovAvg(T.Temp)

fvtool(movAvgOut)

%Gives output the first image I attached

Any help is appreciated.
Thank you.
Ps - DSP is not my field. This work is part of a bigger project i'm involved in, and need some help. Please bear with me if I sound irrational.

A true moving average would have unity gain coefficients but also be divided by the total number of samples (as per the definition of average). This is trivial since it is just a scaling factor, but if the OP divides the output of the moving average by the total number of samples in the filter then the magnitude response will be normalized to 0 on a dB scale, with all other magnitude values given in negative dB.


How to calculate BW ofdm? HI I have these specification for an OFDM program I was wondering, what is the bandwidth and the sampling frequency? I'm not a student of Digital signal processing, hence I am so confused.

ymax=100;
nIFFT = 2048;
nCyclic = nIFFT/4;
pilot_distance = 16;
pilot_amplitude = 2;
offset = nIFFT*3;
% total number of complex samples in the timedomain
nTotal = (nIFFT+nCyclic)*ymax+offset;

I don't think it's possible to deduce the sampling rate from this information. The sampling rate could be anything, since FFT is sampling rate agnostic. Sampling rate is pretty much independent of all other parameters which scale directly with it.

Concerning the bandwidth, it's not clear either. You're specifying pilot distance, but we don't know what the units are and we never see how those variables get used. Is there more to the code you're studying?

Looked at your full code, and the sampling rate is 8 kHz, since that's what frequency is used in the wavwrite command. Will have to get back to bandwidth.


Detecting water level in an image I have designed a surface with a shape like a rounded rectangle and I want to detect its exact boundaries in a photo I take from it. A sample photo is like this:

I have tried thresholding the image but even after a Gaussian smoothing, the threshold won't give me an exact boundary.

Does anyone know of a technique that would give me the precise boundaries, regardless of the issues in lighting, etc?

UPDATE:
I'm sorry I didn't state the problem well, and actually I wasn't aware of the real problem (though I think it hasn't changed too much):

Here is the new photo. The photo is taken looking down inside a plastic bucket with a small amount of water inside, and I need to segment the region inside the water from the walls of the bucket, so I can detect the objects inside the water later. There is some lighting illuminating the bottom of the bucket from underneath (actually a sheet of plexiglass is placed instead of the original bottom).

The problem is that the reflections of light on the wall make it the same color as where the water is, in some areas like the bottom-center in the sample image which makes it hard to find the water level.

Also, I would appreciate any hints on how should I differentiate between the boundary related to the plexiglass, and the boundary related to the water level.

Thanks

Behold the power of shortest path.

The solution roughly follows what I did here: Detecting and isolating part of an image.

In this case we first treat the image by limiting every pixel value above the water-level grey to that same intensity.

Then performing shortest path on the inverse edge image created with a difference of Gaussian filter we get:

Resulting in the following segmentation:

My MATLAB functions that are needed for this solution:

http://imageprocessing.com.au/research/code/polarTransform.m http://imageprocessing.com.au/research/code/circularShortestPath.m http://imageprocessing.com.au/research/code/linearShortestPath.m

Code:

clear all
close all

%% Load image
I = imread('CabQx.jpg');
I = double(rgb2gray(I));
figure
imagesc(I);
colormap gray(256);
title('Image');

%% Adjust intensity
Ia = I;
Ia(Ia > 49) = 49;
figure
imagesc(Ia);
colormap gray(256);
title('Adjusted Image');

%% Bandpass Kernel
H1 = fspecial('gaussian',[100,100],8); % *** Changed filter size
H2 = fspecial('gaussian',[100,100],2);  % ***
H = H2-H1;
H = mat2gray(H);
figure;
imagesc(H);
colorbar;
colormap gray(256);
title('Bandpass kernel');

%% Edge Kernel
Hedge = imag(hilbert(H));
figure;
imagesc(Hedge);
colormap gray(256);
title('Edge Kernel');

%% Edge filtered image
G = sqrt(imfilter(Ia,Hedge,'replicate').^2 + imfilter(Ia,Hedge.','replicate').^2);
G = max(G(:)) - G;
imagesc(G);
colormap gray(256);
title('Inverse edge filtered image');

%% Shortest path
% Mid point
[r,c] = size(G);
r = floor(r/2);
c = floor(c/2);
% Min radius and max radius
min_radius = 100;
max_radius = 1000;
% Shortest path
[ path, energy, lastIndex, pathImage ] = ...
    circularShortestPath(G, 100, [c,r], [min_radius,max_radius], [400,720*2]);
figure;
imagesc(I);
hold on;
plot(path(:,1),path(:,2),'r','LineWidth',2);
hold off;
colormap gray(256);
title('Shortest path');

%% Segment
figure
imagesc(I);
hold on;
fill(path(:,1),path(:,2),'r');
hold off;
colormap gray(256);
title('Mask')

%%
J = zeros(size(I));
mask = roipoly(J,path(:,1),path(:,2));
segmentedI = I .* mask;
figure;
imagesc(segmentedI);
colormap gray(256);
title('Segmented image');


Genetic or evolutionary algorithms for filter design? For what types of filter design problems might genetic or evolutionary algorithms be useful?

What kinds of genetic or evolutionary algorithms are used for DSP problems?

Edit: I expanded the question to include the larger set of evolutionary algorithms, such as differential evolution.

Application of a breeder genetic algorithm for system identification in an adaptive finite impulse response filter: Used a GA to determine the parameters of an adaptive FIR.

Optimizing the design of IIR filter via genetic algorithm: Used a GA to determine the parameters for an IIR based on $H_2$ and $H_{\infty}$ criteria.

A genetic algorithm for optimisation of linear phase FIR filter coefficients: (from the abstract) "The algorithm generates a population of genomes that represents the filter coefficients and compares the amplitude response of each genome to that of the desired amplitude response."

Based on this, it seems like it was a popular way to generate the coefficients for filters in the late 90's early 2000's, but there doesn't seem to be a lot of recent activity on it in the IEEE pubs.


Real Time Object Tracking using image processing How we can track a single selected object from multiple detected moving objects using single fixed camara?

What kind of an object?  If it is a compact object moving in some predictable fashion, you can use the Kalman filter to track based on its motion.  If the object has some texture, and doesn't move very fast relative to the frame rate, you can use the KLT (Kanade, Lucas, Tomasi) algorithm to track points on it.  If it has a distinctive color pattern you can use a Mean-shift tracker to track it.


Pose from Homography matrix - Zhang method This post also contributes to the post at Step by Step Camera Pose Estimation for Visual Tracking and Planar Markers by Jav_Rock (since I cannot add any comment there and I don't know why)

It can be seen that the translation vector 1x3 and rotation matrix 3x3 can be derived from homography matrix.
However, the following question is:  Where are the camera coordinate system and object coordinate system and how are they attached to the camera (or object) ?
There is the relative transformation between the two but the computation from homography or transformation matrix implies nothing about these coordinate systems' location & direction

Then, how to solve the pose estimation problem ?

The coordinate system you choose is completely arbitrary, as no information about real-world coordinates can be inferred. From an image of a table there is no reason to know that one leg is located at any particular $(X, Y, Z)$, or that it is any particular size (you can't tell if it's a doll's table or a giant's table).

Normally you would choose one of your cameras to be located at the origin, looking down the $z$ axis, defined by the matrix:

$$[R|t] = \begin{bmatrix} 1&0&0&0 \\ 0&1&0&0 \\ 0&0&1&0 \end{bmatrix}$$

Then due to the scale ambiguity you would have to choose an arbitrary scale, for example if you are using a stereo camera you could set the distance between the cameras to be one unit distance.


What are the optimal features for text line and table row discrimination in a handwritten document image? I have a handwritten document which contains graphics, tables, and text paragraphs.
I know how to extract graphics from this document. Now I want to know how can I discriminate between text lines and table rows, if I have already segmented the document into lines?

You can see an example of a document image below: in red is the graphic part. That which is framed in blue represents the table row, but there is no ruling line in the table. I want to discriminate between line text and table row using line layout.

You can detect the table from a normal text paragraph using the Hough Transform by which you can differentiate the text from the image. You can also use the bounding box method to detect the table and differentiate from the normal text rows.


Detecting a filled glass object I've been sent here from this question in stackoverflow, please excuse me if the question comes too specific and it's not in the manners in here:)

The task is to find a glass with specific liquid in it. Let me show you the pictures and then describe what i'm trying to achieve and how i was trying to achieve so far in the description below the pictures.

The pictures: (seems i need at least 10 reputation to post pictures and links, so links will have to do :( otherwise you can look at the stack overflow question)

A detailed description:
I was trying to implement an algorithm that would detect a glass of a specific shape in opencv (glass may be transformed by a different camera shot angle/distance). There will be also other glasses of other shapes. The glass i'm searching for will also be filled with some colored liquid that will distinguish it from glasses containing other colors.

So far, i have tried using SIFT feature extractor to try to find some features in the glass and then match them with other photos with the glass in it.

This approach worked only in very specific conditions where i would have glass in a very specific position and the background would be similar to the learning images. The problem also is that the glass is a 3d object and i don't know how to extract features from that (maybe multiple photos from different angles linked-somehow?).

Now i don't know what other approach could i use. I have found some clues on this (here https://stackoverflow.com/questions/10168686/algorithm-improvement-for-coca-cola-can-shape-recognition#answer-10219338 ) but the links seem to be broken.

Another problem would be to detect different "levels of emptiness" in such glass, but i haven't even been able to find the glass itself properly.

What would be your recommendations on the approach in this task? Would it be better to use a different way to find the local 3d object feature? Or would it be better to use other approach altogether? I have heard about algorithms "learning" the object from a set of multiple photos but i have never seen this in practice.

Any advice would be really appreciated

The paper referenced in your link seems to be this one.

Of particular interest there is Table 1 (included below). The accuracy rates aren't great, though they are better than other approaches.


Power spectral density interpretation After reading this question: PSD (Power spectral density) explanation I am still a little confused as to what extra information the PSD gives us over simply taking the magnitude of the fourier transform.

So my understanding is that the PSD is equal to:

$$|X(\omega)|^2$$

but what does that tell me that simply

$$|X(\omega)|$$

does not. I am assuming that squaring the fourier transform's magnitude somehow makes the number more physically (or mathematically) relevant? So how do I interpret the values of the power spectral density, especially when the signal is not a physical one and hence has no actual energy (say for example a financial signal or an image)?

Also, does the power spectral density contain information about the probability of the frequencies in the signal? As in is there a way to use the power spectral density to make a statement about the significance level of making a claim about a harmonic frequency in a stochastic signal? So something kind of like checking the significance of the coefficients from an AR model?

it doesn't tell you new information, but because the abs value function is non-analytic (i.e. not all of its derivatives are continuous), and the magnitude-squared is, then the latter can be manipulated mathematically in ways that the former cannot.  one important property of the power spectrum is that it is the Fourier transform of the autocorrelation function in the time domain.

the PSD gives you the energy of each frequency component for a deterministic signal and the expectation value of the energy for each frequency component for a random (or "stochastic" if you like jargon) signal.


FFT implementation with decimation in frequency vs decimation in time? When using FFT implementations which would be the advantages of using decimation in time (DIT) vs decimation in frequency(DIF) and vise-versa?

I have read here that in practice there is not much difference between DIT and DIF FFTs, but the hardware of some DSPs (Digital Signal Processors) is optimized for certain FFT variants. I would like to have more information about this if possible.

Well, that article equates "FFT" with "Radix-2 FFT", which isn't really true.

Anyway, since both decimation in frequency and decimation in time need exactly the same components running at the same speeds, the only difference I'd expect for radix-2 FFTs would be if you can have a parallel input on one side, and serial output on the other (or vice-versa); in that case, you could do the "large" butterfly on the side where you need to have the data in parallel to avoid latency, and could save on components by re-using a smaller butterfly for the serial end without much of a downside. Whether you want that or not very much depends on space, power, latency and throughput constraints.

For mixed-radix FFTs, where you're not only dealing with a single prime factor in the FFT size, you might not have the freedom to rearrange these layers.


How to remove stains from images? I have this extremely distorted and stained image

Is it possible to remove this stain ? Could image inpainting help ?

Please help

EDIT : Another image

After applying anisotropic diffusion and representing image with imagesc(MATLAB)

I tried inpainting however the result isn't good enough

Is there anyway i could improve this output?

This is not a complete and crisp answer however, i am leaving you with at least some approach for you to fight with. (I would be very glad to know if you have results).

Take a look at these questions:

Removing Glare from Image
How to remove the glare and brightness in an image (Image preprocessing)?

They are essentially trying to solve the same problem.

There are two parts of the problem,

a. Identifying the spot/stain
b. Replacing the stain with what would have been in the place of occlusion.

The nature of the question is trying to solve exact problem (in some sense).

This is not trivial thing. However, in both questions there are some unique pattern that you can exploit.

In all cases, the superimposing element which is required to be removed called here as  (stain, glare, bright spot), overlay has a unique and distinct hue/color which distinguishes itself from regular objects/scene.
In most cases, this hue/color of the overlay fades away into the regular scene. The actual resultant color does change - however, it is better to model overlay as single intensity and color with successively reducing transperency Hence you can say resultant pixel
$$P[x,y]  = (1-\alpha[x,y])*S[x,y] + \alpha[x,y] * OverlayHue $$
$$\tilde S[x,y] = (P[x,y] - OverlayHue * \tilde \alpha[x,y])/(1-\tilde \alpha[x,y]) $$
where $P[x,y]$ is observed image, and $S[x,y]$ is desired occlusion free image.
Note, that alpha can be arbitrarily varied over pixels, but $OverlayHue$ of overlay is considered almost constant. $\tilde S[x,y]$ and $\tilde \alpha[x,y]$ are estimated values by your algorithm for the respective quantities.
The OverlayHue value can be independently estimated by manually segmenting pixel regions where Stain or Flash is clearly dominating.
You can assume that \alpha[x,y] is consistent across all channels (i.e. R,G,B) . Hence you can identify individual components as follows:
$$\tilde S_R[x,y] = (P_R[x,y] - OverlayHue_R * \tilde \alpha[x,y])/(1-\tilde \alpha[x,y]) $$
$$\tilde S_G[x,y] = (P_G[x,y] - OverlayHue_G * \tilde \alpha[x,y])/(1-\tilde \alpha[x,y]) $$
$$\tilde S_B[x,y] = (P_B[x,y] - OverlayHue_B * \tilde \alpha[x,y])/(1-\tilde \alpha[x,y]) $$
You can see that when $\alpha$ is close to 1 implies that overlay is completely occluding the scene and hence no estimate of $\tilde S$ can be good, you should avoid that and keep some reference value to iterate it over time.
You still have more variables than equation unfortunately, and that is due to the physical nature of pixels. A given color could have resulted either due to the pixel's own property or due to stain/glare. Best bet is that you start with identified pixels where you know that $\alpha$ is 1 and then gradually decay down guess for $\alpha$ reducing successively. Over some iteration you can find the patterns.
Also, in order to finally estimate pure black spots, You can also apply smoothing constraint over neighborhood pixels (i.e. $\tilde S[x,y]$ as well as $\tilde \alpha[x,y]$).

This may not be perfect solution, but may be better than most obvious than pixel level clipping or playing around with saturation etc. I sincerely request you do try this in your end and show us results (my workbench is currently in a mess so i couldn't do it!)

Hope this helps.


Determining statistically the bandwidth of a signal from a spectrogram My question is how to determine the parameters of a signal from a spectrogram statistically.  The most important parameters that I am looking for are the lower and higher frequencies of the signal.  I have hundreds of these audio clips and I would like to perform statistics on them to determine the bandwidth of the signal. What is the best method to use to obtain these values statistically?
Find attached an audio clip and the label of signal that I am trying to analyse.

A spectrogram is a tool to analyze non-stationnary signals, and hopefully identify or separate it onto simpler sources. Hence, parts of the spectrograms ought to be "different". Hence, from a statistical point of view, the 2D data is likely to be multimodal, and global statistical measures are likely to be useless.
If you now focus on specific "modes", like identifiable whistles in boxes, then you can retrieve 2D statistical spectral moments from those patches. These could be a 2D center of mass (in time and frequency) and the corresponding 2D spread (equivalent to a standard deviation) that would form an ellipsis around the center of mass. The vertical axis could provide an idea of the frequency span of this mode.
From my experience, this often requires to look at other statistical moments for more robustness.


Computing the decimation ratio between two m-sequences Let's suppose I have an LFSR that generates an m-sequence $y_1[k]$ --- in other words, the LFSR has $N$ bits and $y_1[k]$ has period $m=2^N - 1$.

Now suppose I know someone has decimated this and taken every $j$th element, so $y_2[k] = y_1[jk+b]$. And I know the $y_2[k]$ but the values of $j$ and $b$ are unknown.

Is there any efficient way to compute the decimation ratio $j$ without testing all possible values of $j$?

(Note 1: if $b$ is unknown then there are $N$ possible solutions for $j$ because $2j, 4j, 8j, \ldots\, j^{2^{N-1}}$ can have sequences generated by LFSRs with the same coefficients, differing only by a time shift. I have no idea if the problem is easier to solve for known $b$.)

(Note 2: This smells vaguely like a discrete logarithm problem but I cannot figure out how to restate it as such. I will be satisfied if I can obtain an equivalent discrete logarithm problem, because the values of $N$ that I deal with are small enough that I can use established techniques like Pollard's rho and kangaroo algorithms to compute the value of $j$.)

There are several equivalent ways of stating this problem (I posted one on https://math.stackexchange.com/questions/2480132/determining-decimation-ratio-given-characteristic-polynomials-of-quotient-rings):

$y_1[k] = \operatorname{Tr}_{F_1}(w_1x^k)$ where $F_1 = GF(2)[x]/p_1(x)$ and similarly for $y_2[k], F_2, w_2, p_2(x)$; the traces are related (because of the decimation ratio $j$) by $\operatorname{Tr}_{F_1}(w_1x^{jk+b}) = \operatorname{Tr}_{F_2}(w_2x^k)$; the coefficients of $p_1(x)$ are the LFSR coefficients; $p_2(x)$ can be determined from the Berlekamp-Massey algorithm, and $w_1$ and $w_2$ can be determined from state recovery methods by running the LFSR backwards (I have a blog article in progress on this but it's unfinished and not a quick thing to summarize, sorry)
$p_1(x)$ is the minimal polynomial of $x$ in $F_1 = GF(2)[x]/p_1(x)$; $p_2(x)$ is the minimal polynomial of $x^j$ in $F_1$ and $p_2(x)$ is also the minimal polynomial of $x$ in $F_2 = GF(2)[x]/p_2(x)$; given $p_1(x)$ and $p_2(x)$, what is $j$?

Example:

$p_1(x) = x^9 + x^4 + 1$
$y_1[0..18] = [0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1]$
$j = 107, b = 80$
$y_2[0..18] = [1, 0, 0, 1, 0, 0, 0, 1, 1, 0, 1, 0, 1, 0, 0, 1, 1, 0]$

We can use Berlekamp-Massey on the $y_2$ values to determine $p_2(x) = x^9+x^7+x^5+x+1$, and state recovery techniques to determine that $y_2[k] = HB(v_2x^k)$ with $v_2 = x^8 + x^6 + x^5 + x^4 + x^3$ (equivalently the LFSR has initial state 101111000), where $HB(u)$ is the high bit = coefficient of $x^8$.

But I'm stuck on how to get from there to $j=107$ or other members of its cyclotomic coset $\{107, 214, 428, 345, 179, 358, 205, 410, 309\}$, aside from checking every value of $j$ --- or at least, the smallest coset representative of every value of $j$, and then using Berlekamp-Massey until I get a match of $p_2(x)$.

I finally asked this question on mathoverflow.net and got a very helpful answer.

The short answer is to follow essentially the following procedure:

Find any root $y_1 \in F_1$ of $p_2(y)$ in the polynomial ring $F_1[y]$ --- these are the polynomials $p(y)$ which have coefficients that are elements in $F_1$ and are themselves polynomials.
Take the discrete logarithm to determine $j$ such that $x^j = y_1$
Any decimation ratio in the cyclotomic coset of $j$ is a valid decimation ratio.

I've written up a longer description in my blog post which includes the use of a factoring method that partitions $p(y)$ as $p(y) = \gcd(p(y),S(x^ky))\gcd(p(y),S(x^ky)-1)$ where $S(u) = u + u^2 + u^4 + u^8 + \ldots + u^{2^N-1}$, which is mentioned in Lidl + Niederreiter and also in http://blog.fkraiem.org/2013/12/01/polynomial-factorisation-over-finite-fields-part-3-final-splitting-cantor-zassenhaus-in-odd-characteristic/.


why does the additive synthesis method for a triangle wave require amplitude scaling by 8/pi^2? I had to make a bunch of band limited digital triangle waves recently, so I went to (where else) wikipedia for the equations.
I noticed that there is a constant amplitude scalar of 8/pi^2 (~.8106). And indeed this does bring the amplitude range down to about +/-1. The wikipedia article (https://en.wikipedia.org/wiki/Triangle_wave#Harmonics) doesn't explain the theoretical reasoning for this but I would like to know why this is so.

Are you familiar with how to obtain the Fourier series of an arbitrary periodic function? Because if you are, I'm confused at your confusion.
This is simply the Fourier series expansion of the triangle wave with fundamental frequency $f_0=1/T$. This can be found by taking the inner product of the function with the basis functions. We note that the DC value (the time-average) is zero, and that the function is odd, so it has no cosine terms. So, let's attempt to expand it as the sum of sines or in other words as
$$ x(t) = \sum_{k=0}^{\infty} b_k \sin\left(\frac{2 \pi k t}{T}\right) $$
where
$x(t)$ is the triangle wave defined on the interval $-T/2$ to $T/2$ as
$$ x(t) = \left\{\begin{array}{cc}
\frac{4 t + 2 T}{T} & -\frac{T}{2} < t < -\frac{T}{4} \\
\frac{4 t}{T} & -\frac{T}{4} < t < \frac{T}{4} \\
\frac{4 t - 2 T}{T} & \frac{T}{4} < t < \frac{T}{2}
& \end{array}\right. $$
and elsewhere by the fact that $x(t) = x(t+T)$.
To calculate our Fourier series coefficients, we must compute the integral
$$ b_k = \frac{2}{T} \int_{-\frac{T}{2}}^{\frac{T}{2}} x(t) \sin\left(\frac{2 \pi k t}{T}\right) dt.$$
Because the function we are integrating is a product of two odd functions, it is even, so we can reduce this to
$$ b_k = \frac{4}{T} \int_{0}^{\frac{T}{2}} x(t) \sin\left(\frac{2 \pi k t}{T}\right) dt.$$
We notice that $x(t)$ exhibits even symmetry about the point $\frac{T}{4}$ while even frequency sine waves of the form $\sin(4\pi k t)$ exhibit odd symmetry about the same point. This is quarter wave symmetry and tells us that the even harmonics will all be zero. So, we only need to solve this for odd values of $k$. To that end, we replace $k$ with $2 k + 1$.
We can once again exploit the even symmetry of the odd harmonics to further reduce the range of our integral to
$$ b_{2 k + 1} = \frac{8}{T} \int_{0}^{\frac{T}{4}} \frac{4 t}{T} \sin\left(\frac{2 \pi (2 k + 1) t}{T}\right) dt,$$
and, skipping the calculus (for brevity), we arrive at
$$ b_{2 k + 1} = \frac{8}{\pi^2 (2 k + 1)^2} \cos(\pi k). $$
This reduces to the formula given on Wikipedia. If you don't see it, let me know and I'll be explicit.


As of 2019, which discrete nonlinear, time-invariant systems with memory are considered "easy" to model and identify? There are several types of discrete nonlinear time-invariant systems with memory ("NTIM") which are considered "easy" to model and identify. Any such system can be represented using a Volterra series, but Volterra kernels become very difficult to model for any more than quadratic or perhaps cubic nonlinearities.

However, for certain types of NTIM systems, there exist clear methods to completely measure and replicate some arbitrary unknown system, typically via measurement of certain test signals.

One good example are those systems for which each Volterra kernel has nonzero coefficients only on the diagonal. These are sometimes called "Farina," "diagonal-Volterra," or "Generalized Hammerstein" systems (Fig. 1). Given some suitable smoothness assumptions, these can be perfectly modeled using dynamic convolution, which is a variant of regular convolution that has a different impulse response for each scaled impulse. As a result, an arbitrary system of this type can be perfectly characterized by simply measuring the impulse response for each scaled impulse. Some info on that in this related answer: Dynamic convolution vs Volterra series. (There are also techniques involving synchronized sine sweeping)

Figure 1. A diagonal Volterra series system, which has nonzero Volterra coefficients only on the diagonal of each Volterra kernel, is equivalent (given suitable smoothness assumptions) to a sum of cascades of a memoryless waveshaper (illustrated as a sigmoid) and a linear time-invariant (LTI) system.

I am wondering if there exists any larger set of nonlinear systems that is also considered easy to model. There also exist the Wiener (Fig. 2) and Wiener-Hammerstein systems (Fig. 3) , which are not diagonal-Volterra; are those easy to model? And all of these systems seem to be subsets of "NARMAX" systems; how easy are they to model? Is there any easily modelable type of system which generalizes both the Weiner-Hammerstein and diagonal-Volterra systems?

Figure 2. A Wiener system is a cascade of an LTI system and a memoryless waveshaper.

Figure 3. A Wiener-Hammerstein system is a cascade of an LTI system, a memoryless waveshaper, and another LTI system.

Questions

As of 2019, which nonlinear time-invariant systems with memory are clearly easy to model and profile?
Is there any type of nonlinear system, which includes all the Volterra-diagonal ones as a subset, which is easy to model?
Is there some type of system which includes both the Volterra-diagonal and the Weiner-Hammerstein systems, and is easy to model?
All of these systems are a subset of the "NARMAX" nonlinear models - are those easy to model?

It is very strange phenomena that one object is completely dropped out of attention of researchers. It is Urysohn operator. First of all Urysohn is equivalent to multiple parallel Hammersteins and Urysohn followed by static nonlinearity is a model of any deterministic dynamic object, it maps any given input to any provided output. I obtained Ph.D. in modelling dynamic systems as Urysohn operators in 1991. Since that time I showed how to do that to many people, I published articles, I maintain the site, explaining it, and for this almost 20 years noone started using it, there is no single reference on my articles. People write 10,000 lines to identify Hammerstein, I show how to identify Urysohn by 10 lines of code, write here. Look

for (int i = T - 1; i < N; ++i) {
   double predicted = 0.0;
   for (int j = 0; j < T; ++j) {
      predicted += U[(int)((x[i-j] - xmin) / deltaX), j];
   }
   double error = (y[i] - predicted) / T * learning_rate;
   for (int j = 0; j < T; ++j) {
      U[(int)((x[i-j] - xmin) / deltaX), j] += error;
   }
}

That is it. It replace 10,000 line Hammerstein identification. It published, it showed, why noone use it I don't know. Google show my articles on "Urysohn identification" at the top. People don't search.

Explanation to code. x[i], y[i] elements of input output arrays. U[k,m] elements of model matrix. T is one size of matrix U, another one is determined by deltaX, discretization interval, learning_rate is similar to LMS.


plotting phase transformations of a signal I would like to find the phase variations of the signal with respect to time axes:

S = sin(2*pi*100*t + pi/4) + cos(pi*500*t) + sin(2*pi*100*t + 5*pi/2);

My attempt to get the phase, see figure below:


Confusion in basics of Laplace Transform I have few confusions while starting Laplace Transform. So far I have studied, Fourier series and Fourier Transform. The basic difference which I found from different books is Fourier Transform is only considered the imaginary part whereas the Laplace transform considers both real and imaginary for general values.
i) I want to ask that, is it only difference we have in Laplace and Fourier Transform?
Then i saw the two different equations of Laplace transform which is bilateral Laplace Transform $$ X(s) = \int _{-\infty}^{+\infty} x(t) e^{-st} dt$$
Whereas in the second equation of Laplace Transform is called unilateral Laplace Transform and it defined as:
$$ X(s) = \int _{0}^{+\infty} x(t) e^{-st}dt $$
It omit the negative part and only have for $t>0$

ii) Here I want to ask that what is the reason of omitting the $t<0$ part?

And Lastly, there was an example  $$ x(t)=e^{-at}u(t)$$
After applying Laplace transform it was written The transform exists only if $Re(s+a)$ is possitive.
iii) Now here is am confuse that why it took for possitive?

The unilateral Laplace transform is used for analyzing causal linear time-invariant systems, which have an impulse response $h(t)$ that is zero for $t<0$. The unilateral Laplace transform can be used to solve initial value problems, due to the correspondence

$$x'(t)\Longleftrightarrow sX(s)-x(0)$$

where $x(0)$ is a given initial value for the function $x(t)$. Note that for the bilateral Laplace transform the equivalent correspondence is simply $x'(t)\Longleftrightarrow sX(s)$.

Concerning the signal $x(t)=e^{-at}u(t)$, note that its Laplace transform is

$$X(s)=\int_0^{\infty}e^{-(a+s)t}dt$$

This integral only converges if the exponential decays, which results in the condition $\text{Re}(a+s)>0$. This condition defines the region of convergence (ROC) of the Laplace transform. $X(s)$ only exists for values of the complex variable $s$ satisfying $\text{Re}(s)>-\text{Re}(a)$.


start of packet detection With the current lockdown, I'm spending some time in understanding L1 PHY development  and have some questions on the start of packet detection in an IQ capture

Let's say I have a WiFi IQ capture which consists of 1000 packets and want to extract indexes of start and End of each WiFi Packet.

When I plot abs(IQ_samples), it is like

WiFiPacket|NULL|WiFiPacket|NULL|NULL|NULL|WiFiPacket|NULL|WiFiPacket ....

Since the energy  of NULL is very much close to zero compared to the max_value of abs(IQ_sample of each packet). Can I keep some kind of threshold to find the start and end of the packet.

This is what I found after searching over internet and going through some posts
*) find the minimum value of IQ capture .. that would give me noise value (NVAL)
*) take some threshold which is a factor of NVAL (Threshold = factor x NVAL)

My main confusion is what should be the factor value or how should I find that?

Let us say you detected noise level power at -100dBm. You keep additional 10dB threshold to detect energy. Means if you detect >= -90dBm power, you try to detect preamble of Wifi packet. If you detect -95dBm power, you ignore it assuming it may not be a valid Wifi packet. This really works assuming you are very far from base station/access point(AP). Now suppose say you want a higher threshold so that you don't want too many false alarms (may be in order to save battery). So you raise the threshold to -60dBm. Here, if you are too far from AP, you may not trigger at all. If you are close to say within 10 meters, you start to detect energy more than -60dBm and try to detect Wifi packets.

So, you can decide on threshold on how frequently you want to try to detect Wifi signal.

If your threshold is too low, you may more frequently try to detect packets. Advantage - you can detect packets if you are far away from AP. Disadvantage- you may consume more power trying to detect packets.

If you increase threshold, you will detect packets less frequently. Advantage - you can filter out signals which are weak and can use higher modulation rates. Disadvantage - you cannot move too far away from the AP.


People detection from above I am trying to find some method to detect people using only one camera 3 meters above the ground. This is a frame returned by camera:

UPDATE: Video test -> http://dl.dropbox.com/u/5576334/top_head_shadow.avi

In order to do that, first I understand that I have to perform a background-foreground segmentation. That is the easy part.

With the foreground mask, I am able to make simple operations such Hough transform to find circles, but this way only detects the 60% of heads, including many false positives.

I could use some other simple techniques like color segmentation, but I found that people heads are very different seen from above because of their hairstyle, color, amount of hair,...

Other option I have though about it is the possibility of using HOG Descriptors, or Haar-like features, but I would need an extensive database of people seen from above to train the models. I have not found anything like that.

I thought this would be a very recurrent problem, but I can not find very much about it in the literature or internet. Any help to resolve this task will be appreciated :-)

UPDATE: For more information, The goal is to implement some generic method to make pedestrian flow tracking. The first prototype will be tested in a Mall.

Are there any other objects that can move beside people? If there aren't, you can just find the  blobs (connected components) in your foreground mask, and these are your people.

They can also "collide" one with another, creating one blob instead of two. In this case, you can do a motion tracking and resolve the ambiguity by using the fact that the trajectory and the speed is smooth.

If there are other objects (like dogs, cars), you should create a classifier that gets blob parameters like:

Blob statistics (size, solidity,etc..)
Color
Edge information
Speed (in case of tracking)

And returns the correct class (Human/No Human).


Evaluating the inverse $\mathcal Z$-transform on the unit circle I am trying to understand the math.  The inverse $\mathcal Z$-transform is given by:

$$x[n] = \displaystyle\frac{1}{j2\pi} \int_cX(z)z^{n-1}dz$$

where $\displaystyle \int_c$ is a contour integral. The inverse Fourier transform is given by:

$$x[n] = \displaystyle\frac{1}{2\pi}\int_{-\pi}^\pi X(e^{j\omega})e^{j\omega n}d\omega$$

My textbook claims, and I agree, that evaluating the inverse z-transform at $z=e^{j\omega}$ will result in the inverse Fourier transform.  However, I can't get the math to show this. Substituting $z=e^{j\omega}$, I get as far as:

$$x[n] = \displaystyle \frac{1}{j2\pi}\int_cX(e^{j\omega})e^{j\omega n}e^{-j\omega}d(e^{j\omega})$$

Can somebody tell me intuitively how to simplify this equation to look like the inverse Fourier transform? Is my confusion due to a lack of understanding of the contour integral?

Hint:

With $z = e^{j\omega}$, you have:

\begin{align}
\displaystyle \frac{dz}{d\omega} &= je^{j\omega}\\
\Rightarrow dz &= je^{j\omega} d\omega
\end{align}

With $-\pi \leq \omega \leq \pi$, your contour integral goes from $-\pi$ to $\pi$.

  Can somebody tell me intuitively how to simplify this equation to look
  like the inverse Fourier transform?

Plugging all the above should do the job.

  Is my confusion due to a lack of understanding of the contour
  integral?

It is more the change of variables.

EDIT:

Note that $\omega$ is taken in the fundamental interval $\left[-\pi, \pi\right]$ because $X(\omega)$ is periodic with period $2\pi$. So, $\displaystyle \int_c$ becomes $\displaystyle \int_{2\pi}$


How to Generate False Color Palette How to convert grayscale intensity input to something like the following image?

It seems that no response (0, black) is converted to dark blue, then it goes through other colors to red.

How are these colors derived? Is there any standard gray-to-RGB conversion function for these?

I have tried representing intensity as hue in CIE-LCH model and then converting to RGB through CIE-Lab. But the result weren't as nice as these MRI and PET scans.

Another approach may be a rainbow, i.e. going from high frequency (violet-blue) to lower frequencies (dark red). But is there a conversion formula for this?

What you're looking for is a color map. As its name implies, this process maps an input (single-component, or grayscale) value within a predetermined range to some other color based upon some mapping. The implementations I've seen typically use linear RGB interpolation between a number of control points specified throughout the allowed input range. The color map used in the image you showed is a pretty common one, starting with dark blue, then through cyan, green, yellow, orange, and ending with red. It is referred to as the jet colormap in MATLAB.

So, as an example, for an input intensity on the range $[0,1]$, you might have an RGB lookup table of:

0.0  -> (0, 0, 128)    (dark blue)
0.25 -> (0, 255, 0)    (green)
0.5  -> (255, 255, 0)  (yellow)
0.75 -> (255, 128, 0)  (orange)
1.0  -> (255, 0, 0)    (red)

Then, if you wanted to map the value $0.2$ to a color, then you would linearly interpolate between the first two datapoints as follows:

color = 1 / 5 * (0, 0, 128) / 0.25 + 4 / 5 * (0, 255, 0)
      = (0, 204, 26)

which would yield a greenish cyan, as it's close to the green control point.


Relationship between entropy and SNR In general any form of enropy is defined as uncertainty or randomness. In a noisy environment, with increase in noise, I believe that entropy increases since we are more uncertain about the information content of the desired signal. What is the relationship between entropy and SNR? With increase in signal to noise ration, the noise power decreases but this does not imply that the information content of the signal increases!!
The information content may remain same, so does that mean entropy is unaffected?

When you say that the "information content may remain the same," do you mean the information in the total signal, or the information of the desired signal? Hopefully this will answer both cases. I know Shannon entropy much better than Kolmogorov so I'll use that, but hopefully the logic will translate.

Let's say $X = S + N$ is your total signal ($X$), comprised of the sum of your desired signal $S$ and your noise component $N$. Let's call entropy $H$. As you said, noise adds entropy to a system by increasing its complexity. However, it's not necessarily only because we're more uncertain about the information content of the signal, but because there's more uncertainty in the signal overall. If SNR kind of measures how certain we are of what $S$ is, then $H(X)$ kind of measures how well we can predict future states of $X$ based on the current state of $X$. Entropy is concerned with how complex the whole signal is, regardless of the composition of noise vs. non-noise.

If you increase SNR by removing noise (attenuating $N$), you decrease the total signal $X$'s complexity and thus its entropy. You haven't lost any information carried by $S$, only (presumably meaningless) information carried by $N$. If $N$ is random noise, then obviously it doesn't carry meaningful information, but it takes a certain amount of information to describe $N$'s state, determined by the number of states that N can be in, and the probability of it being in each of those states. That's the entropy.

We can look at two Gaussian distributions with different variances, say one has a variance of $1$ and the other has a variance of $100$. Just looking at the equation for a Gaussian distribution, we see that the $Var=100$ distribution has a maximum probability that is only $\frac {1} {10}$th the value of the $var=1$ distr's probability. Conversely, this means that there is greater probability that the $Var=100$ distr will take values other than the mean, or that there is more certainty that the $Var=1$ distribution will take values near the mean. So, the $Var=1$ distribution has a lower entropy than the $Var=100$ distribution.

We established that higher variance implies higher entropy. Looking at error propagation, it is also true that $Var(X+Y) >= Var(X) + Var(Y)$ (equal for independent $X$, $Y$). If $X = S + N$, then for entropy $H$, $H(X) = H(S+N)$. Since $H$ is (indirectly) a function of variance, we can fudge things a little to say $H(Var[X]) = H(Var[S+N])$. To simplify, we say $S$ and $N$ are independent, so $H(Var[X]) = H(Var[S] + Var[N])$. Improved SNR often means attenuating noise power. This new signal with higher SNR will then be $X = S + (\frac {1} {k})N$, for $k>1$. Entropy then becomes $H(Var[X]) = H(Var[S] + (1/k)^2*Var[N])$. $k$ is greater than $1$, so $Var[N]$ will decrease when N is attenuated. If $Var[N]$ decreases, so does $Var[S+N]$, and therefore $Var[X]$, resulting in a decrease in $H(X)$.

Not very concise, sorry. In short, $X$'s entropy decreases if you increase SNR, but you've done nothing to $S$'s information. I can't find the sources right now, but there's a method to calculate SNR and mutual information (a bivariate measure similar to entropy) from each other. Maybe the main takeaway is that SNR and entropy don't measure the same thing.


How to covert hummming into pitches? I am not sure if this question suits this site, but I want to convert humming into notes or relative values to analyse melodies, word of "Signal Processing" looks suit for this. I am planning to use Flash AS3, and Tortoise and Hare for Cycle Detection. I need very simple script implementation. All I want is find notes more easily with my bass. Is that possible? Or is there any frameworks that already makes this?

You can try using a Pitch Track for monophonic signal, you need track low pitches!

There is a wide range of algorithms to find Pitches in time domain (zero crossing, AMDF, YIN, MPM, ASDF, autocorrelation ....) or Frequency domain (shs, hps, cepstrum ....)

Each one of them possess advantages or disadvantages !

With the Pitch Track you can convert one frequency to note.

I never used Flash AS3 to DSP, one fast and easy way are the autocorrelation based algorithms.


Contours chains matching techniques I have chain contours representation:
- classic Freeman chain
- angle chain with distance (first: angle between neighbor points, second: distance between   them)

What are the best techniques to match them, and If it is possible figure out rotation of element?

For now I tested approach with FFT(isn't it overuse for short chains?), shifting and comparing
For Freeman chains Pairwise Geometrical Histogram, but it figures out rotation only with 45 degree resolution.

Any sugestions, codes, references will be appreciated.

I tried http://www.codeproject.com/Articles/196168/Contour-Analysis-for-Image-Recognition-in-C once. It did not work well for me, but that guy seem to have better results.


Relation between Gaussian derivatives and Gausian-windowed Cosine function I am doing some research on UWB radars which transmit a pulse given by:
$$
g(t) = p(t) cos(2\pi f_{c}t) = A e^{(-\frac{t^{2}}{2\tau ^{2}})} cos(2\pi f_{c}t)
$$
In some other papers, the transmitted pulse is also described as a higher order derivative of the Gaussian term $p(t)$, such that:
$$
g(t) = \frac{d^{n}}{dt^{n}} \left( Ae^{(-\frac{t^{2}}{2\tau^2})} \right)
$$
The Gaussian derivatives can be found recursively, which form the Hermite polynomials. Is there a way of showing how the two are equivalent analytically, i.e. without using graphical methods? I am aware of the fact that Hermite polynomials have asymptotic expansions in terms of a cosine function, but this does not help me for my derivation.

You may be talking about a sine wave (carrier) of frequency $f_c$ modulated by an Hermite polynomial shape function. If a parameter $1/r$ of your UWB pulse waveform is commensurate with a carrier frequency $f_c$, this setup can also be considered an ultra wide band scenario. But if this is the case, the $\cos(2πf_ct)$ wave and the waves, which higher order Hermite polynomials resemble at low values of their argument, are absolutely unrelated things.
The other possibility is that you are talking about an Hermite polynomial approximation. As you are aware of the Hermite polynomial's asymptotic expansions, you should know that the cosine frequency $f_c$ in your first equation is not arbitrary (it must not be) and you should not write down a cosine of independent frequency $f_c$ in your waveform. The approximation of an n-th order Hermite polynomial is written as
$$
\exp(-{\frac {x^{2}}{2}})\cdot H_{n}(x)\sim {\frac {2^{n}}{\sqrt {\pi }}}\Gamma \left({\frac {n+1}{2}}\right)\cos \left(x{\sqrt {2n}}-{\frac {n\pi }{2}}\right) \tag{1}
$$
Because $H_n$ has n roots, eq.1 can be valid only in a bounded range of $x$ values. If you are tolerant of what may appear a circular-reasoning logical fallacy, I refer you to an improved approximation of $H_n$ in the Wikipedia article for you to be able to estimate the bound for $x$:

A better approximation, which accounts for the variation in frequency,
is given by
$$ e^{-{\frac {x^{2}}{2}}}\cdot H_{n}(x) \sim \left({\frac {2n}{e}}\right)^{\frac {n}{2}}{\sqrt {2}}\cos \left(x{\sqrt {2n+1-{\frac {x^{2}}{3}}}}-{\frac {n\pi }{2}}\right)\left(1-{\frac {x^{2}}{2n+1}}\right)^{-{\frac {1}{4}}}. $$

Id est, the eq.1 approximation is valid for $x^2/3 \ll (2n+1)$ .
Taking into account this constraint for $x$, you can prove that the recurrence relation $\widetilde {H_{n}}'(x)=2n\widetilde {H_{n-1}}(x)$ is asymptotically precise for eq.1's functions approximating precise Hermite polynomials $H_{n}, H_{n-1}$. The $\widetilde{H_0}$ approximation expression is a precise Hermite polynomial $H_0$, and one can conclude by induction that eq.1 holds in the indicated range of $x$ values for Hermite polynomials of all orders, with bounds expanding with each transition to higher order polynomials.
For a rigorous derivation of Hermite polynomial approximations, see, for example, The Asymptotic Forms Of The Hermite And Weber Functions By Nathan Schwid, page 358 (20/24), Theorem VIII, and Asymptotic analysis of the Hermite
polynomials from their differential-difference equation by Diego Dominici, page 13/21, Theorem 5.


Signal with Complex Gaussian noise If I have:

$$ y = x_r+jx_i + n_r +j n_i$$

with $n_r$ and $n_i$ Gaussian with mean 0 and variance $\sigma^2$, what is the pdf of the envelope |y| and phase(y)? Is it still Rayleigh-distributed and uniform, respectively?

No, it will be Rician Distribution.

Rayleigh distribution is a special case of Rician Distribution when the normal random variables involved are of zero mean and equal variances.

I have assumed you meant :$$y=x_r + j \cdot x_i + n_r + j \cdot n_i$$ where $x_r, x_i$ are deterministic variables.

Basically, you now have $y_r = (n_r + x_r),\ \sim \mathcal N(x_r, \sigma^2) $ and $y_i = (n_i + x_i), \ \sim \mathcal N(x_i, \sigma^2)$

Now, you can write $y = y_r +j \cdot y_i = r \cdot e^{j\phi}$, where $r,\phi$ are random variables and, $y_r = r \cdot \cos(\phi)$ and $y_i = r \cdot \sin(\phi)$.

Using method of transformations, Jacobian matrix can be given as follows:
$$\begin{pmatrix} y_r \\ y_i \end{pmatrix} = \begin{pmatrix} r \cdot \cos(\phi) \\ r \cdot \sin(\phi)\end{pmatrix}$$
$$\mathbf J = \begin{pmatrix} \cos(\phi)&&-r \cdot \sin(\phi) \\ \sin(\phi)&& r \cdot \cos(\phi)\end{pmatrix}$$
$$|\mathbf J| = r$$
Hence, you see that Jacoboan does not change by changing the mean of the normal random variables from 0 to non-zero.

Joint PDF of derived RVs $r, \phi$ will be given by the following :
$$f_{R,\Phi}(r, \phi) = f_{Y_r, Y_i}(r\cos(\phi), r\sin(\phi)) \cdot |\mathbf J |$$
$$f_{R,\Phi}(r, \phi) = \frac{r}{2\pi \sigma^2} \cdot \exp\left[-\frac{(r\cos(\phi) - x_r)^2+(r\sin(\phi)-x_i)^2)}{2\sigma^2}\right]$$
$$\forall r\in [0,\infty) \\ \phi \in[-\pi, \pi]$$

Solve for Marginal PDF of $R$ and $\Phi$ to see that it is not Rayleigh anymore. But a more generalized Rician Distribution.


Discrete Cosine Transform Is this piece of text (written by me) true about Discrete Cosine Transform ?
Discrete Cosine Transform (DCT) is a lossy data compression algorithm that is used in many compressed image and video formats,  including JPEG,  MJPEG, DV and MPEG. In this algorithm, special DCT coefficients  are calculated for each 8x8 image block.  Then the coefficients are quantized and the image block is represented as a matrix of these quantized coefficients.  The algorithm utilizes the fact that human visual system does not distinguish small changes in color or intensity.

To illustrate Justme's answer:

Discrete Cosine Transform (DCT) is a lossy

The DCT can't be a lossy algorithm, since there's an inverse operation that restores the original input exactly.

data compression algorithm

Also, it's not a compression algorithm: in- and output have the same size.
So, both your central statements are wrong :(

that is used in many compressed image and video formats,

yes.

including JPEG, MJPEG, DV and MPEG.

What's DV? And: MPEG is a huge family of video compression methods. There's not "the MPEG algorithm", there's dozens of different video compression standards under that name.

In this algorithm, special DCT coefficients are calculated for each 8x8 image block.

That applies to JPEG, and probably a few of the many MPEG codecs. It's not true for all MPEG compressors! (For example, MPEG-H Part II, also called H.265, uses blocks of 64×64, 32×32, or 16×16, 8×8 or 4×4, depending on the image content.)

Then the coefficients are quantized

And that's where the lossiness happens: it's not in the DCT, it's what happens to the output of it!

and the image block is represented as a matrix of these quantized coefficients.

Again, only applies to JPEG. Be clear about that!
The matrix is how it's often visualized. The matrix representation is actually non-existent in the memory or storage formats, usually. On the very contrary, the elements are typically stored in a zigzag diagonal ordering if we imagine the matrix. (that's because that puts values that tend to be correlated closer to each other, making the result better to compress using subsequently applied lossless methods like LZW, Huffmann.)

The algorithm utilizes the fact that human visual system does not distinguish small changes in color or intensity.

No, it uses the fact that human perception often cares about high-frequency changes less than about small changes in low-frequency components. Otherwise, the selective quantization would make no sense.

Hm, you've not written the greatest paragraph. But I think you understand a lot of things correctly. Be more careful in really being aware what exactly does what, and you'll be fine!


Confusion with data after cross correlation Surely this question has been answered before in parts but I cannot find one solid answer. Ultimately I am trying to calculate the audio delay between two signals I am recording on USB microphones. With no problem I am recording the audio, reading it back into python, computing the cross correlation/convolving, and then finding the delay. When I manually create the signals (one with a delay) and feed them into my algorithm both methods find the right delay. As a test for real data I am trying to find the speed of sound between two mics (velocity = dis / time). Again, with a man-made delay of 76 frames and mic spacing of .6 meters I get a speed of sound of roughly 348.157...close enough for me.

When I try this with real recorded data no two tests have similar outputs; huge outliers, and no consistency.

What should I try next? filtering?
Add more mics for more accuracy?

The source signal is Guassian white noise. Any help you guys can give would be awesome.

I had done something like this and ran into similar problems a while ago, using Matlab. The audio to the microphones was speech or whistling in my case. I was estimating the delay between two mics to find the direction sound was arriving from.

The audio data from the mics probably has oscillations of some kind (like a sum of sinusoids plus some noise), at least at certain locations. The cross-correlation of one such waveform with a delayed version of it will be a signal with oscillations. This would look like this example, but possibly decaying slowly. I've seen something similar in this question.

The cross-correlation peak should indicate the delay, but there are multiple peaks so you may be detecting the wrong one. Since your input signal is not a pure sine-wave, the largest peak should be the delay but it would be hard to guarantee that in a hardware system. This may be why you're getting inconsistent results.

Since you know the distance between the two microphones (and the system's sampling rate), you can estimate the maximum expected delay (in a number of samples), as the time it would take for sound to travel the full distance between the mics (0.6 meters). Then you should search for the cross-correlation peak only in the first such number samples, rather than in the full waveform.

Low-pass filtering may also be helpful if you're just interested in the delay. The higher frequency components can create peaks in the cross-correlation sooner (oscillations with a smaller period). In my case, I had found that a low-pass cut-off of about 1.0-1.5 kHz was adequate since my signal (speech) was wide-band and contained enough information at low frequencies for the delay estimation.


frequency components of discrete Gaussian-windowed signal When we are do a discrete Gabor transform on a 1D signal, we apply a Gaussian window on a segment of the signal. However, due to the tapering ends of the Gaussian window, the magnitude of the data points will be reduced more towards the ends, This implies that the the frequency components of the segment of the signal will be suppressed when we apply FFT to it. After performing some experiments, I found that at least 2 periods of the signals is required, so that no misjudgement of the behaviour of the signal will happen.

What is the range of frequencies in which the magnitude of the FFT is still representative of the frequencies of windowed signal?

No window will give you zero processing (coherent) gain. Check out the table in fred harris's paper.

The table says there is a processing gain of $0.37$ to $0.51$ depending on your $\alpha$ parameter:


When is a discrete time transfer function unrealizable? I don't understand why the following makes sense:

Given a second-order mass damper system in continuous time:

$H(s) = \frac{1}{ms^{2}+cs}$

Its inverse $H^{-1}(s)$ is unrealizable as a transfer function $G(s) = \frac{Y(s)}{U(s)}$ has a state space realization if and only if the degree of $Y(s)$ is less than or equal to the degree of $U(s)$.

However, the zero-order hold (ZOH) equivalent of $H(s)$ given by:

$H(z^{-1}) = \frac{z^{-1}(b_{0}+b_{1}z^{-1})}{1+a_{1}z^{-1}+a_{2}z^{-2}}$

has a realizable inverse $H^{-1}(z^{-1})$ as its degree in the denominator and numerator is equal.

How does this make sense?

I have read that when designing a disturbance observer, we need a $Q(z)$ filter such that $Q(z)G^{-1}(z)$ is realizable as $G^{-1}(z)$ will not be.

You give a transfer function in $z^{-1}$, but for the purposes of making a state-space model, you need a transfer function in $z$ (more on that later).  So your $H$ becomes
$$H(z) = \frac{b_0 z + b_1}{z^2 + a_1 z + a_2}$$
When you invert that it's improper, which shows that you need to know information one time step into the future.  That's not possible (or, if it is, you should take up life as a stock broker).

(Note that you don't need state space to know this: you can just write the corresponding difference equation for your $H^{-1}$ and you'll see that it's asking for future values of the input to affect the current value of the output.)

Appendix

The reasoning for doing it in $z$ is because the time-domain state-space realization is

$$\begin{split} x_k = A x_{k-1} + B u_k \\ y_k = C x_{k-1} + D u_k\end{split}$$

For $D = 0$, the resulting transfer function is strictly proper in $z$ (but not in $z^{-1}$); for $D \ne 0$ the resulting transfer function in $z$ has equal orders in the numerator and denominator.


Convolution and windowing using a buffer - how do I do overlap add? I'm reading chunks of audio signal (1024 samples) using a buffer, applying a Hanning window, doing an FFT on it, then reading an Impulse Response, doing its FFT and then multiplying the two (convolution) and then taking that resulting signal back to the time domain using an IFFT. I'm zero padding the Impulse Response to match the length of the input signal chunk. However, I still get a noisy click sound after it reads every chunk, because of differences in amplitude I suspect.

At this stage I want to do an Overlap Add but am unclear on how to do it on a buffer. Convolving the input signal chunk with the Impulse Response results on a signal with the same number of samples (1024) so what am I overlapping when it gets to the next chunk?

Should I be using a longer impulse response and zero pad the signal
instead, and then add the rest of that signal after sample 1024 to
the next window and so on?  Or should I use a longer window on the
input signal?

I am doing all this in Python using PyAudio.

Thank you very much for your help

normally i like to define windows to have even symmetry about $t=0$ or $n=0$, but this time i'll do it like MATLAB (except i know how to count from $0$).

Hann Window of width $B$ samples:

$$ w[n] \triangleq \begin{cases} \frac12 \left(1 - \cos\left(\frac{2 \pi}{B} (n+\frac12)\right) \right) \quad \text{for } 0 \le n < B \\
0  \quad \text{for } n < 0 \text{ or } B \le n \\
 \end{cases} $$

another way to write it is

$$ w[n] = \tfrac12 \left(1 - \cos\left(\tfrac{2 \pi}{B} (n+\tfrac12) \right) \right) \cdot \left(u[n] - u[n-B] \right) $$

where $u[n]$ is the discrete unit step function

$$ u[n] \triangleq \begin{cases}
1 \quad \text{for } n \ge 0 \\
0 \quad \text{for } n < 0 \\
\end{cases} $$

all $B$ samples within the window width are non-zero.  let's assume $B$ is even, so that $\frac{B}2$ is an integer.  turns out that the Hann window is "complementary" which means the downslope of one Hann window adds to $1$ with the upslope of the window adjacent to its right if they are staggered by $\frac{B}2$ samples.  in fact this is true:

$$ \sum\limits_{m=-\infty}^{\infty} w\left[n-m\tfrac{B}2 \right] = 1 $$

that means

$$ \begin{align}
x[n] & = x[n] \sum\limits_{m=-\infty}^{\infty} w\left[n-m\tfrac{B}2 \right] \\
     & = \sum\limits_{m=-\infty}^{\infty} x[n]  w\left[n-m\tfrac{B}2 \right] \\
     & = \sum\limits_{m=-\infty}^{\infty} x_m[n-m\tfrac{B}2]   \\
\end{align} $$

where

$$ x_m[n] \triangleq x\left[n+m\tfrac{B}2 \right]  w[n] $$

so we broke up the input into overlapping little "wavelets" (my term, not in the sense of Daubechies) or "grains", $x_m[n]$, each delayed, in the summation, by $\frac{B}2$ samples in relation to their preceding neighbor $x_{m-1}[n]$.

now we know that for the $m$-th delayed grain, $x_m\left[n-m\tfrac{B}2 \right] \ne 0$ only for $m\tfrac{B}2 \le n < m\tfrac{B}2 + B$

FIR convolution

we have an FIR filter with causal finite impulse response, $h[n]$, with length $L$.  so $h[n]=0$ for $n<0$ or $L\le n$.  the output of this FIR filter is

$$\begin{align}
y[n] & = h[n] \ \circledast \ x[n] \\
     & \triangleq \sum\limits_{i=-\infty}^{\infty} h[i] x[n-i] \\
     & = \sum\limits_{i=0}^{L-1} h[i] x[n-i] \\
\end{align}$$

OLA convolution

one way to write this:

$$\begin{align}
y[n] & = h[n] \ \circledast \ x[n] \\
     & = h[n] \ \circledast \ \sum\limits_{m=-\infty}^{\infty} x_m[n-m\tfrac{B}2] \\
     & = \sum\limits_{m=-\infty}^{\infty} h[n] \ \circledast \ x_m[n-m\tfrac{B}2] \\
     & = \sum\limits_{m=-\infty}^{\infty} y_m[n-m\tfrac{B}2] \\
\end{align}$$

where

$$ \begin{align}
 y_m[n] & \triangleq h[n] \ \circledast \ x_m[n] \\
        & = \sum\limits_{i=0}^{L-1} h[i] x_m[n-i] \\
\end{align} $$

another way to say it is

$$\begin{align}
y[n] & = h[n] \ \circledast \ x[n] \\
     & = \sum\limits_{i=0}^{L-1} h[i] x[n-i] \\
     & = \sum\limits_{i=0}^{L-1} h[i] \sum\limits_{m=-\infty}^{\infty} x_m[n-m\tfrac{B}2 - i]  \\
     & = \sum\limits_{m=-\infty}^{\infty} \sum\limits_{i=0}^{L-1} h[i] \ x_m[n-i - m\tfrac{B}2]  \\
     & = \sum\limits_{m=-\infty}^{\infty} y_m[n-m\tfrac{B}2]  \\
\end{align}$$

$y_m[n]$ defined the same way as above.

the potentially non-zero length of the block of input samples $x_m[n]$ is $B$ samples, the length of the FIR is $L$ samples, and that makes the potentially non-zero length of the block of output samples $B+L-1$.  you can see that

$$ y_m[n] = 0 \quad \text{for }n<0 \text{ or }n \ge B+L-1 $$

$y_m[n]$ could be non-zero for $0 \le n < B+L-1$.

unlike the overlapping input segments $x_m[n]$, the length of the output grains $y_m[n]$ need not be the block length $B$, but can be (and are) longer.  so instead of 50% overlap (which is $\frac{B-\tfrac12 B}{B} = \frac12$ or an overlap of $B - \tfrac12 B$ samples) like you might have with the Hann-windowed input, your overlap with the output blocks will be $B+L-1 - \tfrac12 B$ samples.

FFT convolution

the convolution of each block

$$ y_m[n] = \sum\limits_{i=0}^{L-1} h[i] x_m[n-i] $$

can be done with an FFT as long as the FFT size $N$ is at least the length of the output block.

$$ N \ge B+L-1 $$

in this case

$$\begin{align}
X_m[k] & = \mathcal{DFT} \left\{ x_m[n] \right\} \\
       & = \sum\limits_{n=0}^{N-1} x_m[n] \ e^{-j 2 \pi \frac{nk}{N}}  \\
\end{align}$$

$$ Y_m[k] = H[k] \cdot X_m[k] \quad \text{for } 0 \le k < N $$

where

$$\begin{align}
H[k] & = \mathcal{DFT} \left\{ h[n] \right\} \\
       & = \sum\limits_{n=0}^{N-1} h[n] \ e^{-j 2 \pi \frac{nk}{N}}  \\
\end{align}$$

$$\begin{align}
Y_m[k] & = \mathcal{DFT} \left\{ y_m[n] \right\} \\
       & = \sum\limits_{n=0}^{N-1} y_m[n] \ e^{-j 2 \pi \frac{nk}{N}}  \\
\end{align}$$

you need only compute $H[k]$ once for use with all of the blocks. and you get the time-domain samples for each output block with the inverse FFT.

$$ \begin{align}
y_m[n] & = \mathcal{IDFT} \left\{ Y_m[k] \right\} \\
       & = \frac1N \sum\limits_{k=0}^{N-1} Y_m[k] \ e^{+j 2 \pi \frac{nk}{N}}  \\
\end{align} $$


When is aliasing a good thing? In Hamming's book, The Art of Doing Science and Engineering, he relates the following story:

  A group at Naval Postgraduate School was modulating a very high
  frequency signal down to where they could afford to sample, according
  to the sampling theorem as they understood it.  But I realized if they
  cleverly sampled the high frequency then the sampling act itself would
  modulate (alias) it down.  After some days of argument, they removed
  the rack of frequency lowering equipment, and the rest of the
  equipment ran better!

Are there any other ways to use aliasing as a primary technique for processing a signal, as opposed to a side-effect to be avoided?

The quoted text in the question is a case of using bandpass sampling or undersampling.

Here, to avoid aliasing distortion, the signal of interest must be bandpass. That means that the signal's power spectrum is only non-zero between $f_L < |f| < f_H$.

If we sample the signal at a rate $f_s$, then the condition that the subsequent repeated spectra do not overlap means we can avoid aliasing. The repeated spectra happen at every integer multiple of $f_s$.

Mathematically, we can write this condition for avoiding aliasing distortion as

$$\frac{2 f_H}{n} \le f_s \le \frac{2 f_L}{n - 1}$$

where $n$ is an integer that satisfies

$$1 \le n \le \frac{f_H}{f_H - f_L}$$

There are a number of valid frequency ranges you can do this with, as illustrated by the diagram below (taken from the wikipedia link above).

In the above diagram, if the problem lies in the grey areas, then we can avoid aliasing distortion with bandpass sampling --- even though the sampled signal is aliased, we have not distorted the shape of the signal's spectrum.


What's the model for a DAC followed by a sampler? I have a situation where a sampled signal is passed through a DAC and re-sampled by a switched capacitor circuit.

I can't see a model for this in the literature.

The standard model of a DAC is the ZOH. But when it gets re-sampled by the switched cap circuit it seems to me that ZOH effect of the DAC disappears because of the re-sampling.

What does the Signal Processing SE community think?

Assume the re-sampling is at the same rate as the DAC.

What kind of DAC is that? Standard instrumentation DAC (modeled by a ZOH, with a negligible settling time and no reconstruction filter) or something fancier?

What's happening could be described in the following terms: in the short analog land between the DAC and the rest of the circuit, you see your signal with mirror images of its spectrum above Nyquist frequency (amplitude decreasing with frequency, as your ZOH is equivalent to a multiplication of the spectrum by a sinc function). The switched capacitor circuit resamples this signal; but since its sample rate is the same as that of the original signal, there is no aliasing - the mirror images created by the DAC/ZOH will line up perfectly with the original harmonics.