How can I insert message bits into the frequency spectrum of an audio signal using Matlab? I have been trying to use the DSSS spread spectrum technique for audio steganography in MATLAB. I multiplied the PN sequence with the message bits. Now I want to insert this spread message into a cover audio signal. In the time domain, I split the cover audio into blocks and just added the message. But I want to do it in the frequency domain, by taking the FFT of cover audio. How can I do that? Which is the better way?

The FFT is just an implementation of the DFT. There's, mathematically, no difference.

The DFT is a linear operation. So it makes absolutely zero difference whether you add up two signals in time or frequency domain.


Resonance and Pitch classes Suppose we denote the frequency of the note $\textbf{C}_{0}$ by $f$ then we know that $f,2f,3f,...$ are all resonating frequencies. However the frequencies $2f,4f,...,2^{m}f$ belong to the same pitch class $\textbf{C}$ having the "same" sound(higher pitch) whereas the frequency $3f$ is the note $\textbf{G}_{1}$ and belongs to the pitch class $\textbf{G}$. What is the reason(physical) for this ?

I think you gave the reason in your question. A periodic function can be written as a Fourier series, i.e. as a weighted sum of phase shifted sinusoids, where the frequencies of the individual sinusoids are integer multiples of the fundamental frequency. And the odd multiples obviously give you other pitches than the pitch corresponding to the fundamental frequency. This is what causes the same note to sound different when played on different instruments. It's called timbre.

The harmonics of the note C consist of all possible notes (and also in-between notes). The strongest harmonics with pitches different from C are G (the perfect fifth at $3f$), and E (the major third at $5f$).

EDIT: Your comment clarified your question for me. The phenomenon you're addressing is called octave equivalence. The most common argument for explaining octave equivalence is that the perfect octave is the most basic interval after the perfect unison, due to the simple relation between the harmonics of two notes an octave apart. But the assumption of octave equivalence is by no means universally accepted. As an example, every musician knows from experience that certain chord voicings which sound good in some range sound different (usually not as good) in a lower range. For some interesting thoughts on octave equivalence have a look at this page.


How to analyze multiple amplitude modulated signals? What I'm doing is that I send an audio signal consisting of a single frequency into a resonator. To check whether and how the signal has changed I perform a QAM demodulation.
I do not need intend to recover the signal or something like that, I'm just interested in the magnitude of the change/modulation of the input signal.
In the next step I would like to use multiple frequencies simultaneously and analysed how the amplitude of each signal has changed.
Is there a different/better way than just analysing the signal by Fourier transform?

Well, you can of course instead of using the Discrete Fourier Transform as filter bank just employ bandpass filters to extract your individual tones.
(The DFT is really just a filter bank, if you think about it: it correlates an input sample vector with different complex oscillations. Which is convolution with their conjugate time inverse, which is a filter operation.)
This might be advantageous, since the filter shape of a DFT is inherently sinc-shaped, which I guess isn't quite close to what you want in terms of spectral sidelobes.
So, yeah, simply design a narrow low-pass filter, and shift it in frequency to the appropriate position. (Do that for every tone.).
If you've got many tones that you space on a regular grid, you might want to do a polyphase filter bank, but I'd consider this an "advanced computational optimization that you probably don't need in your measurement process, as you don't care whether computations take a milli- or 100 microsecond".


Fractional Bandwidth of a Gaussian Amplitude Modulated Signal a gaussian modulated sinusoidal signal may be expressed as
$$x(t)=A\cdot e^{(j2\pi ft)}\cdot e^{\left[-\frac{1}{2\sigma^{2}}\cdot(t-t_{0})^{2}\right]}$$
Let's consider the case in which the gaussian peak is at 0 and the total amplitude is 1:
$$x(t)=e^{(j2\pi ft)}\cdot e^{\left[-\frac{1}{2\sigma^{2}}\cdot t^{2}\right]}$$

What is the fractional bandwidth of this function?
I'm asking this question since I think it's a quite used parameter. In fact, there are a proper function of MATLAB and a proper function of Python) which can generate this function, and both ask me to specify:

the sine wave frequency (which is $f$);

the signal fractional bandwidth

I think the last one is necessarily linked to the standard deviation $\sigma$. But by which mathematical relationship? It's a bandwidth in time domain (so the FWHM of the gaussian pulse: $2.35\cdot \sigma$) or the bandwidth in frequency domain (so ???)?

The fact that it's modulated with a sinusoid doesn't change the FWHM bandwidth of your pulse – the $e^{jx}$ function has $\left\lvert e^{jx}\right\rvert\equiv 1$ at every point. That doesn't change the amplitude, so the FWHM of a sinusoid-modulated gaussian is just the same as of the unmodulated gaussian.


Blind Carrier Synchronization From my reading, I see that there are two ways people usually do carrier synchronization:

Decision directed, where the transmitted symbols are unknown are the receiver, decisions are made, and then an error is computed. This assumes the receiver knows the type of signal was sent (in order to make the symbol decisions).
Data aided, where there is some known symbol sequence that is used for adjusting the phase and drive the error to zero.

Are there any carrier synchronization techniques which are non-data-aided and blind in the sense that all they require are IQ samples? I'm thinking of a scenario where I could give the technique a QPSK signal or a 16-QAM signal and the frequency offset would be fixed for both signals.

The best idea I have so far is to use a combination of signal identification and one of the two methods above (IQ samples $\rightarrow$ Signal ID alg. $\rightarrow$ Choose correct traditional technique). Maybe this is a good way to go about it, but if there is some more formal method, I am unsure.

Edit #1

Things that I know or have existing ways to estimate:

Symbol rate
Pulse shape
Bandwidth
Equally likely symbols
Modulation could be PSK or QAM
Assuming a single carrier system

Edit #2

The fourth power method is something that I known of but not in the way that is being described by Dan and Marcus. Apparently can be used for general constellations:

% make some random constellation by choosing 10 points within a 10-by-10 grid
constSize = 10;
iComponent = 10*(rand(constSize , 1)-0.5);
qComponent = 10*(rand(constSize , 1)-0.5);
const = iComponent + 1j*qComponent;
symbols = const(randi([1, length(const)], 1000, 1));

pulseFilter = rcosdesign(rolloff, span, sps, 'sqrt');
basebandSig = upfirdn(symbols, pulseFilter, sps);
basebandSig = exp(1j*2*pi*f0*[0:length(basebandSig)-1]'/sampleRate) .* basebandSig; % adding frequency offset f0 Hertz
basebandSig4 = basebandSig.^4;

% Fourth power estimate
fax = -sampleRate/2:sampleRate/nfft:sampleRate/2-1/nfft;
fftBasebandSig4 = fftshift(fft(basebandSig4));
[~, maxInd] = max(abs(fftBasebandSig4));
f0Hat = fax(maxInd)/4;

The problem is I don't fully understand what is happening. I understand for BPSK or QPSK how raising to the fourth power, all the points go to a single phase and then you can see the frequency offset but for the random constellation I generated, this doesn't happen.

To try and understand what is happening I plotted what happens to the signal when raising to the fourth power:

I see more points around (0, 0), but this also happens if I raise to the fifth power, as well as a peak at 5x the true frequency offset. What is special about the fourth power and do raising to other powers ever offer any benefit? (For a classic M-PSK constellation, you'd raise to the $M^{th}$ power)

Edit #3

In case of future readers, between the accepted answer and the paper, Blind estimation of frequency offset in the presence of unknown multipath, I was able to gain a better grasp of the "raising the oversampled signal to a power" method.

One approach not yet mentioned is frequency multiplication by a factor of the number of phase positions used followed by a PLL for noise reduction and then a frequency divider by then same factor.  The reason for this is when you multiply the frequency by N, you also multiply the phase by N. If the phase positions are selections of $2\pi/N$ (or those with any fixed offset), then multiplying the phases by N results in a consistent phase angle versus time, which then results in a fixed frequency tone that a PLL can lock to.

For example with BPSK (2 phase positions) a simple squaring will produce a specular component at 2x the carrier frequency since changing between 0° and 180° will result in changing between 0° and 0° (at twice the frequency) after doubling.  With QPSK (4 phase positions) you need to raise the signal to the 4th power (which will quadruple the frequency, and the phase).  The reason for the PLL is that any AM components (through pulse shaping or QAM modulation) will result in other lower level frequency components that the PLL will essentially filter out.

To answer your question as to what is so special about raising to the Nth power is that you basically raise to a sufficient enough power such that the modulo phase variation of the result is equivalent to small angle PM (which has a dominant carrier) in contrast to BPSK for example which has complete carrier suppression with equiprobable data.


Linkwitz-Riley Crossover Sum as Allpass Filter I was curious about crossover filter design, so I did some reading on Linkwitz-Riley filters.  Seems to me that the general idea is that if you add HP and LP filters and they are properly designed, you get an all-pass filter, which made sense on its face.  However, I cannot seem to work this out mathematically.  Using an LR 2, my understanding is that the filter transfer functions are simple biquads with a Q of 0.5.
$H_{LP}\left ( s \right )=\frac{\omega _{0}^{2}}{\omega _{0}^{2} + \frac{\omega _{0}}{Q}s + s^{2}} $
$H_{HP}\left ( s \right )=\frac{s^{2}}{\omega _{0}^{2} + \frac{\omega _{0}}{Q}s + s^{2}} $
And if they're summed you'd get:
$H_{LP+HP}\left ( s \right )=\frac{\omega _{0}^{2}+s^{2}}{\omega _{0}^{2} + \frac{\omega _{0}}{Q}s + s^{2}} $
Which is not an all pass filter unless Q is infinite.  Am I totally missing something here?

Just adding to the answer. The sign for recombining L/R filter alternates, so it's '-' for second order '+' for 4th order and so forth.
The algebra is isn't all that pleasant so I will only go through the 2nd order. A second order L/R lowpass is simply the cascade of two first order Butterworth lowpass. The first order butterworths are
$$L = \frac{1}{1+jx}, H = \frac{jx}{1+jx}$$
where x is the normalized cutoff frequency $x = \omega/\omega_0$
That makes the difference of the L/R filters
$$H = L^2-H^2 = \frac{1^2}{(1+jx)^2} - \frac{(jx)^2}{(1+jx)^2} = \frac{1+x^2}{(1+jx)^2}$$
The numerator is real, so we just have to proof that the magnitude of the denominator is the same.
$$|(1+jx)^2| = |(1+jx)|^2 = (1+jx) \cdot (1-jx) = 1 + jx =ix + x^2 = 1+x^2$$
which indeed equals the numerator.
Another interesting class of "recombining to allpass" are odd order Butterworth filters. Here the lowpass and the highpass are 90 degrees out of phase so both the sum AND the difference of the filters will be create an allpass, although one has significantly more group delay than the other and the "optimum" choice of the sign also alternates with order.
Odd order Butterworth also have only -3dB of gain at the crossover whereas L/R have -6dB.


Conclusions of sampling around Nyquist Rate I'm trying to understand some results of playing around with sampling around a signal's Nyquist sampling rate. For my example, I'm sampling a $B=5\mathrm{Hz}$ wave over a 1 second period.

In the first image, I use a very high sampling rate of 1000Hz to emulate the exact signal (at least for graphing), $\sin(2 \pi 5 t)$. This sampling is subsequently used in the next graphs as a reference.
In the next, I use a rate just below Nyquist's, 9Hz. The resulting wave looks like a sine wave, but doesn't represent the highs and lows of the wave. I'm still tracking with why hitting the Nyquist rate is important from this image.
Next, I get confused with Nyquist: I sample at exactly $2B=10\mathrm{Hz}$, the lower bound. This produces nearly zeros in the data:

[ 0.00000000e+00  1.22464680e-16 -2.44929360e-16 -1.40896280e-15
-4.89858720e-16  6.12323400e-16  2.81792560e-15 -2.69546092e-15
-9.79717439e-16  1.10218212e-15]

(I'm not even sure if these are floating point errors or not). Yes, there is some oscillation, but it is very minor. It seems like we need to further increase the sample rate... Is the rate of $2B$ exclusive?
In the fourth image, I sample above at 11Hz, which above the Nyquist rate. However, this sampling looks very much like image #2 flipped about the x axis.  Again, there is some oscillation, but I don't see how we could reconstruct a 5 Hz wave from this, especially when it has nearly the same qualities as image #2.
Finally, I sample at $4B=10\mathrm{Hz}$ and at last, I can see the 5Hz wave I expect. Why did it take 4x to get here when the theory only states 2x?

I feel like the conclusions I've made above are incorrect (in probably a variety of ways). Would someone please help me explain where I'm wrong?

Jupyter Notebook available here: https://gist.github.com/t-mart/8a46e87938904c7f6a5a102ad6a2ef0e

Is the rate of 2B exclusive?

Yes. The sampling theorem states that the signal must be band limited to half the sample rate. That implies that the energy at the Nyquist frequency must be zero. In practice you need a healthy margin between the highest usable frequency and the Nyquist frequency. There is always some "transition band" that you need to get the energy from "useful" to "nothing"

  I don't see how we could reconstruct a 5 Hz wave from this

You can using Shannon Whittaker interpolation. https://en.wikipedia.org/wiki/Whittaker%E2%80%93Shannon_interpolation_formula It doesn't matter how it looks to your eye, it just matters that all the information is there, which it is.

  Why did it take 4x to get here when the theory only states 2x?

It doesn't take 4x. You need to use a quantitative mathematical criteria to determine whether a sampling is "good" or not. Visual inspection doesn't help here. The most straight forward would be: if you send this sequence into a D/A converter and compare the output with the original input, how different are the two signals? That's actually what the Shannon interpolation does: it mimics an DAC converter with an infinitely steep anti-aliasing filter.


How can I properly detect a ball moving in a film capture? I have problem in detecting a ball moving in a flume (image below).

since converting the frames into grayscale format didn't work I converted it into rgb and make it more saturated. problems are:
1- in the first frames which ball is in the right position, light disturbed the ball detection (image subtraction didn't work well)

So I did it without subtraction.

this is my problem in detecting the ball in the first frames. in the rest frames the problem is its shadow and the white dots on the ball. I can't see the exact ball to find it's exact position in each frame (its shadow lead to some uncertainties)

anyway, I'm using matlab. I shared the code below (subtraction is not working well because of light, note that ball moves since water make it move in the flume).

Any Ideas about how can I detect the ball more precisely? handle the light in the first frames? handle the ball's shadow?

thanks in advance.

% the code starts

I=imread('frame3860.png');
I=imadjust(I,stretchlim(I));
I=lensdistort(I,-0.3);
I2=imcrop(I, [263.5 85.5 1468 940]);
I2=imrotate(I2,-1,'bilinear','crop');
I2=rgb2hsv(I2);
I2(:,:,2)=I2(:,:,2)*1.9;
I3=I2(:,:,1)>0.57 & I2(:,:,1)<0.75;
I4=bwareaopen(I3,200);
I4=medfilt2(I4);
I4=imfill(I4,'holes');
str=strel('disk',1);
I5=imdilate(I4,str);
I5=I5-I4;
e = edge(I5, 'canny');
radii = 30:1:70;
h = circle_hough(e, radii, 'same', 'normalise');
[~,maxIndex] = max(h(:)); %find the global maximum of your accumulator array
[i,j,k] = ind2sub(size(h), maxIndex);
radius = radii(k);
center.x = j;
center.y = i;
N = 200;
theta=linspace(0,2*pi,N);
rho=ones(1,N)*radius;
[X,Y] = pol2cart(theta,rho);
cent(1,:)=[j,i]; % centroid of the ball in each frame

% code finished

If you have any idea for a better ball detection in my situation I will be  grateful.

P.S. I removed the images since I didn't have the reputation needed to post more than 2 images

thanks in advance

First, for removing shadows try to normilize your RGB image prior to converting it to HSV color model. For more information on this regard see :
http://aishack.in/tutorials/normalized-rgb/

In case the ball is always moving, Optical flow would be a good option I guess. Simply calculate optical flow and mark the region with highest flow as the ball.

The other option would be emplying a tracking algorithm, which, I feel MeanShift would be a good option here. However, a tracking algorithm requires you to pinpoint the object of the interest in first few frames. If the object is fixed (it's always the same ball) then tracking would diffinitly work.


DFT of a cosine/sine wave? I'm seeing everywhere that the DFT of a sinusoid (i.e. a sine that has been windowed with a rectangular window) are just two deltas. What I don't understand is why the DFT is not the convolution of the transform of the infinite sinusoid (two deltas) and the transform of the rectangular window (a sinc). Wouldn't this be the right thing to do?

If a pure sinusoid is exactly integer periodic in the DFT's aperture, then the Sinc function, resulting from the rectangular window, will pass exactly thru zero at all but two points of the DFT results, thus appear to be invisible except for the two deltas.  If you vary the DFT window size by fractions of a period of the sinusoid, then it will show up, as you expect, since the DFT will no longer sample the Sinc when it crosses zero.


Could someone please explain to me the relationship between frequency bands and channels? Specifically, I am a tech writer working with WiFi hardware, and I need to understand why the 2.4GHz frequency band has only three channels, while the 5.0 GHz band has so many more.

I would presume that, since they're both "frequency bands", they would both have about the same bandwidth, and so would have about the same number of channels.   Yet it appears that the 5.0GHz spectrum is divided into channels about 5MHz wide, and occupies the 4.9GHz to 5.9GHz spectrum, while the 2.4 GHz band has only three or four usable channels (1, 6, 11, and 14), and occupies only the 2.4~2.5GHz spectrum.

Were these allocations strictly arbitrary, or is there some sort of physical explanation for these allocations?  Does anyone know why the decision was made to allocate these frequency bands in the way they have been set up?

The answer is mostly political: spectrum allocation is handled by the regional authorities based on a large number of economical, technical, political factors, and not the least on how much money can me made selling or leasing the rights. In the US that's done by the FCC.

This being said, most bands tend to have constant relative bandwidth, not absolute bandwidth. Bands that are centered in the GHz range are typically 100s of MHz wide. Bands that are centered in the 100kHZ range are maybe 10s of kHZ wide (if at all).

Using this reasoning you would expect 5 GHz band to have roughly twice the bandwidths than the 2.4 GHz band. Since the WIFI channels are constant absolute bandwidth (20 MHz), there will be more channels the higher the center frequency of the band is. But again, it's most regulatory determined.


How to eliminate negative frequencies from IQ signal I have a 192kHz IQ signal from an RF receiver, and i'm trying to remove signals in the negative (or positive) frequency spectrum.

I see that the negative frequency signals are -90 degrees phase shifted from I, where positive frequencies have the usual +90 degrees shift I vs Q.

Here's the problem:

If i set the center frequency to f, and i have two signals, one at f+1000, one at f-1000, they would both be audible on the same frequency, if i mix the whole thing into audio, either interfering with each other, or cancelling each other out.

I'm missing the step where i can decide whether i want to hear anything below or above center.

How is this usually being filtered?

I managed to do this now by having a 3 step process after converting the I/Q data to complex numbers. This is far from me understanding what i'm doing, but it's eliminating negative frequencies.

I shift the the frequencies by multiplying each i/q sample with the complex output of an oscillator at frequency PI/2. The oscillator outputs a complex number with Real being sin(pi*0.5*t), and Imaginary being cos() of the same values. Then I multiply that complex number with each sample (t advances +1 with each sample).
Next step is to high-pass filter at a quarter sample rate. So if our sample rate is 192000, my high pass cutoff is at 48000.
Lastly, we shift the whole thing back by doing the same multiplication as in 1), except with PI*1.5.

The result is I/Q data that's got everything in the negative spectrum filtered out.

I read about Hilbert transforms etc etc, but got lost halfway through that. Maybe someone has an explanation of what i did or why this works? :)


Removal the DC Component of Gabor Filters I recently was studying gabor filters and came to know that they have a shortcoming that they do have a minimum but some reasonable dc component.
now my question is why is it necessary to remove the dc component (using log gabor).

If there are any other limitation of gabor filters then please do tell me.
thanks..

gabor filter formed a sine/cosine wave with a gaussian, where the center of frequency of the filter is specified by the frequency of sine-cosine and the band width of the filter is specified by the width of the gaussian (valid perfectly inside a certain range and outside it, attenuates the frequency). In other words, more than one gabor filter to cover the spectrum.
the gabor filter can be designed for a band-width of 1 octave maximum.

you can find more details about gabor and log-gabor


Filter results from force measurements I was provided with several force measurements. Precisely:
The vertical compression force exerted from the milling tool on a working peace. The plot of the force looks reasonable and of course it seems also reasonable to me that the higher the force the higher the noise (because of the milling).
I am looking for some criteria to filter the measures. Any suggestion about how I could filter the data?

Signal processing becomes more valuable and useful when applied into the problems of different fields such as mining, seismology, oceaonography, climatology etc (apart from the obvious applications in multimedia, communications and control).

There are devices used to test and measure strength, stress and fracture behaviour of materials that's to be encountered in nature or in man made structures so that a good understanding of the materials stress strain behaviour enables stronger, durable and safer constructions.

One method of obtaining such important information of materials comes from specific lab equipments that test a sample material by providing a controlled force and measure resulting displacements. And the material characteristics is analysed according to the recorded force and displacement measurements.

Therefore the method requires these measurements to reflect the true conditions as accurate as possible in order to make reliable and precise conclusions from them.

There are many reasons why the measurements can be imperfect however. The distortions must have been analysed considering the sample-device mechanical dynamic interactions and any electronic noise or device nonlinearities that can also interfere with the measurements.

Under these circumstances, consider the following imperfect measurement results of a given sample test run:

As can be seen, this is quite an unexpected result from seemingly such a stable (slowly moving) device. One would for example expect a more monotonic curve (though not necessarily a strictly monotonic curve but at least not so much oscillatory either...)

Looking in detail to the first 8000 samples gives the following plot:

A Fourier analysis at this point could be of great help, revealed by the following plot:

Now we are enlightened by the presence of those harmonic spikes that are present at about 8 Hz, 16 Hz and 24 Hz respectively. In addition, there is also a DC looking spike at the origin, but that's actually the signal of interest that we consider as the true measurment which is very slowly (hence DC looking) moving down from -0.5 mV to -1.5 mV on the previous figure.

To verify the presence of 8, 16 and 24 Hz oscillations we further zoom into the signal to observe them directly in time domain rather than indirectly provided by the Fourier analysis computation...

In this very clear cropped portion of length about 100 samples, from the sample 1460 to sample 1560, there are approximately 8 periods of a fundamental oscillation. (And it's also understood from its waveshape that this is not a pure sine wave but involves harmonics of it which create the unsymmetric up down appearence).

Given the sampling period of $Ts=0.01$ seconds approximateley, this makes the fundamental period of oscillations to be $T_0 = 0.01*100/8 = 0.125$ seconds and therefore the fundamental frequency of the oscillations is found about $F_0 = 1/T_0 = 1/0.125 = 8$ Hz, Just as predicted from the Fourier analysis plot.

So therefore we are convinced both directly from the time domain measurement and also indirectly from the frequency domain computation that there are three quasi-oscillatory components at about 8Hz, 16Hz and 24Hz in addition to the signal itself around DC frequencies.

So in this case most probably due to mechanical interaction of the force piston and the sample material, there are uncontrolled vibrations occuring, which are casuing the oscillatory measurements. Or it could as well be electronic noise picked up by some reasons undecidable from the measurements alone.

Now, the are a number of solutions which would enable us to get rid of the unwanted oscillations (provied that these oscillations are indeed unwanted !) One of the simplest of them is to use a lowpass filter with an appropriate cutoff frequency. In the following I'm using a simple FIR lowpass filter of cutoff frequency below 8 Hz so that those oscillaitons are removed, but the useful signal is retained. The resulting plot is the following:

Clearly (or not!) this is an improvement, but a number of extra trials must have been performed in order to decide on the best solutions. The number of taps used was quite high approaching 256. This would induce a group delay of 128 samples, hence care must be taken to interpret the exact timing of the resulting samples. But this is not at all a problem, as long as you can perform non-real time (offline , non-causal) filering and that you know precisely how much delay there will be in the output which is the case here.

A much shother filter could have been designed which would yield much less sample shifts. Also instead of a lowpass filter a harmonic notch filter could also be designed, the decision of which is guided by the useful signal charactheristics.


direct Frequency domain FIR filtering vs Overlap-add method I'm trying to do bandpass filtering of a EEG signal samples at 250Hz and benchmarking the following 4 methods of FIR filtering for different filter orders. The length of the signal is 15000 samples.

Time Domain Approach (convolution)
Direct frequency domain approach
Overlap-add
Overlap-save

For follow the procedure listed in the following Wiki pages for overlap-add and overlap-save
http://en.wikipedia.org/wiki/Overlap%E2%80%93add_method
http://en.wikipedia.org/wiki/Overlap%E2%80%93save_method

For frequency domain approach I do FFT for the entire length of the signal, multiply with the frequency response of filter H and then IFFT the result.

For overlap and add, I chose M(overlap) to be 1+length_of_response and L to be twice of M.

Although one would expect Overlap-add to outperform Direct Frequency domain approach, this is not what I see in my benchmark results(see attached figure).

Please help me out in understanding where I might be going wrong.

An FFT of twice the filter order is kind of short.  Take a look at the "Choice of FFT size" section from this article I wrote a while back.
Also your whole signal will fit into a reasonable size FFT, so I think the single FFT approach is fine in this case.


Hilbert transform of low pass and high pass signal Why is the Hilbert transform of the product of a low pass signal and a high pass signal equal to the product of the low pass signal itself with the Hilbert of the high pass signal?

If $s(t)$ is a baseband signal and $x(t)$ is highpass, then under certain conditions, the Hilbert transform of $s(t)x(t)$ is $s(t)\hat{x}(t)$, where $\hat{x}(t)$ is the Hilbert transform of $x(t)$. This is known as Bedrosian's identiy (or theorem).

The easiest case to consider is when the high-pass signal is a sinusoid.

The frequency response of a Hilber transformer is
$$H(j\omega)=\begin{cases}-j,\text{ if }\omega>0 \\+j,\text{ if }\omega<0.\end{cases}$$

Note that the Hilbert transform of $\cos(\omega_0t)$ is $\sin(\omega_0t)$.

Now consider a band-limited baseband signal $s(t)$ with Fourier transform $S(j\omega)$. The spectrum of $s(t)\cos(\omega_0t)$ is $$\frac12 S(\omega-\omega_0) + \frac12 S(\omega+\omega_0).$$ The Hilbert transform of $s(t)\cos(\omega_0t)$ is then $$\frac{-j}{2} S(\omega-\omega_0) + \frac{j}{2} S(\omega-\omega_0)=\frac{-j}{2}\left(S(\omega-\omega_0) - S(\omega+\omega_0)\right),$$ which hapens to be the spectrum of $s(t)\sin(\omega_0t)$.

In conclusion, the Hilbert transform of the baseband signal $s(t)$ times the passband signal $\cos(\omega_0t)$ is the original baseband signal times the Hilbert transform of the passband signal.

This idea applies whenever you're finding the Hilbert transform of a signal whose spectrum is divided into two distinct bands which do not overlap; in other words, one band consists of only positive frequencies, and the other band is completely negative frequencies. This means that the lowest frequency component of the high-pass signal has to be larger than the bandwidth of $s(t)$. You can find more about the Bedrosian identity here, and a couple of proofs here.


2D deconvolution of recorded electron beam data I'm currently working on a project that involves using an electron gun and it would be really nice to know the spot shape of the electrons coming out of the gun (the frequency of electrons at some x,y away from the center of the beam).  Our current assumption is that the beam's spot shape is somewhat gaussian.

We have a picoammeter hooked up to a detector that outputs that current corresponding to the number of electrons that hit the detector.  The detector is circular with some radius, R.

If I'm not mistaken, sweeping the beam and recording readings from the picoammeter will result in an matrix that is the convolution of the beam spot shape and a 2D step function that's 0 outside of R and non-zero inside R.  I have looked through a few previous questions about deconvolutions related to images with a gaussian blur and it seems like this is very similar except the function is a step function instead of a gaussian.

I have tried using Matlab's conv2, FFT2, and iFFT2 functions to get the deconvolution as this thread suggests, but it produces something looking like this:

The code that produced this is:

mu = [0 0];
Sigma = [1 0; 0 1];
x = -3:.2:3;
[X1,X2] = meshgrid(x,x);
F = mvnpdf([X1(:) X2(:)],mu,Sigma);
F = reshape(F,length(x),length(x));
for i=1:31
    for j=1:31
        if x1(i)^2+x2(i)^2 >= 4
            circ(i,j)=0;
        else
            circ(i,j)=1;
        end
    end
end
conv = conv2(F,circ,'same');
D = ifft2( fft2(conv) ./ fft2(circ) );
surf(x1,x2,D);

Is there something inherently wrong with trying to deconvolves a step function?  Also, I tried using the full convolution, but that seemed to produce an even weirder deconvolution...

Am I just going about this the wrong way? (I've included python because I would prefer to use python if a solution exists, but matlab is fine too)

Your method of deconvolving seems to work.  Below is some Gnu Octave (more or less Matlab) code that implements your deconvolution.

Feed it two files that represent your target and your beam and it will perform a convolution of the two, then use your method to reverse the convolution.
The result shows the original two images, the convolved image, and the recovered shape of the beam.
The attached picture shows the results from my test.

Top left is the target (50x50px)

Lower left is the beam (14x16px)

Top right is the convolution of the two (the map of your electron counts)

Lower right is the recovered beam shape

The padding is a little off, so there's some wrap-around in the recovered image.  The result is clear enough, though.

I would check the images that you used for beam and target, and also check the convolved image to see if it is what you expected.  Part of the problem might also be the "same" that you used in the convolution.  That returns an image the same size as one of the convolved images, which will cause you to lose a lot of information.  If I use "same" then I don't get a recovered shape at all.  The convolved image is larger than the larger of the two images, so you have to use "full" then pad the target image to do the deconvolution.

ETA:
I got so busy looking at the details of the deconvolution that I forgot to mention your thoughts on the effect of scanning the beam across the target.

I agree that this is equivalent to convolving the beam shape with the target shape, so that I think you are on the right track to determining the beam shape - if I weren't convinced of that I wouldn't have bothered to check the deconvolution.

There's a few things that you will need to keep an eye on when you go to actually do the scan to determine the beam shape:

Scan area - you will need to overscan the target.  Take the radius of the target as R and the maximum radius (estimated) of the beam as r, then you will need to scan a square of at least R*2+2*(r*2) - actually more, but that is the absolute minimum.  Since you don't know the beam shape (and therefore its true size,) it may be better to go overboard here rather than have to repeat the measurements
The image you use to do the deconvolution must match (as perfectly as possible) the true shape and size of your target.  The position of the target shape within the image isn't as important.
When you pad your target shape image, center the image in the padding.  As I've done it, the padding is all on the top and left side.  This puts the right and bottom edges of the target shape right at the edges of the image.  Rounding errors can (and do, as you can see) cause part of the image to wrap around.  Since you will be creating the images, you could just create your image with the appropriate border instead of padding it in code.

target=imread("target.png");
beam=imread("beam.png");
convolved=real(fftconv2(target,beam,"full"));
widthc=size(convolved,2);
heightc=size(convolved,1);
widtht=size(target,2);
heightt=size(target,1);
padw=int32(widthc-widtht);
padh=int32(heightc-heightt);
paddedtarget=padarray(target,[padh,padw],0,'pre');
recovered=real(ifft2(fft2(convolved)./fft2(paddedtarget)));
subplot(2,2,1)
imshow(single(target));
subplot(2,2,2)
imshow(single(convolved));
subplot(2,2,3)
imshow(single(beam));
subplot(2,2,4)
imshow(single(recovered));


What Is the Difference between Difference of Gaussian, Laplace of Gaussian, and Mexican Hat Wavelet? There are three techniques used in CV that seem very similar to each other, but with subtle differences:

Laplacian of Gaussian: $\nabla^2\left[g(x,y,t)\ast f(x,y)\right]$
Difference of Gaussians: $ \left[g_1(x,y,t)\ast f(x,y)\right] - \left[g_2(x,y,t)\ast f(x,y)\right]$
Convolution with Ricker wavelet: $\textrm{Ricker}(x,y,t)\ast f(x,y)$

As I understand it currently: DoG is an approximation of LoG. Both are used in blob detection, and both perform essentially as band-pass filters. Convolution with a Mexican Hat/Ricker wavelet seems to achieve very much the same effect.

I've applied all three techniques to a pulse signal (with requisite scaling to get the magnitudes similar) and the results are pretty darn close. In fact, LoG and Ricker look nearly identical. The only real difference I noticed is with DoG, I had 2 free parameters to tune ($\sigma_1$ and $\sigma_1$) vs 1 for LoG and Ricker. I also found the wavelet was the easiest/fastest, as it could be done with a single convolution (done via multiplication in Fourier space with FT of a kernel) vs 2 for DoG, and a convolution plus a Laplacian for LoG.

What are the comparative advantages/disadvantages of each technique?
Are there different use-cases where one outshines the other?

I also have the intuitive thought that on discrete samples, LoG and Ricker degenerate to the same operation, since $\nabla^2$ can be implemented as the kernel
$$\begin{bmatrix}-1,& 2,& -1\end{bmatrix}\quad\text{or}\quad\begin{bmatrix}
0 & -1 & 0 \\
-1 & 4 & -1 \\
0 & -1 & 0
\end{bmatrix}\quad\text{for 2D images}$$.

Applying that operation to a gaussian gives rise to the Ricker/Hat wavelet. Furthermore, since LoG and DoG are related to the heat diffusion equation, I reckon that I could get both to match with enough parameter fiddling.

(I'm still getting my feet wet with this stuff to feel free to correct/clarify any of this!)

Laplace of Gaussian

The Laplace of Gaussian (LoG) of image $f$ can be written as

$$
\nabla^2 (f * g) = f * \nabla^2 g
$$

with $g$ the Gaussian kernel and $*$ the convolution. That is, the Laplace of the image smoothed by a Gaussian kernel is identical to the image convolved with the Laplace of the Gaussian kernel. This convolution can be further expanded, in the 2D case, as

$$
f * \nabla^2 g = f * \left(\frac{\partial^2}{\partial x^2}g+\frac{\partial^2}{\partial y^2}g\right) = f * \frac{\partial^2}{\partial x^2}g + f * \frac{\partial^2}{\partial y^2}g
$$

Thus, it is possible to compute it as the addition of two convolutions of the input image with second derivatives of the Gaussian kernel (in 3D this is 3 convolutions, etc.). This is interesting because the Gaussian kernel is separable, as are its derivatives. That is,

$$
f(x,y) * g(x,y) = f(x,y) * \left( g(x) * g(y) \right) = \left( f(x,y) * g(x) \right) * g(y)
$$

meaning that instead of a 2D convolution, we can compute the same thing using two 1D convolutions. This saves a lot of computations. For the smallest thinkable Gaussian kernel you'd have 5 samples along each dimension. A 2D convolution requires 25 multiplications and additions, two 1D convolutions require 10. The larger the kernel, or the more dimensions in the image, the more significant these computational savings are.

Thus, the LoG can be computed using four 1D convolutions. The LoG kernel itself, though, is not separable.

There is an approximation where the image is first convolved with a Gaussian kernel and then $\nabla^2$ is implemented using finite differences, leading to the 3x3 kernel with -4 in the middle and 1 in its four edge neighbors.

The Ricker wavelet or Mexican hat operator are identical to the LoG, up to scaling and normalization.

Difference of Gaussians

The difference of Gaussians (DoG) of image $f$ can be written as

$$
f * g_{(1)} - f * g_{(2)} = f * (g_{(1)} - g_{(2)})
$$

So, just as with the LoG, the DoG can be seen as a single non-separable 2D convolution or the sum (difference in this case) of two separable convolutions. Seeing it this way, it looks like there is no computational advantage to using the DoG over the LoG. However, the DoG is a tunable band-pass filter, the LoG is not tunable in that same way, and should be seen as the derivative operator it is. The DoG also appears naturally in the scale-space setting, where the image is filtered at many scales (Gaussians with different sigmas), the difference between subsequent scales is a DoG.

There is an approximation to the DoG kernel that is separable, reducing computational cost by half, though that approximation is not isotropic, leading to rotational dependence of the filter.

I once showed (for myself) the equivalence of the LoG and DoG, for a DoG where the difference in sigma between the two Gaussian kernels is infinitesimally small (up to scaling). I don't have records of this, but it was not difficult to show.

Other forms of computing these filters

Laurent's answer mentions recursive filtering, and the OP mentions computation in the Fourier domain. These concepts apply to both the LoG and the DoG.

The Gaussian and its derivatives can be computed using a causal and anti-causal IIR filter. So all 1D convolutions mentioned above can be applied in constant time w.r.t. the sigma. Note that this is only efficient for larger sigmas.

Likewise, any convolution can be computed in the Fourier domain, so both the DoG and LoG 2D kernels can be transformed to the Fourier domain (or rather computed there) and applied by multiplication.

In conclusion

There are no significant differences in the computational complexity of these two approaches. I have yet to find a good reason to approximate the LoG using the DoG.


How to distinguish between received quadratures How can you distinguish between two quadratures in QAM when the phase is not known at the receiver? For example, if the quadratures are sin and cos, how can the receiver tell which one is which, since both look the same except for the phases?

I've found the idea of random-phase carrier but I'm not sure if this is related to this question, or how it can be used if it is.

You would need to recover some timing information from the transmitted signal.

A common way is to transmit some sort of preamble before the rest of the data in your frame that is transmitted. The preamble would contain known symbols, so it could be used for multiple purposes, including channel estimation and timing recovery (also known as synchronization).

Once you can recover timing, e.g., from using the received preamble, you can attempt to demodulate your QAM symbols. And besides the phase information that you recover, you also want to avoid any difference in carrier frequency of the signal coming in to the receiver, and the oscillator in the receiver. That's where carrier frequency offset estimation comes into the picture too.


Pattern recognition using 2D fft in matlab I have a dataset with two columns (width and length, both in milimeter). When I make a scatter plot (please refer to the first figure), I see location of defects in a plane. I want to find if the location of the defect has some kind of pattern or it is just a noise. For that I was suggested to use 2D fft to recognize pattern. I don't have any background in 2D fft but I read some articles and used matlab function fft2. I get a plot (please refer to the second figure) which I dont I understand so much at this moment. Can any one has hint to get me going? A nice article could also be of help. I read the following article as well.

http://www.qsimaging.com/ccd_noise_interpret_ffts.html

Thanks in advance!

DFT is a great tool to find periodic patterns in signals.  However, you should understand what kind of period you want to find. If you want to find a series of points that lie on a regular grid (with even spacing), or more mathematically speaking, a subset of points $(x_i,y_i)$, such that there exists $(dx,dy)$ for which the following equations hold:

$x_i = x_1 + dx*i $
   $ y_i =y_1 + dy*j $

(The patterns can also be rotated, and DFT will still find them).

You should compute a 2D histogram of your data, with large number of bins, and run DFT on it. Afterwards, you should look for a peak in any location expect the DC. An ideal noise should have no patters, thus the contrast of the frequencies responses should be very low.


Synchronize a software waveform generator I'm having some trouble with a software waveform generator I'm attempting to put together. Especially with synchronizing the waveform.

Goal

Generate a mains synchronized (50 Hz) sinusoidal waveform to output on the DAC of a STM32F103.

Current Method

I developed a NCO (numerical controller oscillator) based on the accumulator method from this wiki.
This works, but my first attempt did not have precise enough calculations for the accumulator increment to get a reasonable close to the original frequency, so I had to move to q7.24 for this.
My samplerate is 10 kHz, with an interrupt every 2 ms.

Problem

I'm measuring the mains frequency with a 4 MHz free-running timer, and calculate the frequency to provide a q15.16 variable. This provides the frequency with reasonable precision. It is satisfactory for all other products that require frequency measurement.
It is however, not accurate enough for this purpose. There will be error due to the deviation of the 12 MHz crystal.

This is causing a problem, since the signal will need to synchronize to this 50 Hz reference. Which it will need some feedback for.

I developed phase measurement using the square wave subtraction method. This can calculate the phase within the limits of the samplerate.
With this result I feed a PI control loop that adds or subtracts from the accumulator increment in the NCO.
This keeps it in phase, but it still drifts around significantly on changes of mains. (and it changes a lot)

Is there a better way to synchronize a software numerical controlled oscillator? Or can I improve my current method?

Code

// Calculate increment for ac waveform generator accumulator
uint64_t f = (FREQ()->freq) << 8;  // Convert frequency q15.16 to q7.24
uint64_t t = (f * 1678);           // Multiply by sampleperiod, f * (1/fs) = f * 0.0001
inc = t / (1<<24);                 // q24*q24=q48, need to shift back to q24
// 50 hz approx 83886 q24

// Generate AC Waveform
// Add to (f * (1/fs)) to accumulator
accumlator24 += inc;
accumlator24 += phase_correction;
// Simulate 24 bit unsigned integer
if(accumlator24 & 0xFF000000u){
  accumlator24 = accumlator24 & 0x00FFFFFFu;
}
// Convert accumulator from q24 to q16 for fixed point sine
uint32_t accumulator16 = accumlator24 / 256;
// Multiply q16 accumulator "0.0 to 1.0" with 2*PI
x = fix16_mul(accumulator, fix16_2pi);
x = fix16_cos(x);
// Fill phase detect buffer (bitbanded)
if( x > 0 ){
  WFG_phase_bb[WFG_phase_dac_i++] = 1;
}else{
  WFG_phase_bb[WFG_phase_dac_i++] = 0;
}
// Multiply process value
x = fix16_smul(x, amplitude);
// Move to DAC range (0.0-1.0)
x = fix16_sadd(x, F16C(1,0));
x = fix16_smul(x, F16C(0,5));
x = fix16_clamp(x, 0, 0xFFFF);
// Copy lsb of fix16 (0.0-1.0) to dac buffer
dac1_data[i] = x;

Solution

In my implementation above the phase measuring of the generated vs reference waveform was lacking.
Because the resolution of this is equal to the samplerate due to the reference being a binary signal. I have improved the phase measurement from the binary comparison to the multiplication method. This allows for the mentioned topologys by Dan Boschen to be used.

I've added input normalization. This is actually remarkable easy to do if the signal is clean. Multiply by (1/peak).

if( input > peak ){
    peak = input;
}else{
    peak *= 0.995; // decay the peak, simulates a rectifier
}
factor = 1 / peak;
normalized = input * factor.

This normalized signal can be multiplied with the generated signal. This gives the phase as offset the multiplied signal. A low pass filter is used to get a stable value.

The phase detection with the multiplication method is proven to be accurate and provides a stable phase control loop.

It is not clear to me how you are determining your phase error so there may be issues there, and it also looks like  you are adding the phase correction to the accumulator directly. Consider adding the corrections to the frequency control word instead which would be analogous to a traditional analog PLL such as that shown in the graphic below (with the NCO being the equivalent to a VCO); this would allow the NCO to be an additional integrator in the loop which is necessary to implement a desired second order type 2 PI loop with your use of a PI Loop filter.  If you feed the corrections to the phase accumulator directly, you are treating the NCO as a phase control device instead of a frequency control device.

A type 2 implementation will track a frequency offset with zero error, while a type 1 implementation (which is what would result if the control was fed to the phase accumulator directly) would have an offset that is proportional to the inverse of the loop gain.

For the phase detector implementation, consider using a simple multiplier followed by a low pass filter to remove the double frequency output. (This multiplier can be analog or digital-- in your case since you are dealing with an analog signal an analog multiplier may be simplest. If you can be ensured the duty cycle is 50%, an XOR gate is perfect for this purpose. Otherwise a passive mixer is often a good choice).   The average value will be proportional the phase error with a quadrature offset (so the Loop would therefore lock in quadrature to the input signal).  Care should be taken that the low pass filter simply removes the double frequency component and is not too tight in bandwidth so as to compete with the loop bandwidth (it is preferred to set the loop bandwidth with the loop filter alone). Depending on implementation (see diagrams below) a "Phase/Frequency Detector" may be a preferred choice for the phase detector given that it won't lock to alias frequency offsets. If a multiplier type phase detector is used, the VCO/NCO output must be constrained to a narrow frequency range to avoid such false locks.

Three variants of a PLL approach are shown below. Depending on your need for digital and analog signal types in the rest of your system and programmability of loop dynamics one may be preferred over the other.

For the last implementation shown (which I think would be simplest) there are plenty of integrated chip solutions out there to make this even simpler. For example this chip: https://www.analog.com/en/products/hmc984.html#product-overview incorporates a phase/frequency detector together with a high range divider. The divider is 14 bits but 18 bits are required to get to 240,000, so the approach in this case would be to "pre-scale" the 12 MHz VCXO output with a divide by 16 to create a 750KHz clock, which would then be divided by 15,000 with this chip. There are many other integrated solutions to consider if you went down this path.


Grayscale/Brightness Resolution vs Contrast A bit of noob question, but can't find a definite answer anywhere:

What is the difference between Gray-scale resolution of an image and the Image contrast.
Are they different terminology for the same thing?

Gray-scale resolution and image contrast are different.

Gray-scale resolution looks at how many gray scales are being used to present the image.

Image contrast, on the other hand, looks at the difference between the maximum and minimum values of the gray scales used to present an image.

For example, a low gray-scale resolution, low contrast image might only use #000000 and #010101 (#RRGGBB where RR = red value, GG = green value, BB = blue value).

A low gray-scale resolution, high contrast image might only use #000000 and #FFFFFF (and only those two gray-scale values).

A high gray-scale resolution, low contrast image might have a maximum gray scale value of #1F1F1F (5 bits) and a minimum of #000000, but use all gray scales in between.

A high gray-scale resolution, high contrast image might have a maximum gray scale value of #FFFFFF (8 bits) and a minimum of #000000 (and use all gray scales in between).


Taylor series expansion in mean shift tracking I am trying to understand the first order Taylor series expansion when deriving the math behind mean shift tracking. Given the Battacharya coefficient between the target descriptor and candidate descriptor $\rho[\textbf{p}(y),\textbf{q}]$, how does
$$\rho[\textbf{p}(y),\textbf{q}] \approx \frac{1}{2}\sum_m\sqrt{p_m(y_0)q_m} + \frac{1}{2}\sum_mp_m(y)\sqrt{\frac{q_m}{p_m(y_0)}}$$
I am trying to see how each term can be related to the taylor series expansion about a point $a$ given by $$f(x) \approx f(a) + f'(a)(x-a)$$
What does $f$, $x$ and $a$ correspond to in the mean shift problem ?
Would appreciate any tips and pointers regarding this. Many of the lectures online do not really explain the details.

For $y\approx y_0$ you have
$$\begin{align}\sqrt{p_m(y)q_m}&\approx\sqrt{p_m(y_0)q_m}+\frac{p_m(y)-p_m(y_0)}{2\sqrt{p_m(y_0)}}\sqrt{q_m}\\&=\frac12\sqrt{p_m(y_0)q_m}+\frac12 p_m(y)\sqrt{\frac{q_m}{p_m(y_0)}}\tag{1}\end{align}$$
So the function you're approximating is $f(x)=\sqrt{x}$, with $x=p_m(y)$, and you do this in the vicinity of the value $a=p_m(y_0)$.


Find the nyquist sampling rate? I have a signal $x(t)$ for which I want to find the Nyquist frequency :
$$ x(t) = \frac{\sin{\pi t/2}}{\pi t/2} \ast \sum^\infty_{n=-\infty}\delta(t-10n)$$

I am trying to solve this in the time domain like so :
$$ \begin{split} x(t) &= \frac{\sin{\pi t/2}}{\pi t/2} \ast \sum^\infty_{n=-\infty}\delta(t-10n)\\
 &=\frac{\sin{\pi t/2}}{\pi t/2} \ast[...\delta(t+20)+\delta(t+10)+\delta(t)+\delta(t-10)+\delta(t-20)+...]\\
&=\frac{\sin{\pi (t+20)/2}}{\pi t/2}
+\frac{\sin{\pi (t+10)/2}}{\pi t/2}
+\frac{\sin{\pi t/2}}{\pi t/2}
+\frac{\sin{\pi (t-10)/2}}{\pi t/2}
+\frac{\sin{\pi (t-20)/2}}{\pi t/2}
+...\end{split}$$

Now shifting will not change the frequency so :

$$\begin{split} \omega_m &= \pi/2\\
2\pi f_m &= \pi/2\\
f_m &= 1/4\\
F_{\text{Nyquist}} &= 2\times 1/4 = 0.5\end{split}$$

But when I solved this quesiton in the frequency domain I obtained $0.4$ which is the correct answer.

What is the mistake in this method?

The primary problem in your time domain analysis is the assumption that if both $x(t)$ and $x(t-d)$ have the same Nyquist sampling rate so will $x(t)+x(t-d)$.

You can see that the summation may alter the signal in such a way that resulting signal's bandwidth may be different than the individual ones.

As an example consider the composite signal $x(t) = x_a(t) + x_b(t)$ where $x_a(t)$ is a low-pass signal with individual Nyquist sampling rate of $w_a$ and $x_b(t) = \sin(w_0 t)$, is a high frequency sine wave with a Nyquist sampling rate of $2w_0$. Therefore both $x(t)$ and $x(t-d)$ will have the same Nyquist sampling rate of $2 w_0$.

Then we have $x(t-d) = x_a(t-d) + \sin(w_0(t-d)) = x_a(t-d) + \sin(w_0 t - w_0 d)$. Therefore the sum of those two composite signals, $x(t) + x(t-d)$, will be $$x(t)+x(t-d) = x_a(t) + x_a(t-d) + \sin(w_0 t) + \sin(w_0 t - w_0 d)$$

Now if $d$ is chosen such that $w_0 d = \pi$ then $\sin(w_0 t) + \sin(w_0 t - \pi) = 0$, for all $t$, therefore $x(t) + x(t-d) = x_a(t) + x_a(t-d)$. There will be cancellation of the high frequency terms so that the Nyquist sampling rate of the sum will be $w_a$ (assuming the sum of low-pass signals will not further alter it). Hence the sum will have a different Nyquist sampling rate eventhough $x(t)$ and $x(t-d)$ individually have the same Nyquist sampling rate.


Temperature invariant color space I'd like to build a color difference metric that will be invariant to small changes in lighting conditions. In other words, I should be able to tell that a red shirt is not the same as a blue shirt, even if one image is lit differently (brighter / darker, outdoors / incandescent).

Moving to YUV space and ignoring the luminance component I can deal with changing brightness conditions by considering the metric:
$$d = \sqrt{\Delta U^2  + \Delta V^2}$$
But this metric still gives large "distances" when different light sources are used, due to the change in chrominance along the Planckian locus.

My question is:

Is there a color difference metric that is invariant to temperature differences? in a sense, I'm looking for a color space in which one of the three axes is temperature. Does such a space exist?

Other than searching for a color space, why don't you just first apply a color correction and then get the hue space?

If illumination is an issue, you might consider an illumination invariant color correction:
http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.108.6859&rep=rep1&type=pdf

Then simple convert your image to HSV and get H channel.

Additionally, if you have some clue about what you are looking at (such as roads, football fields etc) then you might use that too. Such an algorithm is mentioned here:
https://stackoverflow.com/questions/10922467/illumination-invariant-image


Find state space model from transfer function Let's suppose we have:

G(s) = (s+1)/(s^2-2s+1)

how can we find the state space representation of the transfer function:

x_dot = x2

x2_dot = 2*x2-x1+u

where u is an arbitrary input.

I am very new to this topic, so a detailed answere would be great ! :)

So, assuming that the transfer function is between the output $Y(s)$ and the input $U(s)$, namely
$$
\dfrac{Y(s)}{U(s)} = G(s)\qquad\qquad(1)
$$
multiplying (1) by $U(s)\cdot(s^2-2s+1))$ yields
$$
s^2 Y(s) -2s Y(s) + Y(s) = sU(s) + U(s)
$$
Going back in the time domain we obtain
$$
\ddot y -2\dot y + y = \dot u + u\qquad (2)
$$
Now, we look for a realization of the kind
\begin{align*}
\dot x &= Ax+Bu\qquad (3)\\
y &=Cx
\end{align*}
with $x=(x_1,x_2)$ and
\begin{align*}
A &= \begin{pmatrix}0 & 1\\a_1 & a_2\end{pmatrix} & B&=\begin{pmatrix}0\\1\end{pmatrix}, & C&=\begin{pmatrix}c_1&c_2\end{pmatrix}\end{align*}
The next step is to find the values of $(a_1,a_2,c_1,c_2)$ for which (3) has the same input-output behaviour of (2). From (3) we have
\begin{align*}
 y&=c_1x_1+c_2x_2\\
\dot y&= c_1\dot x_1 + c_2\dot x_2 = c_1x_2 + c_2a_1x_1 + c_2a_2 x_2 + c_2u\\
\ddot y&= c_1\ddot x_1 + c_2\ddot x_2 = c_1a_1x_1+c_1a_2x_2 +c_1u+ c_2a_1x_2+ c_2a_2a_1x_1+c_2a_2^2 x_2 + c_2a_2u + c_2\dot u
\end{align*}
substituting into (2) yields
\begin{align*}
0&=(c_1+c_1a_1+c_2a_2a_1-2c_2a_1)x_1 + (c_2+c_1a_2+c_2a_1+c_2a_2^2-2c_1-2c_2a_2)x_2\\
&+(c_1+a_2c_2-2c_2-1)u + (c_2-1)\dot u
\end{align*}
since that equality must hold for all $(x_1,x_2,u,\dot u)$ that' equivalent to ask
\begin{align*}
c_1+c_1a_1+c_2a_2a_1-2c_2a_1&=0\\c_2+c_1a_2+c_2a_1+c_2a_2^2-2c_1-2c_2a_2&=0\\
c_1+a_2c_2-2c_2-1&=0\\c_2-1&=0
\end{align*}
and a solution is
$$
(a_1,a_2,c_1,c_2)=(-1,2,1,1)
$$
Therefore
\begin{align*}
\begin{pmatrix}\dot x_1\\\dot x_2\end{pmatrix} &=\begin{pmatrix}0&1\\-1&2\end{pmatrix}\begin{pmatrix}x_1\\x_2\end{pmatrix} + \begin{pmatrix}0\\1\end{pmatrix}\\y&=\begin{pmatrix}1&1\end{pmatrix}\begin{pmatrix}x_1\\x_2\end{pmatrix}
\end{align*}


Highest frequency component in 2D Fourier transform I was trying to understand the 2d dft(discrete Fourier transform) and discrete spatial Fourier transform (dsft) of images.
But while understanding the graph i didn't understand that how come ($\pi,\pi$) is high frequency component and
not ($2\pi,2\pi$) is the highest frequency component.

Note that the terms $e^{j\omega_in_i}$ in the formula of the Fourier transform are $2\pi$-periodic. Consequently, the Fourier transform $H(\omega_1,\omega_2)$ is also periodic in both dimensions with period $2\pi$. If we choose $\omega=0$ as the center of our fundamental interval (as shown in the figure), the fundamental interval extends from $-\pi$ to $\pi$. The frequency $\omega=\pi$ is referred to as the Nyquist frequency, and it corresponds to the maximum possible frequency of a discrete signal.


Plotting convoluted and original signals on the same graph So I have a frequency response (got it by windowed fft-ing) of a note E4 played on a guitar and now I'm trying to implement a dynamic threshold for finding the peaks (fundamental frequency of the note and its harmonics).

I found an article that describes doing it by convolving normalized Hamming window with the frequency response. Now, I've done it, but I just don't know how to plot it correctly as the convoluted signal is delayed for a few samples like this (the red one):

Here's the code in Scilab:

// plot the frequency response
x_axis = linspace(0, Fs, row);
plot(x_axis(1 : $/2), win_response);

// calculate convolution on each window, siggram_n is just normalized
// matrix of windowed freq. response
[rows_sig, cols] = size(siggram_n);
conv = [];
for i = 1 : cols
    conv = [conv (convol(siggram_n(:, i), win_response))'];
end

// calculate the average of all convolutions
[rows, cols] = size(ko);
conv_n = [];
for i = 1 : rows
    conv_n = [conv_n sum(conv(i, :))/cols];
end

// raise to the power of 0.7 to flatten the threshold - from the article
conv_n = conv_n .^ 0.7;

// plot the convoluted signal
x_axis = linspace(0, Fs, 2*length(conv_n));
plot(x_axis(1 : $/2), conv_n, 'red');

It's not only important that I know how to plot this, because I'll need to find the peaks, so I have to somehow "synchronize" both signals. What am I doing wrong?

Looks like you need to do a couple of things:

Use conv instead of convol
Change the order of the parameters of convol like this: conv( win_response,siggram_n(:, i))
Tell conv that it should return an array the same size as the first parameter conv( win_response,siggram_n(:, i),"same")

If you don't use "same" then you get a convolution result back that is the length of the first array plus twice the length of the second array.  Convol has the parameters in the opposite order, and doesn't have a "same" option.

Also, change the name of your variable conv to something else.  It is bad practice to use keywords as variable names.  In this case, you have to change it else you can't use the conv method properly.

See the documentation for conv


Laplacian Operator with and without Diagonal Direction Elements in the Kernel This is a general question on the laplacian operator, which has two different versions. The first version is :

\begin{matrix}
  0 & 1 & 0 \\
  1 & -4 & 1 \\
  0 & 1 & 0
 \end{matrix}

The second version includes the diagonal:

\begin{matrix}
  1 & 1 & 1 \\
  1 & -8 &1 \\
  1 & 1 & 1
 \end{matrix}
Including the diagonal direction will make the laplacian isotropic. But what exactly is isotropic? How would that affect the resulting image if the diagonal components are included? What happen if I don't include the diagonal components?

All of those are Discrete Approximation of the operator - Laplace Operator.
You can discretize it in any logical manner.
In the cases above, Istotropic means if you rotate it it looks the same.
The diagonal elements means the filter would be sensitive to changes which are in 45, 135, 225 and 315 degrees, it it means it will also amplify more noise.


Scaling of a rect function,drawing of signals Im trying to understand on how the drawing of signals (time shifting,scaling...) works. This one example is confusing me, and a bit of help would be much appreciated. This is the starting Signal

We're susposed to give an $x(t)$ expression using rect, dirac and I've done so correctly;

$$x(t) = (t+2) \operatorname{rect}(t/2) + 2\delta(t-2)$$

Now we're susposed to draw a signal $f(t) = x(t/2)$;

And the rect function that is from -1 to 1 only "moves"to 2 and -2, the dirac moves to 4 that makes sense to me, but why doesn't the rect function move to 3 and -3. Thank you

$f(t)$ is just time-stretched $x(t)$ by a factor of 2. Meaning $f(t)$ will take the same value at $t$, which $x(t)$ was taking at $\frac{t}{2}$. So, if the original $rect(t)$ was extending from $-1$ to $1$, then the time-stretched $rect$ will extend from $-2$ to $2$. It can not go from $-3$ to $3$. One thing to keep in mind is that all  time operations (shifting and scaling) will be applied on the independent variable $t$ in $x(t)$. So,
$$f(t) = x(\frac{t}{2}) = (\frac{t}{2} + 2).rect(\frac{t}{4}) + 2\delta(\frac{t}{2} - 2)$$


Fourier transform of $\cos(n\omega t)$ My question is probably very stupid, but I've been strugling for a while on it now...

In need to find the Fourier transform of $1+\cos^3(2\pi ft)$.

I wrote that :
$$\cos^3(2\pi ft)=\frac{\cos(6\pi ft)+3\cos(2\pi ft)}{4}$$
And so I have:
$$\delta(f) +\frac 18 \bigg[\delta(f-3f_0)+\delta(f-3f_0)\bigg] +\frac 38 \bigg[\delta(f-f_0)+\delta(f-f_0)\bigg]$$

So, on my spectrum, I should have a dirac at $0$, a smaller one at $f_0$ and a smaller at $3f_0$...

But when I process it with matlab (using fast fourier transform), I get this :

(With a frequency of $10\textrm{ kHz}$).

So the dirac I thought would be at $3f_0$ is in fact at $\frac{f_0}{2}$. What am I missing ?

Thank you so much everyone!

I am really not confident in my mathematics skills and so I was focusing on the Fourier transform, but as you all guessed, the problem was my sampling frequency way too small (I was using 25kHz). With 60kHz it works like a charm.
Still feels like an idiot though ;-)

Thank you !

PS : To sum up if anyone was facing the same kind of problem

The theoric part was wright. We have :
$$\delta(f) +\frac 18 \bigg[\delta(f-3f_0)+\delta(f-3f_0)\bigg] +\frac 38 \bigg[\delta(f-f_0)+\delta(f-f_0)\bigg]$$

The problem was with my simulation: I forgot to check if my sampling frequency was meeting the Nyquist-Shannon criteria ($Fs > 2F_{max}$).

So, with this code :

Fs = 7e4;            % Sampling frequency
T = 1/Fs;             % Sampling period
L = 100000;             % Length of signal
t = (0:L-1)*T;        % Time vector
fm = 10e3;

X = 1+(cos(2*pi*fm*t)).^3;
%X = 1 + 0.25*cos(3*2*pi*10000*t)+(3/4)*cos(2*pi*10000*t);
Y = abs(fft(X));

f = 0:Fs/L:Fs/2;
plot(f,Y(1:L/2+1))
xlabel('f (Hz)')
ylabel('|P1(f)|')

We get the expected result:


How does delta modulation work practically? How does delta modulation work in practice? And where is it used? I find it confusing that, in 1-bit quantization we use only 2 level; if we apply this quantization, there will be loss in information. So, why do we use this technique?

One key point to note about delta modulation schemes is that they quantize the change between the current value and the previously quantized value.

The other key point to note is that, in general, delta modulation schemes are run at sampling rates MUCH higher than the Nyquist rate would suggest. Early portable CD players used a scheme called MASH (multi-stage noise shaping) on their DACs to allow cheaper, more accurate reproduction of sound.

I can't find a reference quickly, but I seem to recall that for the 44.1kHz CD sampling rate, MASH 1-bit DACs ran at sampling rates in the order of 5-6 MHz.


Synchrosqueezing Wavelet Transform explanation? How does Synchrosqueezing Wavelet Transform work, intuitively? What does the "synchrosqueezed" part do, and how is it different from simply the (continuous) Wavelet Transform?

Synchrosqueezing is a powerful reassignment method. To grasp its mechanisms, we dissect the (continuous) Wavelet Transform, and how its pitfalls can be remedied. Physical and statistical interpretations are provided.
If unfamiliar with CWT, I recommend this tutorial. SSWT is implemented in MATLAB as wsst, and in Python, ssqueezepy. (-- All answer code)

Begin with CWT of a pure tone:

A straight line in the time-frequency (rather, time-scale) plane, for our fixed-frequency sinusoid over all time - fair. ... except is it a straight line? No, it's a band of lines, seemingly centered about some maximum, likely the "true scale". Zooming,

makes this more pronounced. Let's plot rows within this zoomed band, one by one:

and all superimposed, each for samples 0 to 127 (horizontal zoom):

Notice anything interesting? They all have the same frequency. It isn't particular to this sinusoid, but is how CWT works in correlating wavelets with signals.
It appears much of information "repeats"; there is redundancy. Can we take advantage of this? Well, if we just assume that all these adjacent bands actually stem from one and the same band, then we can merge them into one - and this, in a nutshell, is what synchrosqueezing does. Naturally it's more nuanced, but the underlying idea is that we sum components of the same instantaneous frequency to obtain a sharper, focused time-frequency representation.
Here's that same CWT, synchrosqueezed:

Now that is a straight line.

How's it work, exactly?
We have an idea, but how exactly is this mathematically formulated? Motivated by speaker identification and Empirical Mode Decomposition, SSWT builds upon the modulation model:
$$
f(t) = \sum_{k=1}^{K} A_k(t) \cos(\phi_k (t)), \tag{1}
$$
where $A_k(t)$ is the instantaneous amplitude and
$$
\omega_k(t) = \frac{d}{dt}(\phi_k(t)) \tag{2}
$$
the instantaneous frequency of component $k$, where we seek to find $K$ such "components" that sum to the original signal. More on this below, "MM vs FT".
At this stage, we only have the CWT, $W_f(a, b)$ (a=scale, b=timeshift); how do we extract $\omega$ from it? Revisit the zoomed pure tone plots; again, the $b$-dependence preserves the original harmonic oscillations at the correct frequency, regardless of $a$. This suggests we compute, for any $(a, b)$, the instantaneous frequency via
$$
\omega(a, b) = -j[W_f(a, b)]^{-1} \frac{\partial}{\partial b}W_f(a, b), \tag{3}
$$
where we've taken the log-derivative, $f' / f$. To see why, we can show that CWT of $f(t)=A_0 \cos (\omega_0 t)$ is:
$$
W_f(a, b) = \frac{A_0}{4 \pi} \sqrt{a} \overline{\hat{\psi}(a \omega_0)} e^{j b \omega_0} \tag{4}
$$
and thus partial-diffing w.r.t. $b$, we extract $\omega_0$, and the rest in (3) gets divided out. ("But what if $f$ is less nice?" - see caveats).
Finally, equipped with $\omega (a, b)$, we transfer the information from the $(a, b)$-plane to a $(\omega, b)$ plane:
$$
\boxed{ S_f (\omega_l, b) = \sum_{a_k\text{ such that } |\omega(a_k, b) - w_l| \leq \Delta \omega / 2} W_f (a_k, b) a_k^{-3/2}} \tag{5}
$$
with $w_l$ spaced apart by $\Delta w$, and $a^{-3/2}$ for normalization (see "Notes").
And that's about it. Essentially, take our CWT, and reassign it, intelligently.

So where are the "components"? -- Extracted from high-valued (ridge) curves in the SSWT plane; in the pure tone case, it's one line, and $K=1$. More examples; we select a part of the plane and invert over it as many times as needed.

Modulation Model vs Fourier Transform:
What's $(1)$ all about, and why not just use FT? Consider a pendulum oscillating with fixed period and constant damping, and its FT:
$$
s(t) = e^{-t} \cos (25t) u(t)\ \Leftrightarrow\ S(\omega) = \frac{1 + j\omega}{(1 + j\omega)^2 + 625}
$$

What does the Fourier Transform tell us? Infinitely many frequencies, but at least peaking at the pendulum's actual frequency. Is this a sensible physical description? Hardly (only in certain indirect senses); the problem is,  FT uses fixed-amplitude complex sinusoid frequencies as its building blocks (basis functions, or "bases"), whereas here we have a variable amplitude that cannot be easily represented by constant frequencies, so FT is forced to "compensate" with all these additional "frequencies".
This isn't limited to amplitude modulation; the less sinusoidal or non-periodic the function, the less meaningful its FT spectrum (though not always). Simple example: 1Hz triangle wave, multiple FT frequencies. Frequency-modulation suffers likewise; more intuition here.
These are the pitfalls the Modulation Model aims to address - by decoupling amplitude and frequency over time from the global signal, rather than assuming the same (and constant!) amplitude and frequency for all time.
Meanwhile, SSWT - perfection:

Is synchrosqueezing magic?
We seem to gain a lot by ssqueezing - an apparently perfect frequency resolution, violating Heisenberg's uncertainty, and partial noise cancellation ("Notes"). How can this be?
A prior. We assume $f(t)$ is well-captured by the $A_k(t) \cos(\phi_k (t))$ components, e.g. based on our knowledge of the underlying physical process. In fact we assume much more than that, shown bit later, but the idea is, this works well on a subset of all possible signals:

Indeed, there are many ways synchrosqueezing can go awry, and the more the input obeys SSWT's assumptions (which aren't too restrictive, and many signals naturally comply), the better the results.

What are SSWT's assumptions? (when will it fail?)
This is a topic of its own (which I may post on later), but briefly, the formulation's as follows. Firstly note that we must somehow restrict what $A(t)$ and $\psi(t)$ can be, else, for example, $A(t)$ can simply cancel out the cosine and become any other function. More precisely, the components are to be such that:

More info in ref 2.

How would it be implemented? There's now Python code, clean & commented. Regardless, worth noting:

For very small CWT coefficients, phase is unstable (just like for DFT), which we work around by zeroing all such coefficients below a given threshold.
For any frequency row/bin $w_l$ in SSWT plane, we reassign from $W_f(a, b)$ based on what's closest to $w_l$ according to $\omega (a, b)$, and for log-scaled CWT we use log-distance.

Summary:
SSWT is a time-frequency analysis tool. CWT extracts the time-frequency information, and synchrosqueezing intelligently reassigns it - providing a sparser, sharper, noise-robust, and partly denoised representation. The success of synchrosqueezing is based in and explanied by its prior; the more the input obeys assumptions, the better the results.

Notes & caveats:

What if $f$ isn't nice in $\omega(a,  b)$ example? Valid question; in practice, the more the function satisfies aforementioned assumptions, the less of a problem this is, as authors demonstrate through various lemmas.
In the SSWT of damped pendulum, I cheated a little by extending signal's time to $(-2, 6)$; this is only to prevent boundary effects, which is a CWT phenomenon that can be remedied; here's directly 0 to 6.
Partial noise cancellation? Indeed; see pg 536 of ref 1.
What's the $a^{-3/2}$ in $(5)$? Synchrosqueezing effectively inverts $W_f$ onto the reassigned plane, using one-integral iCWT.
"Fourier bad?" My earlier comparison is prone to criticism. To be  clear, FT is the most solid and general-purpose basis that we have for a signals framework. But it's not an all-purpose-best; depending on context, other constructions are more meaningful and more useful.

Where to learn more?
The refernced papers are a good source, so are MATLAB's wsst and cwt docs and ssqueezepy's source code. I may also write further Q&A's, which you can be notified of by subbing this thread.

References:

A Nonlinear Squeezing of the CWT Based on Auditory Nerve Models - I. Daubechies, S. Maes. Excellent origin paper with succinct intuitions.
Synchrosqueezed Wavelet Transforms: a tool for Empirical Mode Decomposition - I. Daubechies, J. Lu, H.T. Wu. Good followup paper with examples.
The Synchrosqueezing algorithm for time-varying spectral analysis: robustness properties and new paleoclimate applications - G. Thakur, E. Brevdo, et al. Further exploration of robustness properties and implementation details (including threshold-setting).


Wav file scale factor I seek a scaled audio file format for my dsp output.

I'd like to produce "scaled wav files", when I've done my pre processing. Eg. if the amplitude is 15.5 (of whatever unit), you should be able to open a standard wav/whatever file and see "15.5".

If you produce a normal 16 bit wav file, you have to tell the user, to scale all values with value xx.xx. This would be ok, if the scale factor could be included in the file, in a standardized fashion. But so far I haven't been able to find any "scale factor" field in the wav format?

Another way could be to produce a 32 bit float wav file. This could actually contain the value "15.5". Problem is that wav float values are defined to be in the range [ -1.0 : 1.0 ]. I've tried opening such a file in different freewares like eg. Audacity. But they don't take well to the larger scale. They still only display the [ -1.0 : 1.0 ].
Do you know a program that will display large scale float wav files?

(Btw I do realize that the "scaled audio" is not "audio" anymore.)

This link implies that .wav files have several standards for storing 32 bit floating point.  It looks like foramt 32-bit 16.8 float (type 1 - 32-bit) would do what you want.

As best I can see, SOX (reference) can work (read it at least) with this type of format, so you could look at the source code to see how to generate/read these files.  CoolEdit claims to be able read and write this format (see first link.)


How Do Steerable Pyramids Work? I'm reading a paper about texture synthesis using pyramids, and the very first section mention two methods for constructing pyramids, one is the Laplacian Pyramid (which I'm quite familiar with), the other one is "steerable pyramid". I got the idea of what the latter cope with in comparison with the former, but I do struggle to understand how they work (if I wanted to implement them from scratch for example).

Is there a sufficient simple explanation of how they work?

There's a diagram extrapolated from one of the papers I've been reading through, the intuition I have behind it is that I need to perform these operation in the fourier domain, but I might be wrong. Can anyone expound it?

Thank you

We give a precise description of the pyramid algorithm in this paper: http://www.ipol.im/pub/art/2014/79/

It's a wavelet transform.

If you like to use the exact filters you can use Fourier transform but the filters are usually well approximated in spatial domain on a small support (8 x 8 px^2). Note that the direct convolution is faster than the Fourier implementation when the filters have few pixels (less than 16 x 16 px).


Guitar fundamental frequency vs pitch When a note is played on the guitar, is the fundamental frequency always the strongest one or could a harmonic be the strongest one (in terms of amplitude). When could that happen? Also, when playing a note of, for example, 440Hz, is 220Hz also a harmonic or only larger multiples?

When I produce a spectrogram of a B3 guitar note I get this:

As you can see, the frequency of B3 (246Hz) does not seem to be the strongest frequency. Its approximate half (near 130) seems to be the strongest one.

This is the audio file (the B3 note is clear).

To generate the spectrogram, I use a STFT taken from the Librosa library:

import plot as plt

def transcribe(filename):
  y, sr = librosa.load(filename, sr=40000)

  D = librosa.stft(y)
  plt.plot_spectrogram(D)

The plot.py module:

import matplotlib.pylab as plt

def plot_spectrogram(stft_output, xlabel='Time', ylabel='Hz'):
  plt.figure()
  librosa.display.specshow(librosa.amplitude_to_db(stft_output,
    ref=np.max), y_axis='log', x_axis='time')
  plt.title('Power Spectrogram')
  plt.colorbar(format='%+2.0f dB')
  plt.tight_layout()
  plt.savefig('static/plots/' + 'spectrogram.png')

Pitch is psychoacoustic.  Thus, a harmonic frequency can be stronger than the fundamental, without changing the pitch heard.  This is common in some stringed instruments.  And a note's pitch frequency could even be missing from a notes frequency spectrum.  For instance: An old telephone circuit can cut off frequencies below 200 Hz, but people still have no problem recognizing a low male voice pitched at around 100 Hz.

For many timbres, humans will hear the harmonic spacing as the note pitch.  In your plot, the harmonic spacing appears to be around 130 Hz, so that will be the note pitch perceived (if played at the same sample rate at the plot axis suggest).  Not 260.


Increasing the number of points in the frequency spectrum I have an image with few pixels in length and height. For this image I calculated the two dimensional Fourier transformation. What I got for the frequency spectrum in one direction was a very frequency spectrum with very few points (no continuous first derivation).
To increase the number of points of the frequency spectrum, I used a window function, in my case a Hanning-filter/rectangular filter with zeros outside of the picture so that the filtered picture has much more pixels in length (horizontal) and height (vertical). This resulted in a frequency spectrum with more points so that the frequnecy spectrum in e. g. the horizontal direction looked more continuous but there were in contrast to the unfiltered image higher peaks for low frequencies. What I want to do now is to find a possibility to reduce that.
My approach was the convolution theorem (but I have no idea how to use it in two dimensions) which says, that $F(f) * F(g)$ = const $F(gf)$. I have $f$ (original image) and $g$ (multiplication by rectangular window: some $1$ in the centre, $0$ outside it) and in that way the fourier spectrum that I get is $F(gf)$. To get $F(f)$ I have to deconvolve $F(f)$ and $F(g)$.
But the problems are at first:

Holds the convolution theorem for this?
Will the matrix $F(f)$ that I will finally get, the same that I got by Fourier transformation without the window function of the rectangular but with more points in the frequency spectrum? In other words, will it reduce effectively the influence of the filter in the frequency spectrum?
Is there an effective way to calculate the deconvolution with e.g. matlab?

If there are problems, is there an alternative way to increase the number of points in the frequency spectrum without changing it at all? If not, how can I fit the points in the frequency spectrum in the best way?

I would be very grateful for any helpful hints.

You can also interpolate as many points as you need using a Sinc or windowed Sinc interpolation kernel.


what is the difference between buffering and queuing? I've some problems with definition of queuing and buffering so
what is the difference between buffering and queuing ?

I think that, especially in the context of  real-time DSP the terms are really talking about the same topic.

Perhaps the most general Wikipedia reference about buffers would include both first-in-first-out (FIFO) and last-in-first-out (LIFO) but a LIFO buffer is usually called a "stack".

If it's a FIFO buffer, we usually call it a "queue" and, particularly if input and output to the buffer are asynchronous and controlled by independent source and destination, then there is a corner of the science we call "queueing theory" that worries about the size of the FIFO buffer, where the input and output pointers are, and what needs to be done to avoid "buffer starvation" (i think i am not misusing the term).

As far as real-time DSP is concerned, with uniform sampling, the only thing related to this topic that I worry about is, simply double buffering, of which I can't find a good on-line reference right away.  By FIFO buffering both the input and output to a DSP process, one can write more efficient code where the processing block time is the time it takes to fill or empty the buffer (i.e. $B$ times the sample time, where $B$ is the number of samples in the buffer block).  At the beginning of that block time the process is guaranteed to have $B$ valid samples waiting for the input and the process need not guarantee the correctness of any of the output samples until the very end of the block period.  This can save processing overhead by loading states and coefficients and other parameters just once (rather than $B$ times, once for each sample) and storing the states just once.  So the cost of this overhead is amortized over the $B$ samples.  The price paid is a delay of $2B$ samples.


Alternative for pcshow() function in matlab I'm using Matlab R2015a and I want estimate depth
but an error occurs and I wonder that it doesn't know pcshow() function. Why it doesn't work and is there another alternative for that?

pcshow() was implemented in 2015b and in R2015a the equivalent function is called showPointCloud.


Separate key-related frequencies from audio signal I have various samples of various instruments sound, and I'm trying to classify them using pyAudioAnalysis and various Machine Learning algorithms. The samples record the instrument playing different notes, which obviously alters their frequency spectrum. I am looking for a way to separate the information about the note, ideally by reducing the harmonics of the relevant frequency, so that the feature extraction step does not give relevance to this information. Does it make sense? Is there a "simple" way to achieve this?

Thank you in advance.

Attenuating the harmonics of your instruments isn't really a viable option, because you'll end up annihilating the signal itself. The harmonics here are what are sampling the behaviour of your instrument.

I'm not sure if I'd call it simple, but a possible solution to separating the contribution of pitch in such a classification system might be cepstral liftering.

This is, calculate a cepstral representation of the recording of each instrument (the Fourier transform of the logarithm of the frequency response) and separate the periodic components (pitch) from the aperiodic (perhaps timbre?). Such separation is known as liftering - the analog of filtering in the cepstral domain.

The separated signals should roughly correspond to a periodic signal and an impulse response better characterising the instrument.

This is an approach used in to separate the harmonics of speech (carrying pitch information) from the transfer function of the vocal tract (which generates a filter forming vowels amongst other characteristics).

I'm not sure if this will work in your case, but I think it would be worth exploring.


Does lowpass filtering and decimating improve SNR? I did not manage to find clear answers to this question: does decimation after an appropriate lowpass filter improve the signal to noise ratio?

To investigate this question, I considered data generated by a sine plus some random noise. The FFT shows a very clear peak (two with the symmetry). If I naively downsample without filtering, the signal gets aliased and the SNR decreases, which makes sense: I removed relevant information about the sine while keeping the noise unchanged.
Then, I applied a lowpass filter to avoid aliasing when subsampling. The SNR did not seem to significantly change.

This question is of practical importance: when measuring a phenomenon, say at 1Hz, is it useful to use a very high sampling frequency such as 1kHz or more in order to improve the SNR by filtering and decimating (given that the memory and computational cost is not a problem, of course)?

So, let's not forget what SNR is: it's a relation of powers present.

The thing that improves SNR is a propoer low-pass: it leaves the signal power alone and reduces the power of the noise.

An ideal low-pass filter will leave zero noise outside its specified bandwidth – so it doesn't matter whether you "cut off" these bandwidths using decimation or not, from a pure SNR perspective.

  This question is of practical importance: when measuring a phenomenon, say at 1Hz, is it useful to use a very high sampling frequency such as 1kHz or more in order to improve the SNR by filtering and decimating (given that the memory and computational cost is not a problem, of course)?

Yes! That is called oversampling, and it's useful for a couple of reasons, foremost for doing so-called noise shaping, which allows you to get a lower noise floor in your band of interest. With appropriate filtering, you do get a lot of SNR gain.

This is actually done pretty often; Software Defined Radio devices often run at dozens to hundreds of Megasamples, no matter how much bandwidth you want to observe, and do the filtering and decimation in-device. That way, you get the oversampling benefits, and also, you don't waste money and potential phase linearity etc by implementing steep-edged analog anti-aliasing filters, but can get away with much more relaxed filter specs.


How can I reconstruct a Time series using it AR coefficients in MATLAB? I have estimated AR coefficients of a time series using "aryule" function in MATLAB. Now I want to obtain the error of the estimated model. I think at first I must reconstruct it. so How can I reconstruct a Time series using it AR coefficients in MATLAB?

If you drive your AR model with white noise, you can generate a time series that has the same power spectral density as the original time series.  But it won't be the same time series.


Why are analytic FFT and computational FFT not matching? I have an exponentially decaying function that I want to Fourier transform. My question is why I am not getting the same result with the analytical form and the computational one with np.fft?
I have a function r(t) defined as follow:

    import numpy as np
    fs =5e6
    f0 = fs/2/102
    dt = 1./fs
    t = np.arange(0,0.0021544,dt)
    tau = 30e-6
    r = lambda t : (1/tau)*np.exp(-t/(tau))

that I want to Fourier-transform. According to the definition of the Fourier transform of this function r I am expecting:

    tp =t[0:102] # so half of a period
    freq = np.fft.fftfreq(np.size(tp), d = dt)
    Rf_def = np.divide(1,1+1j*2*np.pi*freq*tau)  # This is the analytical definition of the fourier transform of my function r(tp)

But when I do :

    Rf = np.fft.fft(r(tp))

I don't get the same result as the analytical FFT that I defined above. Why is that ?

This is a plot of the analytical form of the fourier transform so plt.plot(Rf_def).

This is a plot of the compututional fourier transfrom so plt.plot(Rf).

The fourier transform is defined for both positive and negative frequencies.  It looks like your first plot is just the positive frequencies, whereas the second plot contains both the positive and negative frequencies.  It's  common in computational fft routines to have the positive frequencies occupy elements 1 to N/2 in the results vector, and then elements N/2 to N are the negative frequencies.  Also, computational results are often unscaled.  You'll need to divide the results by the length of the window. Here is a short article I wrote explaining this for matlab http://mechanicalvibration.com/Making_matlab_s_fft_functio.html  Your results for python should be similar.


What Is the Best First Order IIR (AR Filter) Approximation to a Moving Average Filter (FIR Filter)? Assume the following first order IIR Filter:

$$ y[n] = \alpha x[n] + (1 - \alpha) y[n - 1] $$

How can I choose the parameter $ \alpha $ s.t. the IIR approximates as good as possible the FIR which is the arithmetic mean of the last $ k $ samples:

$$ z[n] = \frac{1}{k}x[n] + \frac{1}{k}x[n-1] + \ldots + \frac{1}{k}x[n-k+1] $$

Where $ n \in [k, \infty) $, meaning the input for the IIR might be longer than $ k $ and yet I'd like to have the best approximation of the mean of last $ k $ inputs.

I know the IIR has infinite impulse response, hence I'm looking for the best approximation.
I'd be happy for analytic solution whether it is for $ {L}_{2} $ or $ {L}_{1} $ cost function.

How could this optimization problems can be solved given only 1st order IIR.

Thanks.

OK, let's try to derive the best:
$$
\begin{array}{lcl}
y[n] &=& \alpha x[n] + (1 - \alpha) y[n - 1] \\
&=&  \alpha x[n] + (1 - \alpha) \alpha x[n-1] + (1 - \alpha)^2 y[n - 2]\\
&=&  \alpha x[n] + (1 - \alpha) \alpha x[n-1] + (1 - \alpha)^2 \alpha x[n-2] + (1 - \alpha)^3 y[n - 3]\\
\end{array}
$$
so that the coefficient of $x[n-m]$ is $\alpha(1-\alpha)^m$.

The best mean-square approximation will minimize:
$$
\begin{array}{lcl}
J(\alpha) &=& \sum_{m=0}^{k-1} (\alpha(1-\alpha)^m - \frac{1}{k})^2 + \sum_{m=k}^\infty \alpha^2(1-\alpha)^{2m}\\
&=& \sum_{m=0}^{k-1} \left(\alpha^2(1-\alpha)^{2m} - \frac{2}{k}\alpha(1-\alpha)^m + \frac{1}{k^2}\right) + \alpha^2 (1-\alpha)^{2k} \sum_{m=0}^\infty (1-\alpha)^{2m} \\
&=& \alpha^2\frac{1- (1-\alpha)^{2k}}{1 - (1-\alpha)^2} + \frac{2\alpha}{k} \frac{1 - (1 - \alpha)^k}{1 - (1 - \alpha)} + \frac{\alpha^2(1-\alpha)^{2k}}{1 - (1 - \alpha)^2}+ \frac{1}{k}\\
&=& \frac{\alpha^2}{1 - (1 - \alpha)^2} + \frac{2}{k} (1-(1-\alpha)^k) + \frac{1}{k}\\
&=& \frac{\alpha^2}{2\alpha - \alpha^2 }+ \frac{2}{k} (1-(1-\alpha)^k)+ \frac{1}{k}\\
&=& \frac{\alpha}{2 - \alpha }+ \frac{2}{k} (1-(1-\alpha)^k)+ \frac{1}{k}\\
\end{array}
$$
because the FIR coefficients are zero for $m > k - 1$.

Next step is to take derivatives and equate to zero.

Looking at a plot of the derived $J$ for $K = 1000$ and $\alpha$ from 0 to 1, it looks like the problem (as I've set it up) is ill-posed, because the best answer is $\alpha = 0$.

I think there's a mistake here.
The way it should be according to my calculations is:

$$
\begin{array}{lcl}
J(\alpha) &=& \sum_{m=0}^{k-1} (\alpha(1-\alpha)^m - \frac{1}{k})^2 + \sum_{m=k}^\infty \alpha^2(1-\alpha)^{2m} \\
&=& \sum_{m=0}^{k-1} \left(\alpha^2(1-\alpha)^{2m} - \frac{2}{k}\alpha(1-\alpha)^m + \frac{1}{k^2}\right) + \alpha^2 (1-\alpha)^{2k} \sum_{m=0}^\infty (1-\alpha)^{2m} \\
&=& \alpha^2\frac{1- (1-\alpha)^{2k}}{1 - (1-\alpha)^2} - \frac{2\alpha}{k} \frac{1 - (1 - \alpha)^k}{1 - (1 - \alpha)} + \frac{1}{k} + \frac{\alpha^2(1-\alpha)^{2k}}{1 - (1 - \alpha)^2}
\end{array}
$$

Simplifying it according to Mathematica yields:

$$ J(\alpha) = \frac{\alpha}{2 - \alpha} + \frac{2 {(1 - \alpha)}^{k} -1}{k} $$

Using the following code on MATLAB yields something equivalent though different:

syms a k;

expr1 = (a ^ 2) * ((1 - ((1 - a) ^ (2 * k))) / (1 - ((1 - a) ^ 2)));
expr2 = ((2 * a) / k) * ((1 - ((1 - a) ^ (k))) / (1 - (1 - a)));
expr3 = (1 / k);
expr4 = ((a ^ 2) * ((1 - a) ^ (2 * k))) / (1 - ((1 - a) ^ (2)));

simpExpr = simplify(expr1 - expr2 + expr3 + expr4);

$$ J(\alpha) = \frac{-2}{\alpha - 2} - \frac{k - 2{(1 - \alpha)}^{k} + 1}{k} $$

Anyhow, those functions do have minimum.

So let's assume that we really only care about the approximation over the support (length) of the FIR filter. In that case, the optimization problem is just:
$$
J_2(\alpha) = \sum_{m=0}^{k-1} (\alpha(1-\alpha)^m - \frac{1}{k})^2
$$

Plotting $J_2(\alpha)$ for various values of $K$ versus $\alpha$ results in the date in the plots and table below.

  For $K$ = 8.     $\alpha_{\tt min}$ = 0.1533333
     For $K$ = 16.    $\alpha_{\tt min}$ = 0.08
     For $K$ = 24.    $\alpha_{\tt min}$ = 0.0533333
     For $K$ = 32.    $\alpha_{\tt min}$ = 0.04
     For $K$ = 40.    $\alpha_{\tt min}$ = 0.0333333
     For $K$ = 48.    $\alpha_{\tt min}$ = 0.0266667
     For $K$ = 56.    $\alpha_{\tt min}$ = 0.0233333
     For $K$ = 64.    $\alpha_{\tt min}$ = 0.02
     For $K$ = 72.    $\alpha_{\tt min}$ = 0.0166667

The red dashed lines are $1/K$ and the green lines are $\alpha_{\tt min}$, the value of $\alpha$ that minimizes $J_2(\alpha)$ (chosen from $\tt alpha = [0:.01:1]/3;$).


Finding the discrete time fourier coefficients to this problem I'm trying to find the fourier series to this discrete time signal.
$$x_1[n] =\begin{cases}
+\frac72&\text{if }0\le n \le 4\\
-\frac72&\text{if }5\le n \le 9
\end{cases}$$
My approach:
We see that the discrete period, $N$, is 10, so $\Omega_0 = \frac{2\pi}{10} = \frac{\pi}{5}$
Now using the formula to get the DT fourier coefficients, $X[k]$
\begin{align}
X[k] &= \frac{1}{N}\sum_{n = 0}^{n = 9}x_1[n]e^{-j\frac{\pi}{5}kn} \tag{1}\\
&= \frac{1}{10}\left( \sum_{n = 0}^{n = 4}\frac{7}{2}e^{-j\frac{\pi}{5}kn} + \sum_{n = 5}^{n = 9}\frac{-7}{2}e^{-j\frac{\pi}{5}kn} \right)\tag{2}\\
&= \frac{7j}{10}\sum_{n = 0}^{n = 4}\sin\left(\frac{\pi k n}{5}\right)\tag{3}
\end{align}
However according to the solutions for $k = 2$, we have $X[k = 2] = 0$, but when I substitute $kk = 2$ into the final form of the above expression, I don't get $0$. Thus, my answer must be wrong, but I'm not sure where I went wrong... Can someone please explain where my math is flawed and if there's an easier way to find the coefficients?


Smoothing process of Kalman filter I have a question about the smoothing (backward) process of Kalman filter. Is it correct to say

$E[x_{t|T}] = x_{t|t}$

where $x_{t|t}$ is the estimated result from forward process?

I am struggling to understand the smoothing process of Kalman filter. Any help will be appreciated.

The filtering distribution $p(x_k|y_1,...,y_k)$ is computed via the Kalman filter and given a linear-time invariant system is parametrized by the mean $E[x_k|y_1,...,y_k]$ and covariance matrix.

When doing smoothing we want to compute $p(x_k|y_1,...,y_N)$ with $N \ge k$, in this case the mean we compute is $E[x_k|y_1,...,y_N]$.


How can slope of low-pass biquad be configured? I have a biquad filter, configuread as a low-pass filter (using Zolzer's formulas, similar to these). It does not overshoot and decays nicely in what I assume is -40dB/dec.

How can I device a calculation, that takes not only desired_cutoff_frequency as a parameter, but also desired_decay_slope?

The decay is determined by the filter order. With a biquad (i.e., second order) filter you can't get more than approximately 40dB/decade (unless you're close to the Nyquist frequency where the frequency response has its zero). You can have less if you like, but not more. If you want less, then the only way to do that is to approximate a given slope by using some optimization method.