How do I identify repeated subimages within a larger image? I have a large image of a microscopic die of an integrated circuit. The pixel dimensions are something like 15000 x 10000.

Example image:
Commodore Amiga Paula Image http://tinytransistors.net/images/chips/Amiga_A500/Paula_total.jpg

The image contains a number of repeated structured subimages. Some large blocks might even repeat 3 or 4 times because there was a need for duplicated functionality. (See good candidates on the left of the example image)

I'd like to identify the largest(in pixel dimensions) repeating subimage block, and then highlight the locations of those subimages. I would repeat this process looking for the next largest repeating subimage block, highlight those, and so on.

To be clear, I don't have the subimages in advance. The problem isn't only to find the location of certain subimages within an image. That's Step two. :)

While the image is in fact a color image, it could be easily turned into a gray-scale image without loss of content.

There is obviously some "noise" to an image like this, and there'd have to be some adjustable "fudge factor" on what to consider a match.

I'm not entirely sure where to start on this, or if this is even the correct SE. Would FFT be a potential solution? Using MatLab? Specific one-off image analysis tools?

I could write some simple image analysis software that breaks an image up into smaller blocks, and then searches the image for those blocks. But this sounds incredibly slow to me, and there's got to be a more efficient way to do this.

A naive approach would be to set a maximum size, .25 image width, and compare every possible square with every other possible square, in a "sliding window" fashion.

So take a 5000 by 5000 square as your reference (eg 0,0 to 4999,4999).
now grab all the other possible "windows". Start by comparing your reference square to the square at (5000,0 to 9999,4999), then compare your reference square to the square at (5001,0 to 10000,4999) then to (5002,0 to 10001,4999) etc.

So your compare function could subtract the (greyscale) difference between each pixel of your reference frame and your target frame and give a score.

Next (after comparing your reference to all other possible squares), you select a new reference, the square at (1,0 to 5000,4999) and compare that reference square to all other possible squares.

Once the whole process was done at the 5000,5000 size, you would know if any structure of size 5000x5000 matches any other part of the image (to a threshold).

Then, you would make your reference square smaller, eg 4999x4999  and repeat the whole process at size 4999x4999 to see if any structure of size 4999x4999 is "duplicated"

You would proceed until your minimum size.

Then...you would want to try instead of squares...rectangles of different proportions.

As you can see, this is computationally prohibitive, so you have to make compromises:

make the source image smaller. Only when you have a near match, check the high res map at that location
slide by 10 pixel increments
identify some basic reference structures yourself, and only use those as reference
try a "monte-carlo" approach by trying random locations and getting statistics of promising locations for pixel by pixel checks
implement on gpu
limit the structures to only a few rectangle shapes rather than all
preprocess the image to have higher contrast
start by using small sizes (20x20 pixels) and once you find small matches, grow your "square" there to look for larger matches that might be the "real" duplicated structure

So there's no magic. What there is is domain heuristics.  You tell your script what you think are the likely sizes/shapes/locations/orientations of structures.

However, the question of "why" you need to identify duplicate structures is actually fundamental to how you will optimize the problem.

EDIT:
an entirely different approach would be to use some algorithm that's pre-tuned for "object segmentation" to go through and find anything that seems to have a defined boundary, then run a clustering algorithm on all the "objects", and it should group the most similar ones together. Then you can manually select the matches or do a final threshold step. This is likely much more similar to how the "human" brain does it as is likely to be much faster


How to perform sine generation (in C)? I am interested in how should a listing for a program look like that will generate a sine of f=50hz, using the next arguments for any moment in time.

float phis=162*PI/180;//initial phase of signal
////////////////////////////////////////////////////////////////////////
//  generates a sinusoid of dwLength in dwData
//      and maintains phis value updated
////////////////////////////////////////////////////////////////////////
void rc(short* dwDataOut, DWORD dwLength){
float omegaf=50/(float)dwLength/PI*90/100;
DWORD dwi;
//use phis
    for (dwi=0;dwi<dwLength;dwi++)
        dwDataOut[dwi]=4450*sin(dwi*omegaf+phis)+200;

//update wave phase
    phis-=PI*(1-50/(float)dwLength)+PI/2;
    while(phis<0)
        phis+=2*PI;
}

The program I came up works correctly with a small marj of error only for f=50hz, because right from the start it appears with a small error in getting phis coefficient.

By applying this procedure in a program I'll make use of a cleaner signal from a sound wave, and obtain a signal that will be overlapped as antiphase signal to a more complex 50hz wave, for longer period of time, and so eliminating this harmonic. I’m interested in a more mathematical method to find the value for phis, and probably the arguments of sine function.

(is this function resulting in a true output for f=50hz? what are the errors that can be adjusted?)

I think your recalculation of phis in the function is not correct.

If I do this in scilab:

//25631
N = 100;
phis1=162*%pi/180;//initial phase of signal

[data1,phis2] = rc(N,phis1);
[data2,phis3] = rc(N,phis2);
[data3,phis4] = rc(N,phis3);

clf
plot([data1 data2 data3])

where rc is

////////////////////////////////////////////////////////////////////////
//  generates a sinusoid of dwLength in dwData
//      and maintains phis value updated
////////////////////////////////////////////////////////////////////////
function [dwDataOut,phis] = rc(dwLength, phis)
omegaf=50/dwLength/%pi*90/100;

dwDataOut = zeros(1,dwLength);

for dwi=1:dwLength,
    dwDataOut(dwi) = 4450*sin(dwi*omegaf+phis)+200;
end

//update wave phase
//WAS: phis = phis - %pi*(1-50/dwLength)+%pi/2;
phis = dwLength*omegaf+phis;
while(phis<0)
    phis = phis + 2*%pi;
end
endfunction

Note specifically the line phis = dwLength*omegaf+phis;

Then the plot is as below, which seems to be what you're after.


differences between DIT & DIF algorithms What are the differences between decimation in time and decimation in frequency algorithms of FFT, especially as their names suggest?

How can I see/understand that decimation in time domain is taking place in DIT and decimation in frequency domain is taking place in DIF?

The DFT is (let's say that $N$ is even):

$$ X[k] \ \triangleq  \sum\limits_{n=0}^{N-1} x[n] \, e^{-j 2 \pi \frac{nk}{N}} $$

The radix-2 FFT works by splitting a size-$N$ DFT into two size-$\frac{N}{2}$ DFTs.  (Because the cost of a naive DFT is proportional to $N^2$, cutting the problem in half will cut this cost, maybe, in half.  Two size-$\frac{N}{2}$ DFTs appear to cost less than one size-$N$ DFT.

The Decimation-in-Time FFT splits the two DFTs into even and odd-indexed input samples:

$$\begin{align}
 X[k] \ &\triangleq  \sum\limits_{n=0}^{N-1} x[n] \, e^{-j 2 \pi \frac{nk}{N}} \\
\\
&=  \sum\limits_{n\text{ even}} x[n] \, e^{-j 2 \pi \frac{nk}{N}} + \sum\limits_{n\text{ odd}} x[n] \, e^{-j 2 \pi \frac{nk}{N}} \\
\\
&=  \sum\limits_{m=0}^{\frac{N}{2}-1} x[2m] \, e^{-j 2 \pi \frac{(2m)k}{N}} + \sum\limits_{m=0}^{\frac{N}{2}-1} x[2m+1] \, e^{-j 2 \pi \frac{(2m+1)k}{N}} \\
\\
&=  \sum\limits_{m=0}^{\frac{N}{2}-1} x[2m] \, e^{-j 2 \pi \frac{mk}{N/2}} + e^{-j 2 \pi \frac{k}{N}} \sum\limits_{m=0}^{\frac{N}{2}-1} x[2m+1] \, e^{-j 2 \pi \frac{mk}{N/2}} \\
\\
\end{align}$$

This shows that if the evens are stuffed into one DFT and the odds are stuffed into another DFT, that the results of the two DFTs can be combined by multiplying the result of the right-hand (odd) DFT with a single factor $ e^{-j 2 \pi \frac{k}{N}} $ and added to the left-hand (even) DFT.  The result for the two frequencies $X[k]$ and $X[k+\tfrac{N}{2}]$ are the same except the right hand term is subtracted than added.

$$\begin{align}
X[k] &= \sum\limits_{m=0}^{\frac{N}{2}-1} x[2m] \, e^{-j 2 \pi \frac{mk}{N/2}} + e^{-j 2 \pi \frac{k}{N}} \sum\limits_{m=0}^{\frac{N}{2}-1} x[2m+1] \, e^{-j 2 \pi \frac{mk}{N/2}} \\
\\
X[k+\tfrac{N}{2}] &= \sum\limits_{m=0}^{\frac{N}{2}-1} x[2m] \, e^{-j 2 \pi \frac{mk}{N/2}} - e^{-j 2 \pi \frac{k}{N}} \sum\limits_{m=0}^{\frac{N}{2}-1} x[2m+1] \, e^{-j 2 \pi \frac{mk}{N/2}} \\
\end{align}$$

This shows how the butterflies are defined on the last pass of the FFT.

The Decimation-in-Frequency FFT splits the two DFTs into the first half and last half of the input samples:

$$\begin{align}
 X[k] \ &\triangleq  \sum\limits_{n=0}^{N-1} x[n] \, e^{-j 2 \pi \frac{nk}{N}} \\
\\
&=  \sum\limits_{n=0}^{\frac{N}{2}-1} x[n] \, e^{-j 2 \pi \frac{nk}{N}} + \sum\limits_{n=\frac{N}{2}}^{N-1} x[n] \, e^{-j 2 \pi \frac{nk}{N}} \\
\\
&=  \sum\limits_{n=0}^{\frac{N}{2}-1} x[n] \, e^{-j 2 \pi \frac{nk}{N}} + \sum\limits_{n=0}^{\frac{N}{2}-1} x[n+\tfrac{N}{2}] \, e^{-j 2 \pi \frac{(n+N/2)k}{N}} \\
\\
&=  \sum\limits_{n=0}^{\frac{N}{2}-1} x[n] \, e^{-j 2 \pi \frac{nk}{N}} +  e^{-j 2 \pi \frac{(N/2)k}{N}} \sum\limits_{n=0}^{\frac{N}{2}-1} x[n+\tfrac{N}{2}] \, e^{-j 2 \pi \frac{nk}{N}} \\
\\
&=  \sum\limits_{n=0}^{\frac{N}{2}-1} x[n] \, e^{-j 2 \pi \frac{nk}{N}} +  e^{-j \pi k} \sum\limits_{n=0}^{\frac{N}{2}-1} x[n+\tfrac{N}{2}] \, e^{-j 2 \pi \frac{nk}{N}} \\
\\
&=  \sum\limits_{n=0}^{\frac{N}{2}-1} x[n] \, e^{-j 2 \pi \frac{nk}{N}} +  (-1)^k \sum\limits_{n=0}^{\frac{N}{2}-1} x[n+\tfrac{N}{2}] \, e^{-j 2 \pi \frac{nk}{N}} \\
\\
&=  \sum\limits_{n=0}^{\frac{N}{2}-1} \Big( x[n]  +  (-1)^k x[n+\tfrac{N}{2}] \Big) \, e^{-j 2 \pi \frac{nk}{N}} \\
\end{align}$$

For even $k$ we add the $x[n+\tfrac{N}{2}]$ to $x[n]$.  For odd $k$, we subtract.

$$\begin{align}
X[2m] &= \sum\limits_{n=0}^{\frac{N}{2}-1} \Big( x[n] + x[n+\tfrac{N}{2}] \Big) \, e^{-j 2 \pi \frac{n(2m)}{N}} \\
\\
 &= \sum\limits_{n=0}^{\frac{N}{2}-1} \Big( x[n] + x[n+\tfrac{N}{2}] \Big) \, e^{-j 2 \pi \frac{nm}{N/2}} \\
\\
X[2m+1] &= \sum\limits_{n=0}^{\frac{N}{2}-1} \Big( x[n] - x[n+\tfrac{N}{2}] \Big) \, e^{-j 2 \pi \frac{n(2m+1)}{N}} \\
\\
 &= \sum\limits_{n=0}^{\frac{N}{2}-1} \Big( x[n] - x[n+\tfrac{N}{2}] \Big) \, e^{-j 2 \pi \frac{n}{N}} e^{-j 2 \pi \frac{nm}{N/2}} \\
\end{align}$$

This shows how the butterflies are defined on the first pass of the FFT.


Difference between PSD estimate and variance of DFT In Bartlett's PSD estimate one averages over L segments of the squared DFT coefficients. From wikipedia I found this formula

$$ \textrm{PSD}(k) = \frac{1}{L}\sum_{l=1}^{L} \frac{1}{M} \lvert X^{[l]}(k)\rvert^2 \tag{1}$$
where $X^{[l]}$ denotes the DFT of the $l^{\rm th}$ segment.

With my (not so standard) DFT Definition

$$
X(k) = \frac{1}{\sqrt{N}} \sum_{n=1}^{N}x[n]e^{-j\omega_kn},  \textrm{with } \omega_k = \frac{2\pi k}{N} \tag{2}
$$

Bartlett's method results in

$$\textrm{PSD}(k) = \frac{1}{L}\sum_{l=1}^{L} \frac{1}{M} \lvert\sqrt{M}\cdot X^{[l]}(k)\rvert^2 = \frac{1}{L}\sum_{l=1}^{L} \lvert X^{[l]}(k)\rvert^2 \tag{3}$$

A similar formula is given in my textbook for the (unbiased) sample variance of the DFT coefficients of a periodic signal:

$$ \hat{\sigma}_x^2(k) = \frac{1}{L-1}\sum_{l=1}^{L} \lvert X^{[l]}(k) - \hat{X}(k)\rvert^2 \tag4$$

with sample mean $\hat{X}(k) = \frac{1}{L} \sum_{l=1}^{L} X^{[l]}(k)$

Assuming that $\hat{X}=0$, equations $(3)$ and $(4)$  seem to be biased & unbiased estimates of the variance $\sigma_X^2(k)$.

Questions:

Always assuming  $\hat{X}=0$, is it generally correct  to interpret the PSD as an biased estimate of the DFT variance? If not, can someone  explain  the difference, please?
Would Bartlett's method become an unbiased estimator for the DFT variance if $(3)$  was scaled by $\frac{1}{L-1}$ instead of $\frac{1}{L}$?

So, as usual, I start by pointing out what the PSD actually is:

The power spectral density is the Fourier Transform of the autocorrelation function of a wide-sense stationary random process.

  Always assuming  $\hat{X}=0$, is it generally correct  to interpret the PSD as an biased estimate of the DFT variance? If not, can someone  explain  the difference, please?

The Wiener-Chintschin¹ theorem points out that this happens to coincide with the expectation value of the absolute of the Fourier transform of a signal.  In other words, to answer your first question:

The PSD is a stochastic (hence, "invisible") property of a random process. The DFT can be used to estimate that; not the other way around!

To answer your second question, a bit of juggling the term of variance is required:

  Would Bartlett's method become an unbiased estimator for the DFT variance if (3) was scaled by $\frac{1}{L-1}$ instead of $\frac{1}{L}$?

Variance is formally defined as the expectation of the square of the deviation of any observation to the mean:

$$Var(X) = \mathbf E\left[\left(X-\mu\right)^2\right]$$

Now, applying Parseval's theorem, the power in a time and frequency domain should be identical. Considering a white random process, where the covariance of two samples is 0, it's easy to see that variance is kind-of-linear:

$$Var\left(\sum\limits_{i=1}^L X_i\right) = \sum\limits_{i=1}^N Var(X_i)$$

Considering the special case of L=1, it becomes clear that $(4)$ isn't applicable here – you kind of get the standard deviation of the process, but not the variance.

¹ there's more latin-lettered spellings to Chintschin than there's letters in in his cyrillic name...


Deep Learning: Classification vs. Convolution for Signal Restoration Assume we have vector $X = [x_1, x_2, x_3, x_4 ,..... x_N], ∈ -1,1$. Therefore the value of $x$ is either 1 or -1.  The vector $X$ is convoluted with random generated vector $Y$ whose length is the same of $X$, so $r = X*Y$ where $*$ denote to the convolution operation.

Is there classification deep learning model which can extract the vector $X$ based on $r$ ?? I mean the input should be the vector $r$ and output the vector $X$.

thank you in advance.

By deep learning, I'll assume you mean neural network. To develop a neural network, you'll need labeled data. This means you need a bunch of example inputs ($r$) and outputs ($X$). Once you have that, you need to build a network. As far as I can tell, currently the guidelines for building these things are very loose and ad-hoc. There is general consensus on good or bad activation functions and general suggestions, but there is no formal theory about how many neurons to choose, how many layers, etc. (at least since the last time I took a dive into this).

Common problems will arise such as overfitting and there are techniques to help this but you need to be smart enough to spot the signs and then apply the fixes. For example, you can read about how observing your training and testing accuracy can help determine if a network is overfitting.

The better part of this answer is this: what you're describing is called deconvolution (https://en.wikipedia.org/wiki/Blind_deconvolution) and has been accomplished using "traditional" signal processing techniques long before the deep learning craze. Personally, any time a problem has a structured solution, I go in that direction. Reliable neural networks can be hard to develop despite all the helpful libraries in the world and when it comes to troubleshooting I think you'll find more traditional "tried-and-true" methods to be much easier to not only get working but also to understand why they work and this will lend itself much better to be able to adapt to your specific problem.


Why is telephone audio sampled at 8 kHz? When did we decide to sample telephone at $8$ kHz? Has this always been the case? Why did we do that? Is it because higher bit rates can't be transferred as quick? And do these reasons still count? And if not, why isn't there a new standard yet? Is it true that $8$ kHz is the lowest possible sampling rate to transfer understandable speech?

I'm trying to find sources for this, but there doesn't seem to be much information about it.

If anyone cares to dig I think they will find that before Bell Telephone started multiplexing voice lines they did a lot of research on the frequency content of the human voice.  They originally used test groups to develop the bel audio unit and the power distribution of the voice as well as the sensitivity of the human ear to various frequencies.  They developed a band pass characteristic that peaked around 2.1 KHz and rolled off below 300 and over 3000 HZ.  That gave a good human sounding voice when done correctly.  All that was analog.
AM radio expanded that to 5 KHz to include music that was acceptable to most folks when we were young and had good ears.  Television flyback transformers were designed to run at ~17.5 KHz because there was a magic number for picture reproduction and most folks could not hear the whine.  Single sideband radio was commercialized in the 1960's and needed very sharp cutoff frequencies.  I used radios with filters at 2.1 and 3.1 KHz. 2.1 had some Donald Duck characteristics.  3.1 sounded fine, again with young ears.  Audio bandpass was increased to 20 KHZ or better with FM because the higher carrier frequencies could handle higher bandwidth for better music reproduction.  Stack up some xylophones or bells or other high pitched instruments and they can get enough harmonic energy into the higher frequencies.  OTOH, as was state, most folks cannot hear it.

The bottom line is that anyone claiming they need 20 KHz bandwidth for voice is not paying attention.  3 KHz will do it, 5 will give you some margin.  If it does not sound right then something other than the bandwidth is the issue.

When digital signaling was being developed folks who know figured out that no matter how weird a waveform looked it could be broken down into a set of sine waves.  The harmonic mixing of those waves produced the typical spikey pattern of voice or music.  Lastly, Nyquist did research on the digital sampling rate needed to reproduce a sine wave at a given frequency.  It turns out that it takes 2 samples to make a sine wave so the highest frequency that will be reproduced is half the sample rate.  You want 5 KHz of audio then sample at 10 KHz.  Fine for voice.  You want higher fidelity music than most folks can hear then sample at 40 KHz or so to get 20+ KHz.

One more tidbit is sampling vs bitrate.  If you sample at a given frequency then multiply that by the word length you will get the minimum bitrate needed to produce the desired signal.  Reduce bitrate and the size of the sample word will get cut to meet the new bit rate at a given sample rate.  That is all "lossless" encoding. This is all from memory and trying to find current data.  It's there if somebody looks for cites.  I'm not going to bother as I am getting too old to care.  I just got tired of wading through a lot of obvious mythical issues when I got interested in doing some audio capture.


Adding measuruments noises to kalman filter I'm implementing navigation system for my robot. There are two ways to get data from it: odometry(encoders from motors) and camera. Both of information sources can give me estimate of robot's position. But camera has very big noise. I've found that it is useful to use Kalman filter. After some literature looking I have a questions:

How to use kalman filter with adding measurement noise $v(k)$?
$$x(k+1) = F(k)x(k-1)+B(k)u(k)+w(k)$$
$$   z(k) = H(k)x(k)+v(k)$$
How to make "coefficient of trust" for camera measurements in Kalman gain very low (because     of very very high camera noise)? (Probably is it needed to make specific Q or R matrix?)

I would appreciate any help. Thanks.

For the first question, in Kalman filter there are two noises namely process and observation noises. Kalman filter assumes the noise to be Gaussian with zero mean and some variance $\mathcal{N}(0, \sigma^{2})$. Let's say we have a sensor that measures a constant quantity but random, for example distance. Let's say $x = 1$m, however the sensor gives us $z=0.9$m. In reality we don't know $x$ for certain so we rely on our sensor. In Matlab, we can simulate this sensor

x = 1; % true value
z = x + randn(); % observed value

randn() is normal distribution with zero mean and variance one $\mathcal{N}(0, 1)$. If you want to increase the accuracy of your sensor, then multiply the normal distribution randn() with some small value that you think it is the error of your sensor.

for i = 1:6
    x = 1; % true value
    z = x + 0.01*randn() % observed value
end

You will notice the observed value varies from the true value however still close to the true value because of 0.01. Remember that Kalman filter is a matter of modeling your system accurately. For kalman filter, you need the motion and observation models given in advanced. Also, your system must be linear and if this is not the case then you should use Extended kalman filter. Moreover, the noise must be Gaussian noise. Now your job to determine these parameters.


Pre-emphasis: What is connection between parameter a and 6db/oct? I know that the standard pre-emphasis function is $y[n] = x[n]-\alpha\cdot x[n-1]$ and $\alpha$ is often set as $0.97$.

However, when I was learning speech processing using Praat, the pre-emphasis factor in it is $6\ \rm dB/oct$ which means the signal emphasise $6\ \rm dB$ per octave increase. That's what confused me.

I've read lots of questions and answers, most of them are explaining how to calculate the value of $\alpha$ or the $\mathcal Z$-transformation. But What I'm confused is the connection between $\alpha$ and $6\ \rm dB/oct$.

I found this equation in the notes of Praat:

  The pre-emphasis factor $\alpha$ is computed as:
  $$
\alpha = \exp(-2\pi F \Delta t)
$$
  where the $F$ above which the spectral slope will increase by $6\ \rm dB/oct$, $\Delta t$ is the sampling period of the sound

According to above notes, it seems that the $\alpha$ only influences the "beginning frequency" of pre-emphasis but not "how much it changes".

Hope I explained myself clearly enough.

You're right that the value $\alpha$ will not influence the slope of the filter (approx. $6$ dB/octave). That slope is determined by the filter order, which is $1$ in this case (1 delay element in the filter implementation). The value $\alpha$ only influences the filter's cut-off frequency.


What is the unit for power spectrum of acceleration signals in log scale I have a signal recorded by an accelerometer in m/s^2. It is basically accelerations over a period of time. I have calculated the power spectrum of the signal over in a certain frequency band and then transferred into a 20log(x) scale for plotting. I am wondering what should be unit for it?

Since it has been transferred into a log scale, I could use dB, but the standard formula for dB calculation is something like 20log(x/x0), where x0 is the reference and cancels out the unit leaving us with just dB. I am not sure if I can the same thing in the above case since I don't have any reference power? Can I just assume something like x0=1 and report the numbers in dB?

If not, then I am not sure what should accompany dB. In one paper, I have seen for accelerations reported in gravities 1g = 9.8m/s^2,dB g^2/2 was used as the unit for the power spectrum? This makes me guess maybe I should use dB m/2s^2 or weirdly dB m^2/2s^4? Any ideas are greatly appreciated.

The typical practice with decibels is to append the unit.  Unfortunately, it's informal, so you'll see things like "dBm" in RF circuits which refers to "dB with a reference of 1mW", not "dB with a reference of 1 meter".

you could:

Put a note in your graph that $0\mathrm{dB} = \mathrm{1 m/s^2}$
Use "dBg" (and use 1g as your reference)
Just keep the y-axis as the log of acceleration squared (and yes, $\mathrm{m^2/s^4}$ is awkward, but it's the correct unit).


Effect of Windowing Function on Autocorrelation Function I need to compute FFT components and autocorrelation of a discrete signal. When I do not use any window function before taking FFT of the signal I get a reasonable AC function from FFT that includes some information about period of the signal:

However when I apply a window function (Hanning) before applying FFT I got an AC function that seems to be distorted:

What is the effect of window function on an AC function that is computed from FFT?

When you are using the FFT to get a measurement of the frequency content of a signal it is appropriate to window the data first to avoid distortions due to a rectangular cutoff.

That is not what you are doing with the FFT, though.  You are using the FFT and inverse FFT to efficiently calculate a correlation (an autocorrelation is just a special case of correlation).  Done properly with padding and such, the FFT will get you the exact same results that you would get by calculating the correlation in the time domain.  If you window, though, you will not get the same results.  In short, you should not window when you are using the FFT to calculate correlations or convolutions.


Choice of cost function in adaptive noise cancellation Why is the least squares cost function used in adaptive noise cancellation with the recursive least squares (RLS) algorithm, but the mean squared error used in signal estimation?

For an ergodic source (where time and ensemble averages should be the same) I would have thought the choice was arbitrary.

From my notes: LMS minimises $\|e[n]\|^{2}$ and RLS minimises $\sum_{i=0}^{N} e^{2}[i], $
where $e[n]$ is the difference between the desired signal and the output of the filter.

You are actually asking 2 separate questions. If I have interpreted your post correctly, they are:

Why is the RLS used in Noise Cancellation, whereas the LMS is used in signal estimation?
(If we assume that the signal is ergodic) Isn't minimising the Ensemble average, $E\{e_n^2 \}$, equivalent to minimising the time average, $ \sum_{i = 1}^{N}e_i^2$; i.e. Shouldn't the LMS perform as good as the RLS?

Answers:

Q1:

You are free to choose ANY adaptive filtering algorithm for ANY application (e.g. Adaptive Noise Cancellation -ANC).  It need not be the RLS. The same goes with signal estimation, or any other application you can think of.

However, if you have a specific application in mind, and you know something about the system or how the algorithm performs in that setting, you may want to choose one adaptive algorithm over another. (E.g. If your application is to do with converging faster as opposed to having a low error at steady state, you might choose Algorithm X over Y)

Once again, this choice is not made based on whether you minimize the mean squared error or the sum of the squared error.

Q2:

This questions still keeps me up on some nights. Fundamentally, Minimum Mean Squared  Error
based algorithms are based on, finding the system parameters which minimise the $ E\{ e_n^2 \}$. Now, the expected value of the error

Both the LMS and the RLS are based on minimising this quantity. How they do it, however, is quite different. Remember, that the Expected value is statistical operator - it's the average over all realisations of the random signal. In real life applications, we only have access to one realisation of this signal.

Example: Let's say you have the price of a stock, all you get is one value at each time (e.g. £10) The stock price never says, "oh well I have a 5% chance of being £10, on the other hand I could be £12, or £14.....etc". But I digress.

SO,
How are we going to minimize the mean squared error, when we can't compute the true mean?
The Least Mean Square (Widrow-Hoff) says, "You know what, I am going to drop the $E\{.\}$ and minimise $e_n^2$. Yes it is a hack, but my god what a hack!

The RLS uses the fact you mentioned: Since the process is ergodic, you can get the time average instead of the ensemble. The sum squared error is equivalent to the average of the squared error. The only thing that's missing is the $1/N$ term in the sum squared error. It is discarded because, it cancels out with another $1/N$ term in the derivation of the RLS. So don't worry.

The interesting thing is how this difference affects the performance of the 2 filters. (This is a whole different question and I suggest that you ask it here - it might get a lot of people interested) The RLS is known to converge faster than the LMS. But when tracking time-varying parameters the LMS can perform better than the RLS in certain conditions.

LMS:
$$ E\{e_n^2 \} \approx e_n^2$$
RLS:
$$ E\{e_n^2 \} \approx \frac{1}{N} \sum_{i =1}^N e_i^2$$

To sum up: Both the RLS and LMS try to approximate $E\{e^2\}$. They way they do it is different - so the way they behave is different. And you are free to choose any filter for any application.


Timing recovery lock in receiver I have done my receiver design with QPSK modulation. I am trying to communication between two systems at L band. My data rate is 1 Mbps. As my data rate is high I have allocated 64 bit for timing recovery.

Timing recovery has farrow interpolation filter followed by Gardner timing recovery, loop filter (2nd order) and interpolation control for the farrow interpolation filter. I am doing 2 samples/symbol analysis.

What I am assuming is the fractional value to the interpolation filter will stabilize once 64 bit allocated for timing recovery is over. But what I observe is the fractional value varies even after the 64 bit allocated is over because of which my data sometimes is wrong.

can the fractional input to the interpolation filter vary?
if yes why is the data sometime is wrong even if it is just direct board
to board communication?

As far as the recovery time and results for the timing recovery, this is a loop implementation and you would need to review the complete loop for stability and gain parameters to balance loop bandwidth & convergence time, stability and noise performance.

Even with direct board to board communication you will have a static timing offset as the receiver needs to determine the optimum point in which to sample the recovered data in each symbol.

Step response testing of the timing error at your loop filter input and output can be insightful; the input to the loop filter will likely be quite noisy but you can filter your observation (as long as the filter bandwidth is much wider than your expected loop bandwidth this will be valid) to be able to see how it responds to a step change in time offset. Even better is if you have an accurate simulation model of your receiver so that you can observe the time offset more directly, and to observe what occurs as you change loop gain parameters.

Also through simulation and with testing on your hardware I would recommend creating an eye-diagram to better understand the core problems you may be having with data recovery. With this you can best observe the predicted sampling location for making a data decision versus the ideal.


Downsampling from 100Hz to 60 Hz I have two sensors measuring a single time varying signal . The first sensor  (reference sensor) is sampling the data at 60 Hz and the second sensor (test sensor) is sampling the data at 60 Hz.  I need to check the accuracy of the test sensor by comparing with the reference sensor. I am confused about how to downsample the data recorded at 100 Hz to 60 Hz.  I see it is quite trivial to downsample from 100 Hz to 50Hz/25 Hz. But I am confused how to downsize to 60 Hz.

After downsampling, is cross-correlation a good technique to compare the data from two sensors?

Can you please help?

Thanks,
Peter

Are the samples adhering to Nyquist (properly pre filtered)?
If so, I would just resample to one common rate, prefereable >= max(rate1, rate2). And align them in time. Check out MATLABs resample() function.
If the samples are not adhering to Nyquist, then you need to think about what they represent, how they are sampled and in what respect you need to compare them.


Matlab Hilbert Transform in C++ I asked this at Stack Overflow, they told me to come here.

First, please excuse my ignorance in this field, I'm a programmer by trade but have been stuck in a situation a little beyond my expertise (in math and signals processing).

I have a Matlab script that I need to port to a C++ program (without compiling the matlab code into a DLL). It uses the hilbert() function with one argument. I'm trying to find a way to implement the same thing in C++ (i.e. have a function that also takes only one argument, and returns the same values).

I have read up on ways of using FFT and IFFT to build it, but can't seem to get anything as simple as the Matlab version. The main thing is that I need it to work on a 128*2000 matrix, and nothing I've found in my search has showed me how to do that.

I would be OK with either a complex value returned, or just the absolute value. The simpler it is to integrate into the code, the better.

Thank you.

If your application can use LGPL open-source code, you could reuse the code for the Remez exchange algorithm in the GNU Octave source.  I've got a standalone library that does that at www.ke5fx.com/dsplib.zip .  See if the function DSP_FIR_PM_Hilbert_transformer_design() does what you want.

The FIR filters implemented by my library are one-dimensional affairs, so if you actually need to do a 2D transform, that would be left as an exercise, so to speak.  Contributions back to dsplib are welcome!


How to reconstruct a sound from magnitude spectrogram? I have an audio magnitude spectrogram but I don't have the phase, try to randomize the phases of each container and then make a reverse fourier, but only pure noise is heard

How can I reconstruct the phase information so that the audio is heard somewhat decent.

Anyone know an algorithm that does this?

In general the whole complex Fourier transform is needed to reconstruct the time domain signal.

However, it's possible to reconstruct a real signal it from its Fourier transform magnitude alone (or phase alone) by using some iterative techniques, or very large matrix solutions.

Have a look at the papers by Monson Hayes...


How to enhance Pedestrian Dead Reckoning using Particle filter? I am working on Pedestrian dead Reckoning using Particle filtering. I am using mobile accelerometer sensor to detect number of steps and stride length and calculating the total distance traveled.

I want to use particle filter to enhance accuracy in my results. I have read the theory of particle filter but stuck at the point that how to apply the model to my case.

In particular my problem, I have the following questions:

How to determine the initial position pdf distribution? I have a Bluetooth beacon to detect the starting position. Particle filter requires a cluster of point at the start, right? I don't think in my case the distribution is Gaussian. Then what should it be?
Should I do the particle filtering on the smartphone or pass the data back and do it on the computer using Matlab? Will the computation very intense.

How to determine the initial position pdf distribution? I have a Bluetooth beacon to detect the starting position. Particle filter requires a cluster of point at the start, right?

No, the distribution of a particle filter does not need to be a cluster when it starts.  A classic case is for a robot to figure out where it is.  At first it does not know where it is.  It takes sensor readings wherever it happens to be at, but that is usually not enough to localize it (if all you needed was a single set of readings, you wouldn't need a particle filter).  Instead, it moves around and takes more readings.  All of the many spread out particles represent possible positions.  As it moves and takes more readings more and more of the possible positions become highly unlikely, until the distribution collapses to a single cluster.  At that point the robot "knows" (or at least thinks it knows) where it's at.

  I don't think in my case the distribution is Gaussian. Then what should it be?

Without giving us more information there is no way to answer that question.  Since it sounds like your initial estimate comes solely from your bluetooth beacon, you will probably need to generate an empirical data-based model of its accuracy.

  Should I do the particle filtering on the smartphone or pass the data back and do it on the computer using Matlab? Will the computation very intense.

Either way.  The computational load will largely depend on how many particles you use.  There is no right answer as to how many particles you need.  More particles means more computation but also lower chance of getting a wrong answer.


Carrier frequency offset Can i track any frequency offset my increasing samples/symbol.
For example if my symbol rate is 1Msps. If my sample rate is 10 Msps.
Can my receiver correct upto 5 Msps carrier offset using FFT
I can follow it by costas carrier recovery to further correct it before timing?

Ok so there are a couple ways to deal with this. As you pointed out there are things like the Costas Loop and PLL but I think band edge recovery is more intuitive. You can also continually apply STFT and track the frequency offset though you have the problem of time vs frequency localization.

I'm not an expert but I would probably square the signal and find the FFT peak, afterward I would do timing recovery and then apply decision directed feedback or a PLL after demodulation to track the remaining offset.

The band-edge filter timing recovery is pretty intuitive though fred harris derivation ignores some important things/obscures them for some reason. EX he tries to argue maximum-liklihood as a way to derive it but his derivation ignores issues like symbols not being time limited so it hides lots of stuff,


Implementation of Block Orthogonal Matching Pursuit (BOMP) Algorithm How would one implement the Block Orthogonal Matching Pursuit (BOMP) Algorithm in MATLAB?

The Block Orthogonal Matching Pursuit (BOMP) Algorithm is basically the Orthogonal Matching Pursuit (OMP) Algorithm with single major difference - Instead of selecting single index which maximizes the correlation we chose a set of indices, sub set of columns of the matrix and the solution vector.

A good reference for the algorithm is given in:

An Optimal Condition for the Block Orthogonal Matching Pursuit Algorithm.
Block Sparsity: Coherence and Efficient Recovery.

The code is given by:

function [ vX ] = SolveLsL0Bomp( mA, vB, numBlocks, paramK, tolVal )
% ----------------------------------------------------------------------------------------------- %
%[ vX ] = SolveLsL0Omp( mA, vB, paramK, tolVal )
% Minimizes Least Squares of Linear System with L0 Constraint Using
% Block Orthogonal Matching Pursuit (OMP) Method.
% \arg \min_{x} {\left\| A x - b \right\|}_{2}^{2} subject to {\left\| x
% \right\|}_{2, 0} \leq K
% Input:
%   - mA                -   Input Matirx.
%                           The model matrix (Fat Matrix). Assumed to be
%                           normlaized. Namely norm(mA(:, ii)) = 1 for any
%                           ii.
%                           Structure: Matrix (m X n).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
%   - vB                -   input Vector.
%                           The model known data.
%                           Structure: Vector (m X 1).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
%   - numBlocks         -   Number of Blocks.
%                           The number of blocks in the problem structure.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: {1, 2, ...}.
%   - paramK            -   Parameter K.
%                           The L0 constraint parameter. Basically the
%                           maximal number of active blocks in the
%                           solution.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: {1, 2, ...}.
%   - tolVal            -   Tolerance Value.
%                           Tolerance value for equality of the Linear
%                           System.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range [0, inf).
% Output:
%   - vX                -   Output Vector.
%                           Structure: Vector (n X 1).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
% References
%   1.  An Optimal Condition for the Block Orthogonal Matching Pursuit
%       Algorithm - https://ieeexplore.ieee.org/document/8404118.
%   2.  Block Sparsity: Coherence and Efficient Recovery - https://ieeexplore.ieee.org/document/4960226.
% Remarks:
%   1.  The algorithm assumes 'mA' is normalized (Each column).
%   2.  The number of columns in matrix 'mA' must be an integer
%       multiplication of the number of blocks.
%   3.  For 'numBlocks = numColumns' (Equivalent of 'numElmBlock = 1') the
%       algorithm becomes the classic OMP.
% Known Issues:
%   1.  A
% TODO:
%   1.  Pre Process 'mA' by normalizing its columns.
% Release Notes:
%   -   1.0.000     19/08/2019
%       *   First realease version.
% ----------------------------------------------------------------------------------------------- %

numRows = size(mA, 1);
numCols = size(mA, 2);

numElmBlock = numCols / numBlocks;
if(round(numElmBlock) ~= numElmBlock)
    error('Number of Blocks Doesn''t Match Size of Arrays');
end

vActiveIdx      = false([numCols, 1]);
vR              = vB;
vX              = zeros([numCols, 1]);
activeBlckIdx   = [];

for ii = 1:paramK

    maxCorr         = 0;

    for jj = 1:numBlocks
        vBlockIdx = (((jj - 1) * numElmBlock) + 1):(jj * numElmBlock);

        currCorr = abs(mA(:, vBlockIdx).' * vR);
        if(currCorr > maxCorr)
            activeBlckIdx = jj;
            maxCorr = currCorr;
        end
    end

    vBlockIdx = (((activeBlckIdx - 1) * numElmBlock) + 1):(activeBlckIdx * numElmBlock);
    vActiveIdx(vBlockIdx) = true();

    vX(vActiveIdx) = mA(:, vActiveIdx) \ vB;
    vR = vB - (mA(:, vActiveIdx) * vX(vActiveIdx));

    resNorm = norm(vR);

    if(resNorm < tolVal)
        break;
    end

end

end

The MATLAB code is available at my StackExchange Signal Processing Q60197 GitHub Repository.
In the full code I compare the Block implementation to OMP to verify the implementation.


List of all available ARM+DSP processors for surveillance I have come across TI's C6000/Integra. I want to know is there any other similar or better ARM+DSPs available. Specially for surveillance/vehicle number plate recognition etc. Having camera input at least for VGA.

Link: TI's-Integra

As you mention other non-ARM options in the comments, here's a few:

Blackfin
Xilinx's Zynq (ARM+FPGA)
Altera also do a similar ARM+FPGA
Some of Freescale's i.MX51 processors might be suitable (ARM+co-processor, not a true DSP) - although I'm not sure if there's out-of-the-box surveillance software
And of course, there's always the PC+GPU option, depending on your power, cost, size, weight constraints


Difference between $\mathbb{E}[\mathbf{x} \mathbf{x}^{\rm{H}}]$ and $\mathbb{E}[(\mathbf{x}-\boldsymbol{\mu}) (\mathbf{x}-\boldsymbol{\mu})^{\rm{H}}]$ Let us have a random vector $\mathbf{x} \sim \mathcal{CN} (\boldsymbol{\mu}, \boldsymbol{\Sigma})$ with $\boldsymbol{\mu} \neq \mathbf{0}$. What can we say about the relationship between the elements of $\mathbf{x}$ in the following two separate cases?

1) If $\boldsymbol{\Sigma} \triangleq \mathbb{E}[(\mathbf{x}-\boldsymbol{\mu}) (\mathbf{x}-\boldsymbol{\mu})^{\rm{H}}]$ is a diagonal matrix;

OR

2) If $\mathbf{R} \triangleq \mathbb{E}[\mathbf{x} \mathbf{x}^{\rm{H}}] = \boldsymbol{\mu} \boldsymbol{\mu}^{\rm{H}} + \boldsymbol{\Sigma}$ is a diagonal matrix.

$ $

In 1), the elements of $\mathbf{x}$ are all independent. But what is the interpretation of 2)?

For the first case, as you wrote, it means the elements are not correlated. Since this is a Gaussian Random Vector it means the elements are independent.
It means that at most only one element of $ \boldsymbol{\mu} $ is not zero.
Since if there were more than 1, the matrix $ \boldsymbol{R} $ wasn't diagonal.

Update

Let's define $ \hat{\boldsymbol{x}} = \boldsymbol{x} - \boldsymbol{\mu} $, namely its centered version.

Looking ta the correlation matrix elements:

$$\begin{align*}
\boldsymbol{R}_{i, j} & = \mathbb{E} \left[ \boldsymbol{x}_{i} \boldsymbol{x}_{j} \right] = \mathbb{E} \left[ \left( \hat{\boldsymbol{x}}_{i} + \boldsymbol{\mu}_{i} \right) \left( \hat{\boldsymbol{x}}_{j} + \boldsymbol{\mu}_{j} \right) \right] \\
& = \mathbb{E} \left[ \hat{\boldsymbol{x}}_{i} \hat{\boldsymbol{x}}_{j} \right] + \boldsymbol{\mu}_{j} \mathbb{E} \left[ \hat{\boldsymbol{x}}_{i} \right] + \boldsymbol{\mu}_{i} \mathbb{E} \left[ \hat{\boldsymbol{x}}_{j} \right] + \boldsymbol{\mu}_{i} \boldsymbol{\mu}_{j} \\
& = \Sigma_{i, j} + \boldsymbol{\mu}_{i} \boldsymbol{\mu}_{j}
\end{align*}$$

Now it is really easy to see what is needed for this element to vanish.


Computing N-point inverse DFT with N/2-point forward DFT I'm trying to implement IDFT using multiple smaller forward DFTs but am confused as to how to do so.

So, I know that the N-point IDFT can be expressed in terms of the forward DFT by:
$$
IDFT = \frac{1}{N}[DFT(X_n^*)]^*
$$

Also, an N-point DFT can be implemented with two N/2-point DFTs by performing a DFT on the even-indexed samples and another on the odd-indexed samples.

How would I go about combining these concepts?

Just follow your recipe as written

Conjugate your frequency domain vector
Apply DFT on even samples
Apply DFT on odd samples
Recombine the two $N/2$ DFTs to create the $N$ point DFT result
Conjugate the DFT result and divide by N.

You can swap the last two steps, if you want.

Keep in mind that there is rarely a good reason to split an $N$ point DFT into two $N/2$ point DFTs. Because of the recombination it's typically more as expensive a single full size DFT. In essence you just take one stage out of the FFT algorithm and compute it manually. The only reason to go there if you have severe memory constraints or cache management issues.

That's different from using a $N/2$ point complex FFTs to compute an $N$ point FFT on real value data. That's generally a considerable savings.


Prove a property using shift theorem and duality I'm reading Lectures on the Fourier Transform and Its Applications and I'm going to prove shift theorem for the inverse Fourier transform using duality. According to the mentioned source, the duality equations are
$$\mathcal{F}f = (\mathcal{F}^{-1}f)^{-} \tag{1}$$
$$\mathcal{F}f^- = \mathcal{F}^{-1}f \tag{2} $$
$$ \mathcal{F}\mathcal{F}f = f^{-} \tag{3}$$
Where $$(\mathcal{F}f)(s) = \int_{-\infty}^{+\infty}e^{-2\pi i st}f(t)dt$$ $$(\mathcal{F}^{-1}f)(t) = \int_{-\infty}^{+\infty}e^{+2\pi i st}f(s)ds$$
And $$f^{-}(t) = f(-t)$$
Using shift theorem which says $$f(t) \leftrightarrow F(s) \implies f(t-b) \leftrightarrow e^{-2\pi isb}F(s) \tag{4}$$ And duality equations I want to show $$F(s) \leftrightarrow f(t) \implies F(s-b) \leftrightarrow e^{2\pi itb}f(t) \tag{5}$$
First of all $\text{(5)}$ seems ambiguous to me but I guess it's $$\mathcal{F}\{ e^{2\pi itb}f(t)\} = F(s - b) \tag{6}$$
It's straightforward to show that directly $$\mathcal{F}\{ e^{2\pi itb}f(t)\} = \int_{-\infty}^{+\infty}e^{-2\pi i st}e^{2\pi itb}f(t)dt = \int_{-\infty}^{+\infty}e^{-2\pi it(s-b)}f(t)dt = F(s-b)$$
But I don't know how to combine duality equations and shift theorem to prove $\text{(5)}$. There is another version of duality which can be found for example here. It says $$x(t) \leftrightarrow X(s) \implies X(t) \leftrightarrow x(-s) \tag{7}$$
Is $\text{(7)}$ same as duality equations i.e. equations $\text{(1)}, \text{(2)}$ and $\text{(3)}$ or they are different from each other?
Edit: Equation $\text{(6)}$ is called modulation theorem. Also it seems in $\text{(5)}$ by $F$ we mean $\mathcal{F}^{-1} f$. So it's completely different from $\text{(6)}$. I'm really confused by these similar equations.
Edit2: I realized that the other version of duality is same as $\text{(3)}$. I mean the equations $\text{(3)}$ and $\text{(7)}$ are identical.

If $X(f)=\mathcal{F}\{x(t)\}$ you have, in analogy with Eq. $(7)$,
$$x(f)\Longleftrightarrow X(-t)\tag{1}$$
where we assume that functions with the independent variable $f$ are Fourier transforms of the corresponding functions with independent variable $t$.
Applying $(1)$ to
$$x(t-b)\Longleftrightarrow e^{-2\pi jbf} X(f)\tag{2}$$
gives
$$x(f-b)\Longleftrightarrow e^{2\pi jbt} X(-t)\tag{3}$$
According to $(1)$, $X(-t)$ is the inverse Fourier transform of $x(f)$, so with $g(t)=X(-t)$ and $G(f)=\mathcal{F}\{g(t)\}=x(f)$ you obtain the desired result:
$$G(f-b)\Longleftrightarrow e^{2\pi jbt} g(t)\tag{4}$$


Root mean square bandwidth of a signal with flat spectrum I am trying to derive a simple expression for root-mean-square bandwidth $B_{\rm rms}$ for a signal that has flat spectral density over its entire bandwidth $B$. The expression for the $B_{\rm rms}$ is

$$
B_{\rm rms} = \sqrt{\frac{\displaystyle\int_{-B/2}^{B/2} \lvert H(\omega)\rvert^2 \omega^2\, \mathrm{d}\omega}{\displaystyle\int_{-B/2}^{B/2} \lvert H(\omega)\rvert^2\, \mathrm{d}\omega}}
$$

where $\lvert H(\omega)\rvert^2 = A$ is the constant PSD over $B_\rm h$.

I read it at some sources that this gets simplified as follows:
$B_{\rm rms} = \frac{B}{\sqrt{3}}$.

However, when I integrate this I do not see how I can get this expression:

$$
B_{\rm rms} = \sqrt{\frac{\displaystyle\int_{-B/2}^{B/2} A \omega^2\, \mathrm{d}\omega}{\displaystyle\int\limits_{-B/2}^{B/2} A\, \mathrm{d}\omega}} = \sqrt{\frac{\displaystyle\frac{1}{3}\cdot\frac{2B^3}{8}}{\displaystyle\frac{2B}{2}}} = \frac{B}{\sqrt{12}}
.$$

Your formula for the RMS bandwidth makes sense for (perfectly band-limited) low pass signals, i.e., for signals with a spectrum centered around $f_0=0$. The bandwidth of low pass signals is defined as the support of their spectrum at positive frequencies. So your integration limits must be $-B$ and $B$. This results in

$$B^2_{\rm rms}=\frac{\displaystyle\int_{-B}^{B}f^2df}{\displaystyle\int_{-B}^Bdf}=\frac{2B^3/3}{2B}=\frac{B^2}{3}\tag{1}\\$$

which is the expression that you're looking for.

Note that the RMS bandwidth is also defined for signals that are not ideally band-limited. In that case you have to integrate from $-\infty$ to $\infty$, as pointed out in MBaz's answer.

If you have a (not necessarily perfectly band-limited) band pass signal centered around $f_0\neq 0$, you must use the following formula for the RMS bandwidth:

$$B^2_{\rm rms}=\frac{4\displaystyle\int_{0}^{\infty}\lvert H(f)\rvert^2(f-f_0)^2df}{\displaystyle\int_{0}^{\infty}\lvert H(f)\rvert^2df}\tag{2}\\$$

If $\lvert H(f)\rvert$ is perfectly band-limited and if it is constant ($\lvert H(f)\rvert=c$) in the interval $[f_0-B/2,f_0+B/2]$, its RMS bandwidth is

$$B^2_{\rm rms}=\frac{4c^2\displaystyle\int_{f_0-B/2}^{f_0+B/2}(f-f_0)^2df}{c^2\displaystyle\int_{f_0-B/2}^{f_0+B/2}df}=\frac{4c^2\displaystyle\int_{-B/2}^{B/2}f^2df}{c^2B}=\frac{B^2}{3}\tag{3}\\$$

which is the same result as for the corresponding low pass signal.

Note the factor $4$ in Eq. $(2)$, which is caused by the fact that the RMS bandwidth (without the correction factor) measures the "width" (variance) of the function $|H(f)|^2$. For low pass signals this width is measured across the whole spectrum (positive and negative frequencies), whereas for band pass signals, this width is measured only at positive frequencies. If this width is the same for a low pass signal and for a band pass signal, then the bandwidth of the band pass signal must be twice the bandwidth of the corresponding low pass signal, because bandwidth is measured at positive frequencies only. This is the reason why the squared bandwidth in the case of a band pass signal needs a factor $4$ to comply with the standard definition of bandwidth for low pass and band pass signals.


Unexpected imaginary part in the fft of a zero-padded cosine I simulated a cosine waveform $ y = \cos(\omega t) $ and I applied the FFT algorithm to it. As expected, I have a frequency peak at $\pm \omega$ only in the real part, and nothing in the imaginary part. As expected, if I zero-pad the array, the peak in the real part becomes a $\mathrm{sinc}(\omega D)$ (where D is the length of the window). However, in the imaginary part, a sort of anti-symmetric $\mathrm{sinc}$ appears. I don't understand why, because I was expecting the imaginary part to stay zero since the Fourier transform of the cosine is purely real.
Why is this happening?
Here a plot of the imaginary and real parts of the FFT of the described signal. In orange: FFT of the original signal. In blue: FFT of the zero-padded signal.

Here is my code, written in MATLAB language.
dt = 0.02;
fs = 1/dt;

start = 0;
stop = 1;
t = start:dt:stop-dt;

w = 2*pi*10;
y = cos(w*t);

N = 10000;
y_pad = [y zeros(1,N)];

Y_PAD = fftshift(fft(y_pad))./length(y);
df_pad = fs / length(y_pad);
f_pad = (0:(length(y_pad)-1))*df_pad;
f_pad(f_pad >= fs/2) = f_pad(f_pad >= fs/2) - fs;
f_pad = fftshift(f_pad);

Y = fftshift(fft(y))./length(y);
df = fs / length(y);
f = (0:(length(y)-1))*df;
f(f >= fs/2) = f(f >= fs/2) - fs;
f = fftshift(f);

figure;
subplot(2,1,1);

plot(f_pad, real(Y_PAD));
hold on
plot(f, real(Y));

ylabel('Re(fft)')
xlabel('f (Hz)')

subplot(2,1,2);

plot(f_pad, imag(Y_PAD));
hold on
plot(f, imag(Y));

ylabel('Im(fft)')
xlabel('f (Hz)')

Edit: By zero-padding only to the right of my signal I am breaking the even parity of the cosine and $y[N-n]=y[n]$ does not hold anymore (see Hilmar's answer). Therefore, I tried to zero pad both to the left and to the right, in order to preserve the parity
N = 10000;
y_pad = [zeros(1,N) y zeros(1,N)];

The result is the following and it looks weird: the imaginary part is still not zero and the real part oscillates between positive and negative values. Why this zero-padding does not work as expected?

A signal that's discrete in frequency is periodic in time. By using an FFT you inherently make this assumption and the FFT calculates the spectrum of a signal that's infinitely repeated in time.
The FFT is ONLY real valued if the time domain signal has even symmetry, i.e. $x(-t) = x(t)$. For a discrete periodic signal of period N this becomes $x[n] = x[N-n]$. That holds true for a pure cosine (at least for some frequencies), but not for the zero padded version.


Why does the amplitude of a discrete convolution depend on the time step? I've just finished covering convolutions in my signals class, and I've been playing around with the conv function in MATLAB, but there's something I don't quite understand. Say I have two discrete signals, $a$ and $b$ sampled at $dt$. As $dt \rightarrow 0$, the amplitude of the result increases. Since the convolution integral is just
$$
y(t)=\int_{-\infty}^\infty f(\tau)h(t-\tau)d\tau,
$$

why does the resulting convolution seem to depend on what my value is for $dt$? By decreasing the sampling rate, shouldn't I just have a better approximation?

Your computer doesn't compute the continuous integral, it does discrete convolution, which is just a sum of products at each time step.  When you increase dt, you get more points in each signal vector, which increases the sum at each time step. You must normalize the result of conv() according to the length of the vectors involved.


Circle detection opencv I Would like an explanation of the parameters for opencv's HoughCircles function.

I'm new to image processing in general and I have recently started using opencv. I'm struggling to find circles. I obtained the image below by converting a source image to greyscale then performing a binary threshold.

To my eye, the leftmost images are almost perfect circles. I have seen demonstrations where, what I assume are very noisy and distorted circles have been identified by opencv's hough transform.

Assuming the image displayed above is saved as "binaryresult.png". If I run the following code I get circles that are poorly placed. I'd like to understand the HoughCircles function's parameters so that I can better detect circles.

import cv2
import numpy as np

image = cv2.imread("binaryresult.png")
image = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)
circles = cv2.HoughCircles(image, cv2.cv.CV_HOUGH_GRADIENT,  2, 32, param1=200, param2=100)

# ensure at least some circles were found
if circles is not None:
    # convert the (x, y) coordinates and radius of the circles to integers
    circles = np.round(circles[0, :]).astype("int")

    # loop over the (x, y) coordinates and radius of the circles
    for (x, y, r) in circles:
        # draw the circle in the output image, then draw a rectangle
        # corresponding to the centre of the circle
        cv2.circle(output, (x, y), r, (0, 255, 0), 4)
        cv2.rectangle(output, (x - 5, y - 5), (x + 5, y + 5), (0, 128, 255), -1)

    # show the output image
    cv2.imshow("output", np.hstack([image]))
cv2.waitKey(0)

I'll start with the Hough line transform.

The hough transform works by creating a buffer (2 dimensional in our case) which represents all the possible lines in the image.

Any possible line can be represented with those 2 numbers $(\rho, \theta)$

From wikipedia

Note this uses polar coordinates to avoid division by zero issues rather than your typical $y = mx + c$ for a line

The hough transform scans an image and casts votes into the 2 dimensional buffer corresponding to $(\rho, \theta)$ values.

The buffer slots that accumulate the most votes are likely to be attributed to the presence of a line for that $(\rho, \theta)$ value.

That's the basic idea of the Hough transform but it can get more complicated from here.

For selecting $(\rho, \theta)$ values from the buffer, methods can get more sophisticated, such as using kernel density estimation. Although I don't think opencv does that.

In the case of the Hough circle transform, we now have a 3 dimensional buffer.
$(x, y, r)$

Where

$x, y $:  coordinates of center of circle

$r $: radius

Now in practice there's more to making the Hough transform work than this.

Take a look at the source code for the icvHoughCirclesGradient function.

Its here

Also cvHoughCircles and cv::HoughCircles

you can see that the function is also doing canny edge detection and sobel filtering.

These operations have their own input parameters which you can also see as inputs to the hough function.

Examining the function signature

void cv::HoughCircles( InputArray _image, OutputArray _circles,
int method, double dp, double min_dist,
double param1, double param2,
int minRadius, int maxRadius )

Method: There seems to be only on method now, probably left this for future development

dp: this is the resolution of the votes buffer

min_dist : is the minimum distance between circles (ie: circles that are too close must be part of the same circle)

minRadius and maxRadius are just that.

param1 and param2 are a bit more cryptic

but if you spot the lines

int canny_threshold = cvRound(param1);
int acc_threshold = cvRound(param2);

You'll see these are related to the canny edge detection and the accumulator threshold for selecting a circle from votes.


Calculating N-point DFT of a signal based on another signal's DFT Recently I have asked this answer.

Now I would like to know a little bit more about expressing N-point DFT's of signals in terms of one another.

Having N-point DFT X(k) of a certain signal x(n), how can I calculate N-point DFT of a signal $x_{s}=x(n)+(-1)^n \cdot{}  x(n)$ . Assuming $N$ is even.

Having thought about it a little bit, I came to a conclusion that we cancel half of the samples out and multiply the value of the rest of samples by a factor of 2. The spectrum will now have greater amplitude of peak but I am not sure how it will look like in terms of the location of frequency peak.

You just need to apply the formula consistently to see what's happening. Noting that $(-1)^n=e^{-j\pi n}$ you get for the DFT of $x[n](-1)^n$

$$\begin{align}\sum_{n=0}^{N-1}x[n](-1)^ne^{-j2\pi nk/N}&=\sum_{n=0}^{N-1}x[n]e^{-j\pi n}e^{-j2\pi nk/N}\\&=\sum_{n=0}^{N-1}x[n]e^{-j2\pi n(k+N/2)/N}\\&=X[k+N/2]\tag{1}\end{align}$$

where $X[k]$ is the DFT of $x[n]$, and where I've assumed that $N$ is even. Consequently, the DFT of $x[n]+(-1)^nx[n]$ is given by

$$\begin{align}\text{DFT}\{x[n]+(-1)^nx[n]\}&=\text{DFT}\{x[n]\}+\text{DFT}\{(-1)^nx[n]\}\\&=X[k]+X[k+N/2]\tag{2}\end{align}$$

where the index $k+N/2$ has to be taken modulo $N$ if $X[k]$ is defined for $k\in[0,N-1]$.


Treating the measured frequency response of a system as a DFT I have frequency response data of an analog system that I would like to turn into a digital filter in Matlab.  Is it kosher to think of my frequency samples as constituting a DFT of the system impulse response, even though there was never a time-domain signal in the first place?  If I do an ifft on the frequency response data, can this be treated as the impulse response of the system?  If so, what is the sample rate?  Do I assume it is twice the highest measured frequency?

Is there a way to build a digital filter from impulse response data without knowing the sample rate?

So, after worrying about whether the frequency response is "low pass in nature" (i don't think it needs to be a LPF, it could be high pass in nature, but all of the interesting features of your frequency response should be below Nyquist and maybe well below Nyquist), you might also wonder about if you have the necessary phase information in your frequency response.  if it is magnitude-only data, then you might have to guess at the phase.

if you iFFT the magnitude-only data (mirroring it about $f=0$), you will come up with an impulse response that is also mirrored about $t=0$.  that's fine, give it a delay of half of the length of the impulse response, and you have a linear-phase FIR.  you get linear phase because you started not knowing anything of the phase.  you can also take that FIR, factor it with some nasty factorization program (maybe MATLAB has that), reflect (using the reciprocal function) all zeros outside the unit circle to inside the unit circle, the re-multiply all of the factors and you will have a minimum-phase FIR with the same magnitude frequency response.


Coding a Variable Delay with LFO (MATLAB/C) I am trying to code a Variable Delay (chorus/flanger) effect in MATLAB using a paradigm that would be friendly for porting to a lower level language like C.  What I currently have is a working echo/delay that uses read/write indices to write into a buffer.  I'm now trying to use an LFO (low freq sine wave) to modulate the delay time but can't seem to figure it out.  Sometimes my errors sound "cool" but they aren't correct or what I was trying to do.

% Create LFO
lfo = abs(sin(2*pi*lfo_rate*[0:length(output)]/fs));

    % Perform the Echo Loop
    for i=1:length(x_t)

    % Write Sample Into Delay Buffer
    buffer(writer) = x_t(i);

    % Get All Delay Into Output
    if i <= length(x_t)
        output(i) = x_t(i) + (b_n * buffer(reader));
    elseif i > length(x_t)
        output(i) = b_n * buffer(reader);
    end

    % Circular Buffer
    writer = writer + 1;
        if writer > length(buffer)
            writer = 1;
        end
    reader = reader + 1;
        if reader > length(buffer)
            reader = 1;
        end
    end

I've been trying to use the LFO Vector to modulate the readers position, but I'm pretty much stuck and the stuff I can find online about this sort of thing isn't too helpful either.

I answer this way because i don't have permission for comments.
I don't have a complete solution. What happens with an chorus flanger is this.

Signal xt goes in split in two(or maybe more) signal paths and mixed together at the end.
The difference is the length of the delay (Flanger <10ms and chorus 20 - 60ms).
First your Lfo discription isn't right.
It should me something like :lfo = A *sin(2*pi*lfo_rate*[0:length(output)]/fs) + Offset; Where A is half the time (Positive and negative) and offset the center time where the LFO is turning around.
And `A =< offset' otherwise the time gets negative

The difficulty is i think a variable buffersize. let say a chorus center time is 50 ms with ±10 ms delay. no Lfo you have 44100 * 20/1000 = 2205 samples. the max buffer = 2646 samples and the min buffer = 1764 Samples.

Hope you can do something with the info


Autocorrelation: numpy versus FFT I have a series a of values (0 and 1) coming from a Brownian process with drift for which I am studying the autocorrelation.

I used two methods:

1) numpy autocorrelation:

corr = np.correlate(a,a,mode='full')/a.size
corr = corr[corr.size//2:]

2) Fourier transforms (Wiener-Khinchin):

A = np.fft.fft(a)
S = np.conj(A)*A/a.size
c_fourier = np.fft.ifft(S)

However, I get different results:

I am not very experienced with stochastics or signal processing, so I have some difficulties understanding where the difference comes from.

As this page https://ccrma.stanford.edu/~jos/ReviewFourier/FFT_Convolution.html mentions, when you perform cyclic convolution (FFT convolution) you have to add zeros to your signal $A$ to the length $$N_Y = N_A + N_A -1,$$ where $N_Y$ is the length of the output convolution and $N_A$ is length of your signal. In your case you are performing the FFT Correlation, which is almost the same operation and you should extend your signal $A$ with zeros as well.

So for example:

a = np.concatenate((a,np.zeros(len(a)-1))) # added zeros to your signal
A = np.fft.fft(a)
S = np.conj(A)*A
c_fourier = np.fft.ifft(S)
c_fourier = c_fourier[:(c_fourier.size//2)+1]

After that you should get the same results as for numpy auto-correlation.


OFDM channel estimation As we know, if we want to get the time domain response of the channel $h[n]$, we can first get the frequency domain response of the channel $H[k]=Y[k]/X[k]$. And then we do IFFT.

But now one OFDM symbol includes some null subcarriers. So we don't need to get the $H[K]$ on these null subcarriers. But under this circumstance, how do we get the time response of the channel $h[n]$?

If some values in $X[k]$ are zero, you cannot determine the impulse response of the complete channel bandwidth ($-f_\mathrm s / 2 \ldots f_\mathrm s / 2$, with sampling frequency $f_\mathrm s$). The IFFT of the estimated $\hat H[k]$ is rather a low pass filtered version of the channel. This might or might not be what you want.

In reality, you can never measure the complete bandwidth of a channel, as it is virtually infinitely large. For a communication system it is usually sufficient to know the channel transfer function of the band in which you intend to transmit a signal. If you transmit only on a fraction of subcarriers of an OFDM system you're effectively narrowing the transmission band and you might be satisfied with knowing the impulse response of this narrower band.


Aperiodic discrete signals What is the fourier transform for finite aperiodic discrete signals and for infinite aperiodic discrete signals?

Are they different?

The discrete-time Fourier transform (DTFT) is just what it is, regardless of the length of the signal:

$$X(\omega)=\sum_{n=-\infty}^{\infty}x[n]e^{-jn\omega}\tag{1}$$

For finite length signals you just have to sum over finitely many elements in $(1)$, which means that the DTFT always exists. This is not necessarily the case when the signal's length is infinite, because then the sum in $(1)$ does not necessarily converge.


Feature extraction from spectrum to aid categorisation i have two sets of signals whose spectrum looks like shown in the figure. The first column of signals belongs to category 1 and the second belongs to category 2.
Looking at the spectrums I want to identify features which will differentiate the two.
Eg:mean,variance etc.

First thing to mention is that you must have some kind of Machine Learning algorithm to perform your recognition. You might want to use the following - what's yours?

Neural Networks (easy to implement)
SVM's (should perform best)
k-NN (eventually, but needs lot's of data to store)

Regarding features, it's hard to determine which to use for your application just by looking at spectra - both of categories have peaks at lower and higher part of range. Also there might be some significant information in time domain or a specific method valid for processing this type of signals. You never said what these signals are and what's the application.

The only thing I notice is that class 2 has some harmonic content occurring at frequencies $160 \mathtt{Hz}$, and $200 - 225 \mathtt{Hz}$. I don't think that simple mean and variance will do. Probably you should try different and choose best performing ones. My suggestions of features you've asked for (to start with), are:

Spectral Slope (what's the gradient of the linear regression of your spectrum)
Audio Spectrum Flatness (calculated in frequency banks - tells you if signal is noisy or harmonic)
Audio Spectrum Centroid (kind of 'mean value' you mentioned)
Audio Spectrum Envelope (can be understood as very general spectrum descriptor - envelope)
MFCC (well defined and described - they are describing your spectrum by using cepstral analysis)

For more info about implementations you can refer to this great book (easy to get): MPEG-7 Audio and Beyond.


Why is the impulse response function of this system 0? Suppose I have an system $ y(t) = t^{2}x(t)$.

The impulse response of this system would be: $h(t) = t^{2} \delta(t)$.

Since $\delta(t) = 0$ for $t \neq 0 , h(t) = 0$ for $t \neq 0$.
And at $t=0, h(t) = 0$.

But since $y(t) = \int_{-\infty}^{\infty}x(r)h(t-r)dr $, wouldn't $y(t) = 0$?

What am I missing here?

This system
$$ y(t) = t^2 x(t) $$

is not LTI and therefore does not have an impulse response of the form $h(t) = \mathcal{T}\{\delta(t)\}$.

So your statement $h(t) = t^2 \delta(t)$ is not correct... Hope this solves your confusion.


Median Filter one after another Operate on an image by performing Median Filtering in a 3x3 window.
Operate on the resulting image by performing, again, Median Filtering in a 3x3 window.
Can the resulting image be obtained from a single Median filtering?

my initial thought is that it can be done with the right mask. maybe a median next to a median. but i'm not sure.

The answer is no:

See the following rearrangements of numbers from 1 to 81.
In the left  case the 3x3 median of 3x3 medians will be

median([77, 72, 67, 62, 57, 52, 47, 37]) = 57,

in right case
median

([5, 15, 23, 32, 41, 50, 59, 68, 77]) = 41

While the median of the whole 9x9 block is 41 for both cases.


Can you calculate the absolute/amplitude values of a Fourier transform without complex/imaginary? The DFT formula is $$X(k) = \sum x(n)e^{-i2πkn/N}$$

Is it possible to obtain absolute values from the DFT without having to calculate the imaginary components? I only need the signal amplitude values, not phase. I'm trying to create a spectrogram so I'm wondering if there's a shortcut or if I have to program complex math into the programming language

Yes.  You can do the calculation.  But, no, it's not a shortcut.

If you want to obtain your DFT magnitudes fast (e.g. using the O(NlogN) FFT algorithm), then you will need the complex math (or its computational equivalent).  This is because the cosine and sine forms of the DFT basis vectors are orthogonal, and you need both to completely be able to represent arbitrary input waveforms.  The sine portions of the basis vectors are represented by the imaginary components of an FFT's complex result.  So the full complex result of an FFT needs to be computed to end up with enough information for a magnitude or amplitude result.

If a slow DFT computation is acceptable, then you could just compute the cosine components using strictly real arithmetic and then separately compute the sine components also using reals.  Then compute the hypotenuse of the real triangle for the amplitude.  No complex numbers needed.  But this method is O(N^2), which would take a lot longer.

For just magnitude results (computed slowly), N Goertzel filters is another possibility, and may save in trig function calls.


Are there any good "conceptual distances" for comparing two images? I want to define some distances to use over an image, which I create.

I think I need to define two different "distances".

The first distance should be the most robustness preserving in compare to the original image.
The second one should would be the most transparency preserving in compare to the
original image.

For example:

Let $I_1, I_2$ and $I_3$ be images, $I_1$ is the original and $I_2, I_3$ are some similar images:

I would like to define a metric or some distance so a and b will holds:

a. $|I1 -I2| > |I_1 -I_3| <=> I_2 $ is more robustness preserving then $I_3$ in compare to
the origin $I_1$.

b. $|I_1 -I_2| > |I_1 -I_3| <=> I_2$ is more transparency preserving then $I_3$               in compare to the origin $I_1$.

as I said the distance in a. and b. don't have to be the same at all.

I saw this post but it didn't help me alot, it may help you to find a good answer:

What distance metric can I use for comparing images?

Ok, I read the question posted in your link, so forget about MSE or SSIM.
What you can do is define a specific distance for the particular type of data you are using.
For example, if you are using video data, depending on what you are most interested in, you can define distance measure based on the average shot length, keyframe colour histogram, average motion in a shot, and so on.
If you are working on a database of face images, you can use e.g. PCA+eigenfaces to project each image on a N-dimensional space and compute an euclidean distance on that space; if you are working with e.g. flower pictures you can define a distance based on edge + colour information...

Here´s my previous answer:

It´s not really clear to me what you mean for robustness/transparency, to be honest.

If you want to check which one of I_2 and I_3 is more similar to I_1 a simple way is to use compute the mean square error, defined as the sum of the squares of all the pixels in the image, normalized for the total # of pixels


what is meant by biased and unbiased Autocovariance sequence estimate? In nonparametric spectral estimation method what is meant by biased and unbiased Autocovariance sequence estimate and why most commonly biased one is used?
I have the answer to the second question from Peter Stoica (SPECTRAL ANALYSIS OF SIGNALS)but not getting what he wants to say.

In spectral estimation, one estimates "spectra" and other values. And one is interested in knowing whether the estimate $\hat{x}$ converge to the "true" but unknown value  $x$. It is unknown, for instance because one has only access to a finite quantity $N$ of samples to compute with, a limited number of realization, because ergodicity is not assured, etc. So we want to know whether  $\hat{x}_N\to x$ as $N\to \infty$. The convergence can be thought in different ways, and the most common is the mean square sense: does

$$\lim_{N\to \infty} E(|\hat{x}_N-x|^2)=0$$
where $E$ is the expectation? This can be really complicated to prove in practice. Thus, often people are satisfied when simpler conditions are met. And there are two necessary conditions for the above:

$$\lim_{N\to \infty} E(\hat{x}_N)=x$$

$$\lim_{N\to \infty} E(|\hat{x}_N|^2)=0$$

The first is the asymptotic unbiasedness, the second the vanishing variance. Together, they define consistent estimators. So, a biased autocorrelation estimate won't converge to the true autocorrelation. And derived estimators (like periodograms) are unlikely to converge as well.


Long Term Average Spectrum of large batch of audio files I'm trying to obtain the overall LTAS of a large batch of audio (wav) files. To do this, I should obtain the LTAS of each file individually, then take the average value across all LTASs for each frequency. I also need the standard deviation for each frequency. I'm doing this in Matlab, but am rather a noob. Can anyone offer any pointers to help me write a Matlab script to do this? Thanks!

use "wavread()" to read the audio files
use "pwelch()" to calculated the power spectrum of the songs
you do energy averagin over multiple songs.

Something like this (just a code sketch)

nx = 8192; % that's about 5 Hz @44.1 kHz
dd = dir('*.wav'); % all files in the current directory
nFiles = length(dd);
LTAS = zeros(nx/2+1,nFiles)
for i = 1:nFiles
  x = wavread(dd(i).name);
  LTAS(:,i) = pwelch(x,hanning(nx),nx/2);
end
LTASOverAllSongs = mean(LTAS,2);


How do ADC/DACs deal with out of sync signals near Nyquist frequency? Let's say a converter is sampling continuously at 48 kHz and the incoming waveform is a sine wave at 24 kHz... How does it faithfully recreate the waveform if the samples are taken at points that are 0?

How does it faithfully recreate the waveform if the samples are taken a points that are 0?

Not at all. Things at $\frac{f_\text{sample}}2$ might just as well be constants. You can't tell [0 0 0 0] (sampled 24 kHz signal at zero-crossings) from [0 0 0 0] (sampled zero constant signal).

Nor can you tell [0.1 -0.1 0.1 -0.1] (sampled 24 kHz signal with amplitude 1 just a bit after each zero-crossing) from [0.1 -0.1 0.1 -0.1] (sampled 24 kHz signal of amplitude 0.1, sampled at its extrema!).

You need the sampling rate to be higher (not only just as high) than two times the highest signal frequency component for this to be unambiguous.


Making images taken under daylight conditions appear as if they were taken at night using image processing elementary techniques I know that these days' style-transfer GANs should be able to style-transfer day images into night images. I however prefer to do it using a set of elementary image processing manipulations, so I can control their random parameter ranges.

I have a small dataset of images taken at night, which are very hard to train on with transfer-learning from an ImageNet pre-trained network. This is kind of making sense, as ImageNet images are usually taken in daylight conditions.

I therefore want to try to create a "night version" of ImageNet, by applying elementary image-processing manipulations on it (e.g., adding noise, shrinking the color histogram and reducing its mean, etc.). That way, I hope I will be able to train a network that will be useful for transfer-learning on night images.

Are there known algorithms/methods for doing that daylight->night conversion using image processing?

That really doesn't make sense: you'd be training a neural network to mimic a conversion algorith, which you already have. That's a waste of electricity ;) I doubt the resulting nets would generalize at all: you're not even training them with actual day/night picture pairs. So, really, as usual: one of the most significant problems in machine learning is getting a good set of input data.

In your case, that's still pretty easy (just take a couple hundred pictures from the same positions at day and night – you can basically do that in a week, if you have a camera that you can reliably bring to the same positions). So, don't skip on that crucial step.


Resampling scattered data with MATLAB's $\tt interp$ and $\tt resample$ I am recording acceleration data with an MPU6050 connected to a Arduino1 and stored on a SD. Here you can find the code.

I need to calculate the FFT of an acceleration signal that was not sampled uniformly, so i have to resample my signal. After this suggestion I tried to study how to resample the signal.

From Mathworks documentation here it seems that I have to know the frequency of my signal and a nominal frequency. I read and studied the function resample, but I don’t understand how it works very well. In my specific case how to define

x= sin(2*pi*f*T)

But in my code below I tried to implement the function (to give you more information for a detailed answer and prove my efforts, but I don't have enough reputation to post all the graph).

Then I studied gettdata  “Scattered interpolant class” finding not very useful in my situation, they seems to work for 2 or more dimensions.

So I followed another way, I used the function interp1. I defined my equispaced vector using the mode and median value of my sample frequency, I discarded the average to avoid the influence of any peaks.

Is a correct approximation?
What kind of errors can occur?

This is my first time I approach to the FFT, and I have a civil engineer background, so is my first time with signal analysis, never studied Signal Theory before.

My interest for the FFT is to define the best low pass filter (example apply a Butterworth but I don’t know how to choose the filter order and cutoff frequency.

I don't have enough reputation to post more than 2 links sorry for that.

filename= uigetfile ('.txt');
fileID = fopen (filename);
logmpu6050 =csvread(filename);
fclose (fileID);
%Starting creating the specific Vectors
%Record Time in millisecond
time=logmpu6050(:,1);
%The x y z are converted to m/s^2
confactor=19.6/32768;
ax=logmpu6050(:,2);
confactor=19.6/32768;
ax=ax*confactor;
ay=logmpu6050(:,3);
ay=ay*confactor;
az=logmpu6050(:,4);
az=az*confactor;
%Define the sample rate subtracting from Sampletime i+1 Sampletime i
n=length(time);
for i = 1:n-1
samplerate(i) = time(i+1)-time(i);
end;

%I try to define the best frequency for "resample" function based on mode and median value
%I will not use average to avoid conditioning due to extreme value
f1=mode(samplerate);
f2=median(samplerate);
x1= sin(2*pi*f1*time);
x2= sin(2*pi*f2*time);
plot(time,x1,'.');
figure
plot(time,x2,'.');

%resampling at time mode record frequency
[test1,test2]= resample(az,time,f1);
%resampling at time median record frequency
[test3,test4]= resample(az,time,f2);
figure
plot(test2, test1)
figure
plot(test4,test3)

f1=int32(mode(samplerate));
f2=int32(median(samplerate));

%Try to implement the matlab function "interp1"
%My input are
%x-> time when the value is recorded "time"
%v-> the acceleration value "az"
%xq-> I tried two define two type of equispaced vectors
%xq1-> Defined by the total recording time divided by the mode
%xq2-> Defined by the total recording time divided by the median
totrectime= int32(logmpu6050(n,1)-logmpu6050(1,1)) %milliseconds
numsamples1= idivide(totrectime,f1); %Rounded toward zero to the nearest integers.
numsamples2= idivide(totrectime,f2);
f1double=double(f1);
f2double=double(f2);
for i = 1: numsamples1
xq1(i)= f1double*i;
end;
for i = 1:numsamples2
xq2(i)= f2double*i;
end;
xq1double=double(xq1);
xq2double=double(xq2);

vq1 = interp1(time,az,xq1double);
vq2 = interp1(time,az,xq2double);
figure
plot(xq1double,vq1);
figure
plot(xq2double,vq2);

In your question, you state that your ultimate goal is to compute the discrete Fourier transform (DFT) from nonuniformly-sampled input data. It is not necessary to resample your data at all.  In fact, there are several fast algorithms, known variously as the unequally spaced fast Fourier transform (USFFT) and the nonuniform fast Fourier transform (NUFFT) that do exactly what you want: take nonuniform samples as input and produce equally spaced samples of the DFT as output. The remainder of your question, dealing with issues related to MATLAB's resampling routines, are a red herring.

Here are some publications describing the USFFT algorithms:

G. Beylkin, On the Fast Fourier Transform of Functions With
Singularities,  Applied and Computational Harmonic Analysis, 2, pp.
363-381, 1995. Link to preprint.
A. Dutt and V. Rokhlin, Fast Fourier Transforms for Nonequispaced Data, SIAM J. Sci. Comput., 14, pp. 1368, 1993.

There are many commercial and free implementations of these algorithms. You can access some of them here.

EDIT: Here is a "quick and dirty" approach that may help you.

Based on your problem description and comments, you are just trying to get a sense of the spectral content of your signal, which does not require high accuracy. Here is one way to formulate your problem so that the spectrum can be easily estimated using any of the USFFT implementations in the above link.

You have a signal $f(t)$ defined on the $T$-duration interval $t\in[-T/2,T/2]$, and you wish to estimate its Fourier transform
$$\hat{f}(\omega)=\int_{-T/2}^{T/2}f(t)e^{-i 2 \pi \omega t}\,dt$$
at a number of values of $\omega$. You have collected $N$ samples of your function
$$f_n=f(t_n),\quad n=1,\ldots,N$$
at the unequally spaced time locations $t_n$, $n=1,\ldots,N$.

In general, you would solve this by formulating a least squares problem for the Fourier transform $\hat{f}$, and using the USFFT to rapidly solve the system of equations, but since you do not require high accuracy, here is a "quick and dirty" solution.

Use the trapezoid rule from Calculus 101 to discretize the above integral, giving the approximation
$$\hat{f}(\omega) \approx \frac{\Delta_1}{2}f_1 e^{-i2\pi\omega t_1}+\left(\sum_{n=2}^{N-1}\frac{\Delta_{n-1}+\Delta_{n}}{2}f_n e^{-i2\pi\omega t_n}\right)+\frac{\Delta_{N-1}}{2}f_N e^{-i2\pi\omega t_N}$$ where
$$\Delta_{n}=t_{n+1}-t_n$$
is the temporal spacing between samples $n+1$ and $n$.

Define the new sequence
$$g_n=\begin{cases}\frac{\Delta_1}{2}f_1, & n=1 \\ \frac{\Delta_{n-1}+\Delta_{n}}{2}f_n, & n=2,\ldots,N-1 \\ \frac{\Delta_{N-1}}{2}f_N, & n=N\end{cases}$$
Now pick some frequency grid you care about. For example, say you want $2000$ equally spaced frequencies between $0$ Hz and $25$ Hz, and define
$$\omega_m = \frac{m-1}{2000-1}25,\quad m=1,\ldots,2000.$$

Now your problem amounts to evaluating the sum
$$\hat{f}(\omega_m) = \sum_{n=1}^N g_n e^{-i2\pi\omega_m t_n},\quad m=1,\ldots,M.$$

This looks like an FFT, except the $t$-values are unequally spaced and the $\omega$-values are unusual. If you do not care about speed, then just evaluate the sum directly using your favorite programming language.  It will take a minute or two of computer time, but so what? If you are just doing this a few times to get a sense of what your signal's spectrum looks like, then it does not matter how fast it is. If at some point in the future you decide that you do care about speed, then the sum can be evaluated extremely rapidly using a single call to a USFFT library, which you can download from the link I posted (or you can use one of the many commercial implementations).  Making the switch to the USFFT will only require you to modify a handful of lines of code.


Why Fourier Transform doesn't show obvious periodicity? I have a clearly 4Hz periodic signal (see image below). However when I apply a fourier transform, there is no peak at 4Hz. Also, when I do a power spectrum density (using Matlab's pwelch algorithm), there is no peak at 4Hz.

My original signal is sampled at 5ms, and is 141 bins (705ms) long.

Why is this & which alogithm will detect the 4Hz peak?

Thanks

An FFT doesn't necessarily show a fundamental periodicity, but instead spectral periodicities, which might be higher harmonics (or overtones) of the periodicity that you might perceive visually.

Also, some spectral peaks can be hidden under (the windowing artifacts of) a relatively  large enough DC offset.  If you subtract the mean of the data before doing the FFT, these frequency peaks might show up more clearly.


How to calculate bandwidth require for OFDM I read LTE PHY profile and do not understand the relationship between number of FFT and bandwidth. For example LTE downlink channel bandwidth 1.25MHz and 5MHz have 128 FFTs and 512 FFTs accordingly. Another example is 802.11 a bandwidth 20MHz, FFT size 52.

I design an OFDM system system at carrier frequency 910MHz, with FFT size for example 64, how much coherent bandwidth that the system need? What is the formula? What other factors should be taken into account?

Add >> 18 Aug (I don't have enough reputaiton to add comment on Jason R reply)

When I test a simple software defined radio OFDM system with parameters as FFT = 256, data carrier per symbol = 125, number of pilot = 25, bandwidth = 5MHz, carrier frequency = 910MHz, the channel estimation wouks fine. If I change bandwidth to 1MHz, then the channel estimation goes wrong. In this case I reduce sampling frequency 5 times, so carrier spacing is reduced 5 times too.

I think the carrier spacing must be greater than a limit so that intersymbol interference (or inter-carrier interference?) is small enough to not cause channel estimation failture. How to find the carrier spacing limit?

Or delay spread $\tau \ll T$ where T is symbol interval may cause the problem. How to calculate the limit $\tau$?

Thanks

The channel bandwidth and FFT size alone don't provide enough information to describe the entire structure of the OFDM signal. Recall the following relationship:

$$
\Delta f = \frac{f_s}{N}
$$

where $\Delta f$ is the subcarrier spacing, $f_s$ is the sample rate used at the modulator input, and $N$ is the FFT size.

In your LTE examples, the two channel bandwidths must use the same subcarrier spacing, as it takes four times as much bandwidth to carry four times the subcarriers (512 versus 128).

For 802.11a, however, the FFT size is typically 64; with a sample rate of 20 MHz, this yields a subcarrier spacing of 312.5 kHz. Of these 64 potential subcarriers, 12 of them are unused (sometimes referred to as "virtual carriers"). Of the remaining 52 subcarriers, 4 contain pilot tones, so only 48 subcarriers carry information at any given time.

For your notional OFDM system, you can choose the amount of bandwidth required by your system, as required based on your throughput requirements. Choose the sample rate that you will input to the modulator, then choose a convenient FFT size. This defines your number of subcarriers and therefore the subcarrier frequency spacing.


Gaussian filter in terms of pixel radius So I'm writing an app that blurs and downscales large quantities of images, but all the image processing libraries I've come across define gaussian blurs in terms of sigma and kernel size rather than pixel radius (like an image editing app would do). Does radius relate to sigma and kernel size in some way, or is it an entirely different thing?

After a bit more digging I found a gaussian kernel calculator that seems to suggest sigma is roughly half the blur radius, and kernel size should be somewhere around (sigma*4)+1 for a full blur.


Convolution of two sets In my systems and signals course I had been asked a question about finding the convolution of two sets.
I was given:
\begin{align}
x[n] &= \{3,2,1\}\\
h[n] &= \{1,-2,3\}\\
\text{Find}\quad y[n] & = x[n]\star h[n]
\end{align}
I have no idea on how to find the convolution of two sets. Any help is appreciated.

If it is difficult for you to remember or calculate the convolution of two sequences then you may try doing it as polynomial multiplication.
Think of x[n] and h[n] as polynomial coefficients. So we have
Px = 3x^2 + 2*x + 1
Ph = 1x^2 - 2*x + 3

Remember that linear convolution of two sequences is polynomial multiplication. Therefore
Py = Px * Ph
Py = (3x^2 + 2*x + 1) * (1x^2 - 2*x + 3)
Py = 3x^4 + 2x^3 - 6x^3 + x^2 + 9x^2 - 4x^2 + 6x - 2x + 3
Py = 3x^4 - 4x^3 +6x^2 + 4x + 3

Now writing the polynomial coefficients back in sequence format
y[n] = {3, -4, 6, 4, 3}


Transmit power of non-orthogonal FDM It is known that the total power of an OFDM signal equals to the sum of the powers on each individual subcarriers. I am a bit puzzled by the total transmit power of a generalized FDM system if subcarriers are not orthogonal and if it should remain the same or changed due to the non-orthogonality.

Let us consider the following OFDM transmit signal:
$$s(t) = \sum_{k = 0}^{N-1} a_k e^{j2\pi f_k t} = \sum_{k = 0}^{N-1} a_k e^{j2\pi \frac{k \: \phi}{T} t}, \quad 0\leq t\leq T,$$

where $k$ is the subcarrier index, $a_k$ are i.i.d. data symbols with zero mean, $N$ is the total number of subcarriers, $T$ is the FDM symbol duration, and $\phi < 1$ is introduced to break the orthogonality.

The transmit power is represented as

$$ P = \frac{1}{T} \int_{0}^{T} s(t) \: s^*(t) dt,$$

which can be further simplified as

$$P = \frac{1}{T} \sum_{k = 0}^{N-1} \sum_{k' = 0}^{N-1} a_k a^*_{k'} \int_{0}^{T} e^{j2\pi \phi (k - k')t/T} dt$$

Now, the value of the integral $\int_{0}^{T} e^{j2\pi \phi (k - k')t/T} dt$ does not equal to $0$ if $k \neq k'$ (due to the non-orthogonality introduced by $\phi$). This means that the transmit power in case of non-orthogonal subcarriers is no longer equals the sum of the power on each subcarrier? However, it is higher than the sum!!!

I feel there is something wrong in this simple derivation, as by intuition, if a power source spends a certain amount of power on each subcarrier, then how come the total power will be greater than the sum of these individual powers?

========================

Update

The following update to make the modeling of non-orthogonality more clear. Let us consider the FDM signal written as

$$s(t) = \sum_{k = 0}^{N-1} a_k e^{j2\pi (f_k + \delta_k) t} = \sum_{k = 0}^{N-1} a_k e^{j2\pi (\frac{k \:}{T} + \delta_k) t}, \quad 0\leq t\leq T,$$

where $\delta_k$ is the frequency offset on subcarrier $k$. Following similar approach, the power $P$ can be found as

$$P = \frac{1}{T} \sum_{k = 0}^{N-1} \sum_{k' = 0}^{N-1} a_k a^*_{k'} \int_{0}^{T} e^{j2\pi (k/T + \delta_k - k'/T - \delta_{k'})t} dt$$

which does not equal to the sum of power on individual subcarriers.

In your example, the sinusoids are still orthogonal to one another; the key is that you must integrate over the correct period when calculating $P$. All of the sinusoids have frequencies that are multiples of $\frac{\phi}{T}$, so their resulting sum will have period $\frac{T}{\phi}$.

You should therefore take this fact into account when calculating the average power over a period of $s(t)$:

$$
P = \frac{\phi}{T} \sum_{k = 0}^{N-1} \sum_{k' = 0}^{N-1} a_k a^*_{k'} \int_{0}^{\frac{T}{\phi}} e^{j2\pi \phi (k - k')t/T} dt
$$

With this formulation, then you do get the useful property that the integral is zero for all $k \ne k'$, which allows you to simplify the expression to be simply the sum of the powers of each of the sinusoids.


Deconvolution Using Complex Division in The Frequency Domain Consider these two signals:

a = [1 1 0 0 0 0 0 0]
b = [1 0 1 0 0 0 0 0]

their convolution is

c = a * b = [1 1 1 1 0 0 0 0]

I am trying to obtain b by using complex division to divide the discrete Fourier transform of c by the discrete Fourier transform of a. I am aware that in general, there may not be a solution when attempting deconvolution in this way, due to division by zero issues etc. -- however, I want to understand why this doesn't work in this specific case.

The discrete Fourier transforms of all three signals are:

F(a) = [
    2                 (0)
    1.707 - 0.707 i
    1 - i
    0.293 - 0.707 i
    0                 (N/2)
    0.293 + 0.707 i
    1 + i
    1.707 + 0.707 i
]

F(b) = [
    2                 (0)
    1 - i
    0
    1 + i
    2                 (N/2)
    1 - i
    0
    1 + i
]

F(c) = [
    4                 (0)
    1 - 2.414 i
    0
    1 - 0.414 i
    0                 (N/2)
    1 + 0.414 i
    0
    1 + 2.414 i
]

The zero-th entry (marked (0)) is the sum of the number of 1 values in the signal, equal to $\Sigma_{i = 0}^7k_i$ for each kernel $k \, \epsilon \{a, b, c\}$; the value of the $N/2$-th entry (marked (N/2)) is equal to $\Sigma_{i = 0}^7k_i(-1)^i$.

The problem arises when F(c) is divided by F(a) at the $N/2$-th entry: this is 0/0, which is undefined, or is sometimes defined as evaluating to 1. However, the correct result for the complete division $F^{-1}(F(c) / F(b)) = F(a)$ can only be obtained if this instance of 0/0 in the $N/2$-th position evaluates to the value 2, not 1, since the $N/2$-th entry of F(b) is 2. Other than this value at the $N/2$-th position, all entries of F(b) can be recovered correctly by simply dividing entries in F(c) by the corresponding entry in F(a), using complex division.

Why is the N/2-th entry of the Fourier transform problematic in this way? Is there a robust way to derive what the correct result should be for the division of Fourier transforms at the N/2-th position? Or is this simply due to the fact that the value 2 is lost due to multiplication by 0 in the pointwise multiplication $F(a) . F(b)$?

Does this problem affect other components of the Fourier transform? Or is it only possible for this to happen at the N/2th position, which will be real-valued and may be zero?

Are there other gotchas to dividing signals in the Fourier domain in this way, or does this only happen in cases where there are zeroes in the Fourier domain of one of the factors?

It's not about the index, it is about the Filter :-).
Think of your coefficients as filters and what would happen to a Sine Signal with frequency of $ \frac{N}{2} {F}_{s} $ that would be filtered (Convolved) with your samples of signal $ a $.
You will sum samples with the same absolute value yet one is negative and the other is positive, namely their sum is zero which means nothing will be left from this signal.
In practice we never use complex division for Deconvolution.
The minimum you can do is use some kind of regularization.
Zeros in Frequency Domain
In the noiseless case, having zero magnitude at some indices means those indices were zero either in the input signal or the filter.
The problem is we can never know hence when we want to recover any of them the data is lost. This is due to convolution being multiplication in frequency domain so at some index $ {C}_{k} = {A}_{k} {B}_{k} $ and if either $ {A}_{k} $ or $ {B}_{k} $ is zero then $ {C}_{k} $ is zero without being able to recover $ {A}_{k} $ or $ {B}_{k} $.
It means, that, potentially if there are no zeros for the data or the filter (At least in their overlapping domain) data can be restored.
In practice, this is very hard problem given noise.