Formula for Bessel low-pass filter coefficients When I am filtering a signal in python, it has a built in function to generate bessel filter coefficients given a cutoff ratio and a filter order (number of poles). I am trying to translate this to C code, but I cannot seem to find the formula that is used to calculate the coefficients.

Can someone point me to the formulas to calculate the IIR filter coefficients given a cutoff ratio and order?

I am trying to read the scipy source on github but I am having a very hard time of it...

Hi, I wrote this.  It's probably more complicated/inefficient than it needs to be. :D

Practical Bessel filter design involves root-finding of a polynomial to generate second-order sections; I don't believe there is any simple formula.  The way SciPy stores filter coefficients is pole and zero locations, so the code finds these locations numerically.

The steps I originally used are:

Generate the coefficients of a reverse Bessel polynomial
Find the roots of the polynomial (= poles of the Bessel filter) using the Aberth-Ehrlich method
Scale the poles inward or outward to normalize either the -3 dB frequency, phase, or group delay.

Step 2 generates approximate roots from Campos-Calderón 2011 as the starting points for the root-finding, but I don't know if that's really necessary.  Most root-finding algorithms just use asymmetrical starting points on a spiral, etc. and still find the answer fine.  I assumed it would be faster this way.

Aberth-Ehrlich is good for quickly finding multiple single roots simultaneously, by modelling them as point charges that repel each other, so multiple test points don't fall into the same holes.

Eventually I ended up eliminating Step 1 altogether.  For the Newton's method part of Aberth-Ehrlich, it actually evaluates $K_{n+\frac 1 2}\left(\frac 1 x \right)$  instead, which is a totally different (compiled) function that has the same roots, and has a relatively simple derivative.  This is possibly crazy and less efficient than it could be, but it works.  (kve is actually zbesk from AMOS Fortran library.)

$\theta_{5}(x)$ vs $K_{5.5}\left(\frac 1 x \right)$:

Then in step 3, I don't even use the Bessel polynomial to normalize the function, since I eliminated the need for all but the last coefficient and it's faster to just generate that one by itself.

There's also this algorithm I found later, which could probably be made more efficient: Orchard 1965 - The Roots of the Maximally Flat-Delay Polynomials (Bessel filters) though it has some numerical error at high N and doesn't give an algorithm for finding the new approximate roots from the old (just "eyeball it on graph paper").

But now you can generate 500th order Bessel filters with 3 different normalizations, instead of just reading the coefficients from a table (not that you would ever need to...):


Derivative of equation containing Gaussian filtered image Let $w$ be our image. For example, consider the following with the vectorized $w$:
$$ E(w) = \frac 12 \|Aw+b\|_2^2 $$
I know the optimal condition of the above equation:
$$ \nabla_wE = A^TA\hat w + b = 0 $$

Now, we consider the filtered image by Gaussian such that
$$ E(w) = \frac 12 \|A (G * w)+b\|_2^2 $$

In this case, how to get the optimal condition? (i.e., how to derivative w.r.t $w$?)

You can represent the convolution as a matrix multiplication with a $N\times N$ Toeplitz matrix, where $N$ is the length of $\mathbf{w}$. Check out this explanation.
So you can write
$$
E(\mathbf{w})=\frac{1}{2}\|\mathbf{AHw}+\mathbf{b}\|^2_2
$$
where $\mathbf{H}$ is your Teoplitz matrix filter. Now if you define
$$
\tilde{\mathbf{A}} = \mathbf{AH}
$$
then
$$
E(\mathbf{w})=\frac{1}{2}\|\tilde{\mathbf{A}}\mathbf{w}+\mathbf{b}\|^2_2
$$

And you can see that the solution is exactly as you wrote but with $\tilde{\mathbf{A}}$ instead of $\mathbf{A}$.


How to difference a loud sound (playing notes), from a soft noise in spectrum. (Pitch detection with bass) I'm trying a program that outputs the bass line that you are playing.
I have some arrays from a big one with the sound in time domain, that are splited (Hann window).
Later, I do the FFT, then I get the peak of each array and apply the harmonic product spectrum algorithm.

Finally, I'm trying to filter all frequencies bellow 39hz and above 440hz, that include all frec. in the bass in E standart.
I've tried to eliminate the noise.

I want to distinguish my FFT output when I play my instrument, from when I don't, taking advantage of the frequency filter, I have made the program assume that when it plays a note that is less than 39hz, it automatically takes it as noise.
The problem appears when the noise is greater than 39hz (it can be from 0hz to 60hz even) and many of the notes I play on the bass are below 60hz, so with this system I cannot know when I make sound because I'am playing, and when the predominant sound is noise and I'm not playing.

So I need to know how to interpret the FFT output to know when I actually play and when I don't.

The FFT output is the energy (I think) of each frequency in the sound wave, and I have tried to eliminate all FFT outputs below 0.5 (The output is normalized from -1 as minimum value to 1 as maximum value of amplitude in the spectrogram), so I really don't know if this works.

Could anyone to tell me how to identify when I play and when I don't?

The intuitive answer is to threshold with respect to the magnitude.

You are performing a short-time Fourier transform (STFT), that is, segmenting the original into short time windows and performing Fourier transform on each of them. The output of STFT is a frequency-time (t-f) bin table of complex values.

There is an inherent tradeoff in STFT. The shorter your windows are the higher is your time resolution and the lower is your frequency resolution. If you want a high-frequency resolution you will have to use longer windows at the cost of poor time resolution. This paper gives a clue on window length optimization. otherwise, you can simply use the common music window length from papers or online examples. For bass, the frequencies are low = long wavelength, hence, you will want longer windows.

Basically, you can get the magnitude and\or the angle of your t-f table. Assuming your bass has much more energy than the noise gathered by the system, you can look at the magnitude of bins in the required frequency range and drop time segments below a certain threshold. Play a bit with the result and present it as an image. It will give you a science of where to cut. This is a MATLAB example of presenting the magnitude as an image (also called spectrogram)


Wireless Communication Let's assume I send a root raised cosine pulse of 1ns at t=0 through a wireless channel.
The power measured @ Receiver from t=0 is shown in image.
Now from this image, I know that there are multipath components(due to spreading of signal power across time).

I want to know that up to which time do I have to measure the power for a single symbol transmitted and how can I know that after a certain amount of time the second symbol will start. (how can I measure that certain amount of time?)
In this received power I know that(by spreading) that there are multipath components. I want to measure the Gain of those multipath components. So how can I sample it and how do I know about the number of multipath components?

Typically a known sequence is transmitted (often as a preamble) that has properties to aid in determining the channel characteristics (by spreading its energy across all frequencies within your signal bandwidth) and in determining the symbol clock (many transitions at the symbol clock boundary). Both can be accomplished by using pseudo-random sequences which also have the desirable property of an autocorrelation that approximates an impulse function (useful for establishing a timing boundaries under low SNR conditions). When shorter sequences are sufficient for  acquisition and channel estimation, Barker codes also have the desirable autocorrelation and signal energy spreading properties mentioned (these go hand in hand: the autocorrelation of white noise is an impulse).

The details on estimating the channel (or equalizing your signal to remove the channel distortions) for single-carrrier modulation schemes are identical to the operations detailed in this post here:  Compensating Loudspeaker frequency response in an audio signal


Adding White Gaussian noise to a voice signal I'm trying to add White Gaussian Noise to an audio file. However; the energy of the noise should be 1/10th of that of the signal. My first attempt is as following:

[y,Fs] = audioread('drum.wav');
%sound(y,Fs);
noisy_sig = addnoise(y,10);
sound(noisy_sig,Fs);

function noisy_sig = addnoise( sig , SNRdb )
   sig_power = norm(sig,2) / length(sig);
   % noise power is equal to sigma^2
   sigma2 = sig_power / 10^(SNRdb/10) ;
   noisy_sig = sig + sqrt(sigma2)*randn(size(sig));
end

When I listen the resulting wav, there is no difference between the original and noisy ones. Am I doing something wrong? Any help would be appreciated.

Your problem is how you calculate the signal power. You are calling norm which calculates $\sqrt{\sum_i |x_i|^2}$ but you want to calculate the sum of the squared values and then divide by the length as you do in your code. Try this instead:
sig_power = norm(sig)^2 / length(sig);. This line will calculate $\frac{1}{N}\sum_{i=1}^N |x_i|^2$.


Tx power vs Error Vector Magnitude I cannot tell how decrease/increase in Tx output power should impact EVM?
What is the mathematical relation of I,0 values on the constellation map and output power?

in the context of 802.11, I am trying to understand what is the limiting factor of the output power of a transmitter? (besides power consumption)

The Protocol does not specify output power number in dBm. It specifies EVM threshold for the different rates.

You need to get the transfer function for AM/AM and AM/PM conversion for your PA to get the mathematical relationship describing the gain and phase shift for the PA as a function of input amplitude (of the envelope of the signal specifically). This describes the non-linear behavior of the PA as you increase output power toward saturation; for low power levels the PA will be linear such that an X increase (in dB) or input power will correspond to a X increase in output power, or in phase there will be no additional phase shift that is dependent on the input power level. However at increased power levels the output power change will be less than the input (AM/AM) and additional phase shift is introduced (AM/PM).

From that it should be very clear how the error vector magnitude can be modified from such a non-linearity, especially with waveforms that are not constant envelope (implying that the input amplitude is changing and this creating the AM/AM and AM/PM distortions).

Thermal heating can be a factor, such that the static AM/AM and AM/PM translations will only approximate the true effects due to the associated memory effects.  My peripheral understanding is that more advanced modeling with Markov processes can be done to more accurately predict the behavior in these cases, particularly with modeling spectral regrowth.


Reconstruction of sampled band-pass signal I am pretty new to signal processing. I am currently trying to reconstruct a sampled band-pass signal created with the filtfilt and cheby2 MATLAB functions. I am trying to do this with sinc interpolation, but I do not know what the discrete reconstruction formula is. I know that

$$h(t) = f_s \cdot \text s \text i \text n \text c(\frac{t \cdot f_s}{2}) \cdot \text c \text o \text s (2\pi t f_s \frac{2n + 1}{4})$$

but I cannot be able to find the corresponding discrete implementation of this - i.e. the equivalent of $$\sum_{n = -\infty}^{\infty} x_n \cdot \text s \text i \text n \text c (t\cdot f_s - n)$$ for the low-pass filter reconstruction that has $$h(t) = f_s \cdot \text s \text i \text n \text c (t\cdot f_s)$$

Can sombody shed some light on this? I have been struggling for a lot of time with this issue (trying different implementations that didn't quite give the right result), but I think the solution should be rather straightforward.

You want to take a signal with sample rate 30 Hz and upsample the signal to 3 kHz rate and then pass it through a low pass filter to perform the sinc interpolation. In MATLAB, you have a couple options:

You could use the built-in interp function. You will need to pass it your signal and the upsample factor (3,000 / 30 = 100). The function will figure out a low pass filter to use and do everything for you. You can also get the filter it used (read the documentation to see how).
You could do it by hand by calling upsample, then design the filter by using the designfilt function for example.

Lets compare the two different methods:

fs1 = 30;
fs2 = 3000;
L = fs2/fs1; % upsample factor
N = 1000; % number of samples in original signal
t = [0:N-1].'/fs1;
sig1 = cos(2*pi*5*t);
upsampledSig1 = L*upsample(sig1, L);
lpFilt = designfilt('lowpassfir', 'PassbandFrequency', 0.25/L, 'StopbandFrequency', 0.5/L, 'PassbandRipple', 0.1, 'StopbandAttenuation', 95, 'DesignMethod', 'kaiserwin');
filterDelay = floor(length(lpFilt.Coefficients)/2);

Now lets use the MATLAB interp function and compare it against our "hand-made" method.

interpSig = interp(sig1, L);
paddedSig1 = [sig1; zeros(filterDelay, 1)];
mySig = filter(lpFilt, paddedSig1);
mySig = mySig(filterDelay+1:end); % strip off filter delay

As you can see, the results are similar with the only difference being a slight scaling difference which can easily be taken care of.


Discrete Real Sinusoid Formula I have a basic question about the formulas represented for discrete sinusoids. In some textbooks we see this form

$$x[n]= A\cos (2\pi f nT + \phi ) $$

and in others we have

$$x[n]= A\cos (2\pi k/N n + \phi ) $$

What $N$ and $T$ stands for? And what is the differences between these two formulas?

In the first formula, $x[n]$ is written as if it were a sampled version of the continuous-time function

$$x_c(t)=A\cos(2\pi ft+\phi)\tag{1}$$

If you sample $(1)$ at $t=nT$, $n\in\mathbb{Z}$, you get the first formula. $T$ is the sampling period, or the inverse of the sampling frequency $f_s=1/T$.

If you compare the two formulas in your question you see that $fT=k/N$. Assuming that $k$ and $N$ are integers, then the second formula is a special case of the first one. The discrete-time signal obtained by sampling $(1)$ is not necessarily periodic. However if $fT$ is rational, as in your second formula, then $x[n]$ is periodic with period $N$.


how do I calculate SNR of noisy signal? I am having problems in understanding how to do it practically

I have a wav file that contains pure speech and another ave file that just contains the background noise (can be various things, such as white noise, crowd noise, a recording of blowing wind, etc). These are either purely speech or purely noise. So I suppose I can get an SNR value from them by simple division of corresponding samples (or the average of a frame of samples) in both files. I then combine them in audacity to get a noisy speech file. I am guessing this file will still have the same SNR.

Now I pass this file through my noise reduction program and get another file as a result. How do I calculate the SNR of this "noise reduced" file?

--- EDIT ---

follow-up question posted HERE

The common definition of SNR is the power of the wanted signal divided by the noise power. Suppose you have obtained the wanted and the noise signal as arrays, calculation of the SNR in Matlab before noise reduction can be done like this:

snr_before = mean( signal .^ 2 ) / mean( noise .^ 2 );
snr_before_db = 10 * log10( snr_before ) % in dB

After noise reduction, the residual noise can be calculated as the difference of the wanted signal and the actual signal. Calculation of SNR is then straightforward:

residual_noise = signal - noise_reduced_signal;
snr_after = mean( signal .^ 2 ) / mean( residual_noise .^ 2 );
snr_after_db = 10 * log10( snr_after )


DAC to Analog Filter to ADC As the order of filter needed for digital information gets higher, it takes longer to process and requires more memory. Maybe analog filtering could be of use here. I know that the words analog and high precision dont often appear in the same room together, but lets assume they do for the sake of this post.

Ok, so I suggest that when a higher order filter is needed the following linkage of devices could be implemented:

DAC > Analog Filter > ADC

My question regarding this are:

Has this linkage of devices been used before?
If so, where can I find literature on it?
If not, why do you think it hasn't been used?

(Assume highly precise analog filters are being used)

You wouldn't normally do that because it isn't worth the cost (hardware.)

FIR filters get computationally very expensive if you need sharp cutoffs.

Digital IIR filters that provide the same cutoff as a long FIR are computationally very cheap.

An analog filter will be equivalent to a digital IIR.

If an analog filter can provide the performance you need, then you are better off implementing a digital IIR filter.  It saves you the DAC/ADC parts and it saves you the mess of an additional analog section with its inherent instabilities and potential to pickup noise.

If the digital IIR that can provide the performance you need is unstable, then it is very likely that the equivalent analog filter will also be very difficult to make work correctly.


FIR analog low-pass filter? Is there a way to construct FIR analog low-pass filter? (By this, I believe I am asking for no-pole analog low-pass filter, in terms of Laplace transform.)

You have to clearly define what you mean by an "analog FIR filter". "No poles" is not correct because (discrete-time) FIR filters do have poles; they are just all at the origin of the $z$-plane (for causal FIR filters). Note that filters without poles do not exist. Take as an example the discrete-time transfer function

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

with a zero at $z=a$. Obviously, $H(z)$ has a pole at $z=0$. Similarly for

$$H(z)=z-a\tag{2}$$

which also has a zero at $z=a$, but a pole at $z=\infty$.

If you mean transversal filters implemented as a tapped delay line, then you can find analog versions of it, see e.g. this paper.

If you literally mean a filter with a finite impulse response (FIR) then, technically speaking, such a thing does not exist in the continuous-time world. But of course, any stable system has an impulse response that decays sufficiently fast such that for practical purposes it may be considered as zero (or unmeasurable) after a certain time. The discussion if such a system should be considered as FIR or not is moot.


FIR filter design for complex signal I have a complex signal generated by an impedance analyzer.

What is the best approach for designing a low pass FIR filter for this?
Is a real filter applied separately to the real and imaginary streams optimal for this or do I need a specialized algorithm for complex filter design?

In fact you have two signals, and it depends on what you want to achieve, but usually you would just filter both signals (the real and the imaginary part) with the same (real-valued) low pass filter. So you either need two (identical) low pass filters, or you filter both signals sequentially with the same filter.

In the general case, filtering a complex signal $x(t)=x_R(t)+jx_I(t)$ with a complex impulse response $h(t)=h_R(t)+jh_I(t)$ requires four real-valued filtering operations (just like complex multiplication requires four real-valued multiplications):

$$\begin{align}y(t)&=(x\star h)(t)\\&=(x_R\star h_R)(t)-(x_I\star h_I)(t)+j\left\{(x_R\star h_I)(t)+(x_I\star h_R)(t)\right\}\end{align}\tag{1}\\$$

In your case $h_I(t)$ is zero, so you're left with just two filtering operations.

One final remark on complex convolution: actually, one can get away with only 3 real-valued convolutions (just like complex multiplication really needs only 3 multiplications if you're smart):

$$\begin{align}y_1(t)&=(x_R\star h_R)(t)\\
y_2(t)&=(x_I\star h_I)(t)\\
y_3(t)&=((x_R+x_I)\star (h_R+h_I))(t)\\
y(t)&=y_1(t)-y_2(t)+j\left\{y_3(t)-y_1(t)-y_2(t)\right\}\end{align}\tag{2}\\$$


Significance of energy and power signals in real world In signals and systems subject we study about calculating energy and power of a signal, but what is their practical significance.

Also what is the significance of a signal which is neither energy nor power signal ?

I mostly agree with @PeterK.'s answer; in practical terms, all signals are energy signals. However, a signal's energy does have important practical significance. Producing a signal of a given energy consumes at least that same amount of energy, for which you have to pay. In a battery-operated device, this becomes very important, because batteries store energy.

If you're using a signal to communicate, then the more energy it has, the farther it will reach (I'm oversimplifying a bit, but in general this is true). The more energy a signal has when it is received, the easier it is for the receiver to recover the information in the signal without (or with few) errors.

Many devices are also limited in the amount of energy they can produce per unit of time, which limits the kinds of signals they can produce.


Reduce the "glare" effect of a receiver within a system I'm working this summer on a very interesting, albeit difficult to understand at 100 %, project and it turns around reducing the glare effect than could be applied on a receiver when he gets signals coming from two emitters sending $X(t)$ (more energy for this signal compared to the other) and $s(t)$ (much weaker).

Let me put it in other words.

The receiver gets a signal mixture $Y(t)=H(t)*X(t)+g(t)*s(t)$ but we don't know $H$ nor $X$ nor $g$ nor $s$ nor the frequencies of the emitted signals. The only thing we know is $Y(t)$ and the fact that their spectra are dissociated in the best cases.

A picture speaks a thousand words, so...

Basically, the signal $X(t)$ overpowers $s(t)$. So I was thinking of a way to get the $s(t)$ and the $X(t)$ and with the few info I've gotten, I couldn't implement the generic solutions. However, I still have one card which is the FASTICA but that would suppose having many samples of "mixed " signal to extract the two important signals. What do you think about it?  I'm open to any software and/or mechanic suggestion.

I've read an interesting paper on this but they were talking about the near-far effect for GPS signals.
Thanks for the help!!

If you want to separate two sources you need some kind of diversity. Your know almost nothing about these signals so the only diversity you can exploit is frequency one. So first both spectrum should not be overlapped. And the second: as you know not all transmitted energy of the strong source (that you called $X(t)$) is allocated in the band you draw. Some out of band emission exists and its power depends on shape filter and power amplifier. So if even $X(t)$ and $s(t)$ are not overlaped at all, $s(t)$ may be totally shaded by out of band emission of $X(t)$. And you can't filter it out. You may refer to the problem of multiple access interference (MAI) problem of multi-user systems for details.

Now suppose frequency diversity is big enough and no out of band emission of $X(t)$ interferes with $s(t)$. As you were told FFT can be used to allocate both of the signals and further we can do bandpass filtering to completely distinguish them. But you probably need a large amount of data to achieve high resolution with FFT approach. Though you can improve FFT performance by windowning data set. The main advantage of FFT in this case is its ability to allocate signal even at very bad SNR (0 or less). You only need to increase its size.

The alternative is using some correlation properties and second order statistics. As I understand Blind source separation rely on it (PCA, ICA, but I can mistake). You may estimate input signal autocorrelation matrix as $R = E[xx^H]$. Further if you make $svd(R)$ you'll get singular eigen values of the signal mixture. Using these methods you can separate sources only if their eigenvalues are distinguishable from each other and noise ones (noise subspace). If the diversity is low or signal is too weak its eigen value decreases and it becomes "invisible" to the algorithm. The main advantage is much less amount of data is needed to allocate sources in comparison with FFT approach.

It is also important how large band of search is. It restricts complexity of the computations required. Furthermore $svd$ based methods (e.g. spectrum estimation) are good for narrowband sources (as I've seen, can mistake) while FFT is well enough for both narrow- and broadband cases.


What is the difference between Constant-Q Transform and Wavelet Transform and which is better i have been seeking for a better Transform than STFT with high overleaping. The Transform should more suitble for the human auditory system.

And i learned that there are at least 3 Methods i may use like MFCC, CQT and CWT.
Just ignore MFCC and considering of CQT and CWT, i have red some paper and thesis, but can't really tell which is better. Some of the papers i red says that CWT is only special case of CQT, and CQT is better. Other says CQT is the special case of CWT. (like this question here :Difference between CQT and WT) Some Paper said that the CQT is not invertible because some samples never get analysed

So i'm now confused, what is the difference between CQT and CWT? Which is better if i need both time and spectral analyse? And if im seeking for a good transform for both speech and music analyse, which transform would be best for me?

And a stupid question: do these Transforms need overleap in order to get better result?

A Constant Q transform is a variation on the DFT. In other words, it is a type of wavelet transform.
I only have a casual understanding of both types of transforms myself, so take what I'm saying with a grain of salt.

A standard DFT uses a constant buffer size throughout all frequencies. This typically leads to a pretty consistent, fully continuous transform. However, the constant bin size for all frequencies leads to some problems when you map frequency on a logarithmic scale. Specifically, peaks on the lower end are incredibly wide (sometimes up to half an octave), lacking any sort of detail.

This is an issue for emulating human perception because humans perceive frequency on a logarithmic scale.

A Constant Q transform seeks to solve this problem by increasing the buffer size for lower frequencies, and alleviate some of the computational strain caused by this by reducing the buffer size used for high frequencies. It's pretty effective at this, but has a few drawbacks.

The O of a Constant Q transform only slightly larger than that of a standard DFT, but because the buffer size changes per frequency, it is impossible to introduce the optimizations an FFT has to a Constant Q transform.

In other words, a Constant Q transform will yield better results where low frequencies and logarithmic frequency mapping are concerned, but it's extremely difficult to make it run realtime and has slightly less detail in the far upper frequencies.

P.S. If I am wrong about any of this, please correct me.


Measure the amount of drift in video and audio frames I setup a stereo system using 2 GoPro Hero Session Cameras, since I need a larger baseline than the one of the Dual Hero System.

This unfortunately means the two cameras are not genlocked and as a result I expect there to be some drifting between the frames of the two video/audio signals. I am using the GoPro remote to start and stop the videos at the same time, but this is obviously not going to solve the problem.

The nominal camera frame rate and audio sampling rate are 60Hz and 32KHz respectively.

I extracted the two .mp4 and .wav files using the python libraries imageio and scipy, and as expected the number of frames don't match.

import imageio
import scipy.io.wavfile as wavereader
video_frames_l = imageio.get_reader(video_l_path,  'ffmpeg')
video_frames_r = imageio.get_reader(video_r_path,  'ffmpeg')
srate_l, audio_frames_l = wavereader.read(audio_l_path)
srate_r, audio_frames_r = wavereader.read(audio_r_path)

print "Left num. video frames : %d"%len(video_frames_l) #prints 20922
print "Right num. video frames: %d"%len(video_frames_r) #prints 20914
print "Left num. audio frames : %d"%len(audio_frames_l) #prints 11168768
print "Right num. audio frames: %d"%len(audio_frames_r) #prints 11164672

I'd like to know how can I quantify:

The synchronization between the two video/audio signals without any post-processing?
How much drift was there over the whole recording?

And more importantly:

How can I improve synchronization in post-processing?

During recording I went into a dark silent room and turned lights on and off and clapped my hands a few times to create some signal.

Thanks for any pointers!

the two cameras are not genlocked and as a result I expect there to be some drifting between the frames of the two video/audio signals.

Drifting would imply that there is frame rate jitter. While jitter is unavoidable, I would not expect it to be that large so as to cause sync problems. This jitter would pose a limit to the longest scene that can be shot before starting to notice the de-syncing of the two streams. Without knowing a little bit more about the characteristics of the jitter, it is difficult to get an estimate of that maximum length.

I am not sure if the "sync cable" of the dual Hero produces a common time reference for them or if it is simply a more accurate way of synchronising video start by triggering both cameras at the same time through the port connector. I suspect that it simply triggers them more accurately at the same time.

Therefore, their frame rate must be fairly accurate.

Accurate frame rate would imply a constant "time difference" between the two streams once video capture has started. Therefore, synchronising the two videos by aligning the "sound peaks" from a clicker or a clap should work adequately.

  I extracted the two .mp4 and .wav files using the python libraries imageio and scipy, and as expected the number of frames don't match.

They wouldn't because the two cameras still operate as two independent ones. Once you have aligned the two streams, you can then extract videos of some length $l$ from the common starting point and of course in that case the two streams would match in length.

  I'd like to know how can I quantify:

  The synchronization between the two video/audio signals without any post processing?
  ...

Synchronisation only makes sense when measuring from a common event.

Bring the two streams into your video editor, measure the time between the two peaks in the sound stream that are produced by the two cameras picking up the sound of the clap. That's your audio delay but only with respect to the start of the recording. Not sure how useful this might be.

If you align the streams, so that they coincide on a frame whose sound clip contains the peak of the clap and then measure their time difference, then that would be your "video" delay. It will be within $\frac{1}{FR}$ of a second, where $FR$ is your frame rate. Again, I am not sure how useful this might be if we assume that the frame rate is fairly accurate.

  ...
  How much drift was there over the whole recording?

Drift is more difficult because you need to capture enough data to reveal its structure. Drift would typically occur at times $\ll T_s$ where $T_s$ is the sampling period. The "trick" is to wait long enough until the jitter errors accumulate so much that they reveal the distortion. For example, in audio, jitter would appear to modulate the captured waveform. To observe this, you need a reference signal (i.e. a sinusoid) and some method of quantifying phase difference (i.e. a mixer or autocorrelation, or other). But, with $T_s$ in the order of tenths of kHz and with jitter expected to occur $\ll T_s$, you would have to accumulate enough errors to make up a "shift" of a few $T_s$ BEFORE you can observe a strong error signal between reference and captured. Without further assumption about jitter, this might well extend into hours of capturing.

Similar considerations apply to video. You can set up the cameras and shoot a very good time reference signal. For example, a rotating radial pattern disc on a brushless hobby RC motor with a common speed controller. In post-processing, align the videos based on their sound streams and then examine the video stream further down the line for any drifts. To do this automatically, you might want to derive a third video stream where each frame is calculated as $I[n+1]-I[n]$ where $I$ is the $n^{th}$ frame from a video stream. This "error" signal should remain constant (and ideally zero) throughout the video.

  And more importantly:

  How can I improve synchronization in post-processing?

You cannot.

However, this company seems to be doing a genlock for GoPros

Hope this helps.


Correlation as measure of similarity Normalized correlation is used to measure how similar are two signals. I can understand similarity visually but I haven't seen a mathematical definition for the similarity of 2 signals anywhere. So I have 2 big questions that are bothering me:

Is there a mathematical definition for how similar two signals are?(similar in same fashion of visually similar)
With that definition, is there a proof that shows correlation is maximum between similar signals and not any other signal?

EDIT:

Okay so I thought about 1 and I have a mathematical definition for similar signals. Basically s1 is similar to s2 if they are proportional, that is s2 = a.s1. So we can focus on question 2 now :). A proof that says out of all functions out there, the correlation of s1 and s2 is maximum only if s1=k.s2

As you correctly noted, similarity is not a rigorously defined mathematical term. However, "distance" can be defined mathematically. Quoting Wikipedia:

  In statistics and related fields, a similarity measure or similarity
  function is a real-valued function that quantifies the similarity
  between two objects. Although no single definition of a similarity
  measure exists, usually such measures are in some sense the inverse of
  distance metrics...

Here's a hand-wavy explanation, the math isn't totally precise, but hopefully it's enough to get some intuition:

Let's say we have real valued signals $x(t)$ and $y(t)$ defined on the unit interval $t \in [0,1]$. We will say that these signals are similar if the distance $||x-y||$ is small. Note that there are various choices of how to measure distance. A commonly used distance measure is the squared difference a.k.a. L2 distance defined as
$$
||x-y||^2_{L^2} \triangleq \int_0^1 (x(t)-y(t))^2 dt \ .
$$

Let's say we have three signals $x$, $y$ and $z$. Normalize all signal energies to 1 so that $\int_0^1 (x(t))^2 dt = 1,$ etc. Suppose $x$ is more similar to $y$ than $z$. This means

$$||x-y||_{L^2} < ||x-z||_{L^2}$$

or,

$$\int_0^1 (x(t)-y(t))^2 dt < \int_0^1 (x(t)-z(t))^2 dt$$

or, rearranging terms,

$$\int_0^1 x(t) y(t) dt > \int_0^1 x(t) z(t) dt \ .$$

In other words, the correlation between $x$ and $y$ is higher than that between $x$ and $z$. This shows that in order to find a "more similar" signal, maximizing the correlation is indeed the right thing to do (assuming L2 distance).


How to detach overlapping contours but not to remove small ones. Alternative to erosion-dilation cycle Background:
I'm trying to detect trees on a sequence of aerial images by calculating a stereo disparity for consequent pairs of images. This is how disparity looks like. White spots are trees:

Next, I am locating contours on the disparity map and getting the following image:

From here I would like to detect the trees.

What I'm doing right now:
First, I'm doing an erosion:

kernel = np.ones(shape=(5, 5),
                 dtype=np.uint8)

image= cv2.erode(image,
                 kernel=kernel,
                 iterations=14)

and getting some overlapping, e.g.:

In order to separate those overlapping contours, I do dilation this time:

image= cv2.dilate(image,
                  kernel=kernel,
                  iterations=12)

Ups, some contours became so small and almost disappeared.
Reversing damage:

image= cv2.erode(image,
                 kernel=kernel,
                 iterations=7)

Ok! Better! But.. all this process seems like a dirty workaround.. Moreover, this won't work for some other images so well (ignore the noise on the left, I should add cropping):

Questions:
 - How do I make sure not to overlap the contours of the trees and not to remove them while dilating/eroding?
 - Is there a better alternative to what I am doing right now?

Edit:
Probably for the full picture I should add that after obtaining a good image with separate contours, I use:

disparity, contours, hierarchy = cv2.findContours(image,
                                                  mode=cv2.RETR_LIST,
                                                  method=cv2.CHAIN_APPROX_SIMPLE)

and plot the contours on the original image:

Edit Nº2:
Adding cross-section of disparity map. Raw data available from here

Plot:

Update on results:

I tried a median filter as it was suggested by @A_A, and got a very good separation of trees from the ground:

ground = medfilt(disparity, kernel_size=125)
trees = disparity - ground
boolean_mask = (trees > 2.) & (trees < 15.)
plt.imshow(boolean_mask, cmap='gray')

The problem with this filter is the speed.
scipy.signal.medfilt took me more than 20 mins to calculate everything. But scipy has also a scipy.ndimage.filters.median_filter. I don't know what's the difference between them but results that I got from it were almost the same apart from some areas close to image borders, and it took somewhere from 2 to 5 mins.

I don't know if this is a corner case or the norm in your dataset but it is a relatively easy situation to deal with. It would be much more difficult to detect trees in an urban environment, for example, where you would probably be looking at more pattern-recognition type of solutions.

The basic idea behind this solution is surface fitting, whether to directly detect the trees or detect the "floating" level of the ground and through that information detect the trees.

The following is based on Octave and using the 1D data but generalisations to other platforms (and 2 dimensions) should be straightforward.

The first thing one can do here is kernel density estimation (KDE). In this case, you would be approximating the total available surface (both the ground and the trees) via a set of Gaussian curves. In the two dimensional case, this will give you either one Gaussian per tree or higher densities of Gaussians around the locations of trees. In one shot. No further processing required. In reality, you might find that this technique approximates your disparity image via kernel density and returns the data of a set of Gaussians that you can process further to determine where there are trees.

But, having had a quick look at the data, it contains sharp discontinuities and that might not be good for KDE convergence. While you could "blur" the disparity with a Gaussian low pass and then try to approximate it with KDE, this means that your original resolution is decreased and it might reduce your ability to discriminate between small and large trees. Anyway, this is an option too (first blur with a suitable FWHM Gaussian, then approximate with KDE).

Again, this seems a relatively easy case where you could even use simple thresholding to detect the trees. The problem there of course is that the ground height now modulates your threshold. A tree on flat ground could be anything above 25 meters but certainly not a 3 meter tall shed on top of a 22 meter hill.

So, we want to threshold but we want to take into account the ground contour. Ideally, you could fit a curve that is closer to what the ground really is (an ellipsoid) but maybe you can get away with simpler fitting schemes.

The main idea here is that we want to "undo" the modulation of the ground to then apply a simple threshold. If we flatten the ground, the 22 meter hill will disappear and then anything above 22 meters is more likely to really be a tree.

If you ignore the zeros (which I am assuming are "missing values" in your case (either high uncertainty or actually missing values)), the signal looks like this:

You can threshold this at $\approx 62$ but you can see that the second tree is lower than the first but further up the slope of what looks like a hill.

Let's flatten the hill. Flattening the hill here is basically equivalent to subtracting the slowly varying component (the hill) from the signal itself and possibly even controlling the amount of removal with a coefficient. Here, I will use simple polynomial fitting but anything that could approximate the slow varying component of the height field would do.

xrange=0:(length(b)-1);
sv = polyval(polyfit(xrange,b,2),xrange); # Fit and evaluate a second order polynomial. You might want to go higher if you have ridges or depressions but not too high because we just want to approximate the slow varying ground changes.
plot(b);hold on;plot(sv,'r');grid on; # Plot everything

This now looks like this:

Obviously, we can simply subtract the red curve from the blue and get back the "flattened" version of the signal, but as you can see, the trees have "pulled the ground up" slightly and this is why if you do this blindingly, you will get this:

In terms of comparing between the trees, we are now definitely in a better condition because we have removed the ground and we can also tell that the second tree is lower than the first, but, we now have a dip from trying to remove the hill. To adjust for situations like this, you can introduce a coefficient by which to control how much subtraction you want to perform, for example:

m=0.43;plot(m.*b-(1-m).*sv);

Here, $m \in [0 . . 1.0]$ controls how much of the real signal and how much of its approximation we are taking into account. I determined 0.43 after trial and error and this results to:

Where basically, the trees start at almost the same height.

So, now you can threshold anything above -8.2 and be more certain that you are selecting trees.

There is a way by which you can flatten the ground almost completely but it requires some knowledge about the size of the tree canopy.

This works on the same premise but approximates the ground using a different method.

Basically, you can use a median filter to "chop the peaks" that the trees are causing to the height field. Anything that remains after that is basically ground variation.
The median filter is applied in a "rolling window" way and substitutes a sample with the median value of $N$ values surrounding it.

Now, to apply this, you need to determine that $N$ which should be almost as wide as your widest tree canopy. I have determined this to be approximately 125 (pixels) in this example, so by doing:

plot(b-medfilt1(b,125)); grid on;

You get something like this:

Which is almost too good. At least the trees start at exactly the same point. Again, it is now the point of applying a simple threshold (anything above 1 (?)) to detect where the trees are.

This process will give you a binary mask which you can then shrink (or even simply cluster, with simple methods) to determine a centroid that you would then associate with the location of a tree. Obviously, processing the size of each "cluster" could probably also tell you how wide is the tree canopy.

In terms of generalising this, you can either fit a 2D polynomial or a 2D NURB or even simply apply a 2D median with medfilt.

Hope this helps.


Simulate packet loss in internet video streaming i have encoded H.264 video using ffmpeg . I want to apply packet loss model like Gilbert Eliot model to depict  video transmission over unreliable network. Like in this thread TV noise is simulated. But noise in internet is different. So is there a way to simulate packet loss/ noise for video streaming using ffmpeg or matlab?

For me your question is not clear enough. This questions come to my mind:

what exactly do you want to transmit?
Do you want to transmit raw H.264 bitstream (NAL units) or maybe MPEG2-TS, or RTP packets? Most video transmission applications would use either MPEG2-TS or RTP. I suppose that you are not interested in HTTP or other TCP based streaming solutions.
you are talking about the noise in the internet. Are you interested in simulation of some specific network topology?
what network protocols do you want to use and how does the transmission channel looks like?

You might take a look at NS-3 network simulator. With it you can simulate quite everything and it's also quite easy to use if you are familiar with C++ and network programming. But it is probably not the easiest way to solve your problem.

I'm pretty sure that you can also stream the video with FFMPEG but I don't think that FFMPEG allows you to drop packets based on some error model. But you could use NetEm to generate packet loss, delay or duplication of packets on your Linux system.


Is Z-transform of $\sin(\omega_0n)$ same as that of $\sin(\omega_0n)u[n]$ The Z-transform tables only mention the transform of $\sin(\omega_0n)u[n]$, e.g. #21 at this link: https://en.wikipedia.org/wiki/Z-transform#Table_of_common_Z-transform_pairs
But how can I find of Z-transform of a signal that is simply $\sin(\omega_0n)$? This signal is an input to a causal system. Should I assume its transform is same as $\sin(\omega_0n)u[n]$?

The two sequences $\sin(\omega_0n)u[n]$ and $\sin(\omega_0n)$ are very different, so their transforms can't be identical. The $\mathcal{Z}$-transform of $x[n]=\sin(\omega_0n)u[n]$ can be computed as

$$X(z)=\sum_{n=0}^{\infty}\sin(\omega_0n)z^{-n}\tag{1}$$

which converges for $|z|>1$ to the expression given in the table you linked to. However, the $\mathcal{Z}$-transform of $x[n]=\sin(\omega_0n)$ would be

$$X(z)\overset{?}{=}\sum_{n=-\infty}^{\infty}\sin(\omega_0n)z^{-n}\tag{2}$$

which doesn't converge anywhere, and hence does not exist.

Consequently, the $\mathcal{Z}$-transform is no suitable tool for computing the response of a discrete-time LTI system to a sinusoidal input. However, note that the discrete-time Fourier transform (DTFT) of $x[n]=\sin(\omega_0n)$ does exist. It contains Dirac delta impulses:

$$\text{DTFT}\{\sin(\omega_0n)\}=-j\pi\left[\delta(\omega-\omega_0)-\delta(\omega+\omega_0)\right]\tag{3}$$


White noise DC component I'm really new to DSP, I'm actually studying Computer Science and took an elective on DSP so my knowledge is pretty limited.

I've learned that pure white noise signal  for example has all possible frequencies, that means that it also has zero frequency so it has a non zero DC component(does that mean it?) but pure white noise by definition has a zero mean and thus has a zero DC component..
What am I missing?

Thanks!

White noise has no component at any specific frequency since it is a noise density and therefore has energy across a range of frequencies. The range of any given specific frequency, including DC is zero, and therefore the energy at DC (and any specific frequency) is indeed zero.

In practice you cannot actually observe DC or any specific frequency as a single point as that would require an infinite observation time. However, the longer you observe, the narrower will be the effective range of frequencies you observe and for this reason the mean will be non-zero but approach zero as your observation time approaches infinity.


Developing DSP software, am I 10 years too late? I'm interested in developing software related to digital signal processing. The most obvious ideas are audio and video editing. There's an overwhelming amount of audio and video editing software out there. I feel like developing such software would be in vain.

What are the most lucrative frontiers in DSP software today?
Do you think audio/video software still have room for improvement?

DSP is a wide field and goes beyond audio/video processing. And, since new digital devices are developed in any field, DSP algorithm implementation, simulation and design is still a very important area.

Some examples that come to my mind beyond audio/video:

wireless communications: for example there is plenty of research going on in 5G cellular (i.e. mobile phones) and 802.11 (i.e. WiFi) standards.
medical devices: more and more devices operate (semi-)automatically, requiring sophisticated signal processing to interpret their measurements.
DSL: as users search for faster internet, the research on DSL improvement is very important.
Big Data and deep learning: On the one hand, it's about designing low-complexity algorithms to solve the problems of huge sizes. On the other hand, it's about providnig hardware that can do these operations quickly.

This list is not at all exhaustive, and there are many more fields, where DSP is absolutely necessary and new techniques are developed on a daily basis.


Why do I have a high peak at the beginning of my FFT (not DC)? I am analyzing a recording from a microelectrode array on a neuronal cell culture and I found something that is disturbing me a lot. Through python, I have taken different chunks of the same signal profile from the original recording, calculated the FFT of each of them and then, the mean among them.
In theory, this makes sense since I am just amplifying common patterns on the wave, however, I am getting a huge peak at the 0 component but if you calculate the mean of the signal, they have zero DC component. I don't know if this issue might be a consequence of a bad use of the FFT on python or maybe I am missing something. Do any of you have any idea of what is going on?

The code is very long but the specific section of the chunks and the fft is here:
chunk = np.zeros((16, 300000))
chunk[0,:] = signal1[320000:620000]
chunk[1,:] = signal1[950000:1250000]
chunk[2,:] = signal1[1640000:1940000]
chunk[3,:] = signal1[2311000:2611000]
chunk[4,:] = signal2[320000:620000]
chunk[5,:] = signal2[950000:1250000]
chunk[6,:] = signal2[1640000:1940000]
chunk[7,:] = signal2[2311000:2611000]
chunk[8,:] = signal3[320000:620000]
chunk[9,:] = signal3[950000:1250000]
chunk[10,:] = signal3[1640000:1940000]
chunk[11,:] = signal3[2311000:2611000]
chunk[12,:] = signal4[320000:620000]
chunk[13,:] = signal4[950000:1250000]
chunk[14,:] = signal4[1640000:1940000]
chunk[15,:] = signal4[2340000:2640000]

myfftchunk0 = np.fft.fft(chunk[0,:])
myfftchunk1 = np.fft.fft(chunk[1,:])
myfftchunk2 = np.fft.fft(chunk[2,:])
myfftchunk3 = np.fft.fft(chunk[3,:])
myfftchunk4 = np.fft.fft(chunk[4,:])
myfftchunk5 = np.fft.fft(chunk[5,:])
myfftchunk6 = np.fft.fft(chunk[6,:])
myfftchunk7 = np.fft.fft(chunk[7,:])
myfftchunk8 = np.fft.fft(chunk[8,:])
myfftchunk9 = np.fft.fft(chunk[9,:])
myfftchunk10 = np.fft.fft(chunk[10,:])
myfftchunk11 = np.fft.fft(chunk[11,:])
myfftchunk12 = np.fft.fft(chunk[12,:])
myfftchunk13 = np.fft.fft(chunk[13,:])
myfftchunk14 = np.fft.fft(chunk[14,:])
myfftchunk15 = np.fft.fft(chunk[15,:])

Meansignal = np.zeros(len(chunk[1, :]))

for i in range(len(chunk[1, :])):
    Meansignal[i] = myfftchunk0[i] + myfftchunk1[i] + myfftchunk2[i] + myfftchunk3[i] + myfftchunk4[i] + myfftchunk5[i] + myfftchunk6[i] + myfftchunk7[i] + myfftchunk8[i] + myfftchunk9[i] + myfftchunk10[i] + myfftchunk11[i] + myfftchunk12[i] + myfftchunk13[i] + myfftchunk14[i] + myfftchunk15[i]
    Meansignal[i] = Meansignal[i] / 16

mynewffthalfMean = Meansignal[:len(Meansignal) // 2]
fnewHalfMean = np.linspace(0, fs / 2, len(mynewffthalfMean))

# Plot figure 2 - FFT
plt.figure(7)
plt.plot(fnewHalfMean, np.abs(abs(mynewffthalfMean))**2)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Magnitude (dB)')

You cannot simultaneously have a non-zero value at the first bin of the DFT (the 0 component) and a zero-mean since the first component is by definition of the DFT the (scaled) mean of the signal:
from the DFT:
$$X[k] = \sum_{n=0}^{N-1}x[n]e^{-j2\pi nk/N}$$
$$X[0] = \sum_{n=0}^{N-1}x[n]e^{j0}= \sum_{n=0}^{N-1}x[n] = N\mu_x$$
So there is a mean in the digital data, and perhaps the OP's confusion is there is not a mean of any significance in the analog waveform being sampled. There are several opportunities to introduce offsets notably in the ADC converter itself, or it may be as simple as the mapped amplitude range of the digital waveform includes an offset. But taking the mean of the data that is being FFT'd should be consistent with what is above, with further scaling by any windowing that may be done according to $\sum_{n=0}^{N-1}w_n$.


Positive and negative frequencies in DFT due to frequency folding, or due to negatively indexed frequencies? When I look for the cause of the mirroring of frequencies in DFT output, I get two types of explanations:

The first one which says the frequencies are mirrored because of the complex exponential which has a positive and negative counterpart. The poritive(real) parts add up and the negative (imaginary) parts cancel out, and that the amplitudes are spilt equally between the positive and negative frequencies.
The second one which says this happens as a consequence of sampling above the Nyquist frequency, i.e., above $\frac{N}2$. But the mirroring would still be present if we take the range $-\frac{N}2$ to $\frac{N}2$.

My question is what is the actual cause for this mirroring, and what actually changes with taking the range as $-\frac{N}2$ to $\frac{N}2$ or $0$ to $N-1$?

I think that considering the DFT from a linear algebraic point of view has some value, so I will attempt to introduce the foundations.

We will assume that our signal is a vector of $N$ complex entries.

$\mathbb{C}^N$ is the vector space of vectors with $N$ complex entries. Let $\mathbf{u}_0,\mathbf{u}_1,\ldots,\mathbf{u}_{N-1}$ be vectors in $\mathbb{C}^{N}$ defined by
\begin{equation}
\begin{split}
\mathbf{u}_k
&=~
\frac{1}{\sqrt{N}}\left(\begin{array}{c}
\exp(2\pi \mathsf{j}\times 0\times(k/N))\\
\exp(2\pi \mathsf{j}\times 1\times(k/N))\\
\exp(2\pi \mathsf{j}\times 2\times(k/N))\\
\exp(2\pi \mathsf{j}\times 3\times(k/N))\\
\vdots\\
\exp(2\pi \mathsf{j}\times (N-1)\times(k/N))
\end{array}\right)\\
&=~
\frac{1}{\sqrt{N}}\left(\begin{array}{c}
1\\
e^{2\pi\mathsf{j}k/N}\\
(e^{2\pi\mathsf{j}k/N})^2\\
\vdots\\
(e^{2\pi\mathsf{j}k/N})^{N-1}
\end{array}\right),
\end{split}
\end{equation}
for $k = 0,1,2,\ldots,N-1$, where $\mathsf{j} = \sqrt{-1}$.

Every entry of $\mathbf{u}_0$ is $1/\sqrt{N}$, so $\mathbf{u}_0$ might be considered as a sampled DC signal.
The entries of $\mathbf{u}_1$ are samples of a complex exponential with fequency $\frac{1}{N}$,
The entries of $\mathbf{u}_2$ are samples of a complex exponential with fequency $\frac{2}{N}$,
and so on, up through frequency $\frac{N-1}{N}$.

$\mathbf{u}_0,\mathbf{u}_1,\ldots,\mathbf{u}_{N-1}$ form an orthonormal basis for $\mathbb{C}^{N}$, which means that each $\mathbf{u}_k$ has norm 1, they are all orthogonal to one another, and each vector in $\mathbb{C}^{N}$ can be represented unambiguously as a linear combination of them. An important result of this is that, if $\mathbf{x}\in\mathbb{C}^{N}$, then there is exactly one list of complex numbers $c_0,c_1,\ldots,c_N$ such that
\begin{equation}
\mathbb{x} = c_0\mathbf{u}_0 + c_1\mathbf{u}_1 + \cdots + c_{N-1}\mathbf{u}_{N-1}.
\end{equation}

The coefficients mentioned above are the entries of the DFT of $\mathbf{x}$:
\begin{equation}
\mathbf{x} = X[0]\mathbf{u}_0 + X[1]\mathbf{u}_1 + \cdots + X[N-1]\mathbf{u}_{N-1}.
\end{equation}
We might interpret $X[0]$ as the strength of the DC component of $\mathbf{x}$, $X[1]$ as the strength of the component of $\mathbf{x}$ with frequency $\frac{1}{N}$, and so on. Since $\mathbf{X} = \mathsf{DFT}\mathbf{x}$ has complex entries, there is some phase information attached to each "strength".

So far, we have considered only components of non-negative frequencies. What if we would rather view $\mathbf{x}$ as a combination of negative and positive frequency components? Consider a component of frequency $-\frac{k}{N}$ for $0< k \leq \frac{N}{2}$:
\begin{equation}
\mathbf{u}_{-k}
~=~
\frac{1}{\sqrt{N}}\left(\begin{array}{c}
\exp(2\pi \mathsf{j}\times 0\times(-k/N))\\
\exp(2\pi \mathsf{j}\times 1\times(-k/N))\\
\exp(2\pi \mathsf{j}\times 2\times(-k/N))\\
\exp(2\pi \mathsf{j}\times 3\times(-k/N))\\
\vdots\\
\exp(2\pi \mathsf{j}\times (N-1)\times(-k/N))
\end{array}\right).
\end{equation}
The $\ell^{\textrm{th}}$ entry of this vector is
\begin{equation}
\begin{split}
u_{-k}[\ell] &=~ \frac{1}{\sqrt{N}}\exp\left(2\pi\mathsf{j}\times\ell\times\frac{-k}{N}\right)\\
&=~ \frac{1}{\sqrt{N}}\exp\left(2\pi\mathsf{j}\times\ell\times\frac{-k}{N}\right)\times\underbrace{\exp\left(2\pi\mathsf{j}\times\ell\times\frac{N}{N}\right)}_{1}\\
&=~ \frac{1}{\sqrt{N}}\exp\left(2\pi\mathsf{j}\times\ell\times\frac{N-k}{N}\right)\\
&=~ u_{N-k}[\ell].
\end{split}
\end{equation}
In other words, the negative-frequency component $\mathbf{u}_{-k}$ is exactly the same as the positive-frequency component $\mathbf{u}_{N-k}$.

Suppose that $N = 2M$ for some positive integer $M$. Then
\begin{equation}
\begin{split}
\mathbf{x} &=~ X[0]\mathbf{u}_0 + \cdots + X[N/2-1]\mathbf{u}_{N/2-1} + X[N/2]\mathbf{u}_{N/2} + \cdots + X[N-1]\mathbf{u}_{N-1}\\
&=~ X[0]\mathbf{u}_0 + \cdots + X[M-1]\mathbf{u}_{M-1} + X[M]\mathbf{u}_{M} + \cdots + X[N-1]\mathbf{u}_{N-1}\\
&=~ \underbrace{X[0]\mathbf{u}_0 + \cdots + X[M-1]\mathbf{u}_{M-1}}_{\textrm{non-negative-frequency components}}
+ \underbrace{X[M]\mathbf{u}_{N-M} + \cdots + X[N-1]\mathbf{u}_{N-1}}_{\textrm{negative-frequency components}}
\end{split}
\end{equation}
For a full decomposition, one can choose the frequency-sets
\begin{equation}
-\frac{N/2}{N},-\frac{N/2-1}{N},\ldots,-\frac{1}{N},0,\frac{1}{N},\ldots,\frac{N/2-1}{N}
\end{equation}
or
\begin{equation}
0,\frac{1}{N},\ldots,\frac{N-1}{N},
\end{equation}
each of which consists of $N$ distinct frequencies. In truth, one can choose other frequency-sets of $N$ frequencies, too, but these are the ones to which we have attached some intuition over the decades.

MATLAB's fft gives the DFT with all non-negative frequencies. To convert the output of fft to the vector of coefficients for negative, zero, and positive frequencies, one applies the fftshift function.

All of this and much more is explained from a linear algebraic point of view in Linear Algebra, Signal Processing, and Wavelets - A Unified Approach (PDF) by Øyvind Ryan of the University of Oslo.


Estimating variance in arbitrary, periodic signal I have a periodic signal $x[m], m \in [0;M-N+1]$ made of modulated templates $s[n],~ n \in [0;N-1],~ N \ll M = NK$ of finite energy and support (i.e. zero outside of its defined interval, which does not extend to infinity). For simplicity, we write thus

$$
x[m] = \sum_{k=1}^{K}A_{k}s[m-n_{k}] + w[m]
$$

with $n_{k}$ the $k$-th repetition of the template modulated by $A_{k}$. The noise term $w[m] \sim \mathcal{N}(\mathbb{0},\sigma^2_{w})$ sums over the whole signal (not just one template) and can be considered i.i.d for each realization.

Is it possible to estimate the variance $\sigma^2_{w}$ of the noise, given this setup? For one, I have seen methods that try to obtain some portion of the signal $x$ where the template $s$ is absent, but in my practical case, that is impossible to obtain. What DSP tricks can be applied in order to obtain this information?

Do a cross correlation with your known $s$; you will get peaks of height $\|s\|A_k$ at the positions of your shifts.

You just sum up these peaks, call that sum $P$, calculate the overall energy in your signal (i.e. $\|x\|$), and calculate $\text{SNR}=\frac P{\|x\|-P}$ if you need the SNR. If you just need the noise: $\|x\|-P$ (quite directly). Done!


Minimum length of windows and separable signals I created a signal with lets say $\omega_1$ that is fixed at 1 bin in the FFT: $$A_1\cos(\omega_1t)$$

I was asked to find the closest $\omega_2$
to $\omega_1$ so that the signals $A_1\cos(\omega_1t)$ and $A_2\cos(\omega_2t)$ are separable.

I only know the ratio $A_1/A_2$ in dB and the $\omega_1$ that I chose.

$$s(t)=A_1\cos(\omega_1t)+A_2\cos(\omega_2t)$$

Now my questions are:

What is the minimum window length?
Number of points for the DFT?
What is the relationship (if any exists) between the ratio $\frac{A_1}{A_2}$ and the DFT points?

I don't know if you can have a formula in terms of window length, number of DFT points and difference between two angular frequencies you are trying to resolve. But you can do some experiments with different types of the windows and determine this in MATLAB. Note you need to consider the following.

For a DFT of size N you get two frequency bins separated by $\frac{2*\pi}{N}$. So you should have N large enough resolve $\omega_2 - \omega_1$
Note that before you take DFT, when you window it the frequency response will be convolution of two frequency response - frequency response of the window and the frequency response of the signal which in your case pretty much a pair of delta function.

Now Hanning window will have a wide main lobe and suppressed side lobes. So if the second signal's frequency falls within the main lobe then you won't be able to detect it unless amplitude of the second signal is large enough compared to the smeared component of the first signal at the location of the second signal. Or if the second signal falls in some of the side lobe location and second signal's amplitude is large enough compared to the first signal's component in side lobe region you should be able to detect it.

If you chose a rectangular window the main lobe will be substantially shorter compared to rectangular window. But the side lobes will be substantially larger compared to the Hanning window case.

There exist the trade-off between using different window sizes.


Recovering a Differential Equation From the Transfer Function of a Cascaded System With respect to the below discussion, consider that we are talking about LTIC systems characterized by constant coefficient ODEs.

Consider a cascaded system whose transfer function H(s) is given by

$$H(s) = \frac {s+2}{(s+2)(s+3)(s+1)} =\frac 1 {(s+1)(s+3)}$$

"Recovering" the defining differential equation from this transfer function yields

$$(D^2 + 4D + 3)y = P(D)x$$

where $P(D)$ is unknown but irrelevant for the purposes of this discussion.
What are the characteristic modes of the system?

Is it acceptable to simply say that the natural response $y_n$ is given as below?

$$y_n = Ae^{-t} + Be^{-3t}$$

It would seem to me that the answer is no. Why can we recover an ODE from a transfer function in general. For instance, what if the cascaded system had've had $(s-2)$ factors rather than $(s+2)$ factors that cancelled (ie. one of the subsystems had a pole at $s=2$). In such a case, would it not be wildly incorrect to recover the ODE above since then we are ignoring an internal characteristic mode that blows up $(e^{2t})$?

You have discovered pole-zero cancellation, which is one of the reasons that state-space systems description was invented.  Another reason is because when you go from a differential equation to a transfer function, you lose information about the structure of a system that you can choose to preserve in a state-space description.

In general, you should be very careful when you see a pole-zero cancellation, because even if the pole is stable, it can still affect the behavior of the system on startup, if there is a transient that pushes the system into nonlinear behavior and excites the extra pole, or if the pole and the zero don't actually match up perfectly in the real world.

And yes, if you have an unstable pole that's cancelled by a zero then Bad Things will happen to the system that simply aren't described by the transfer function.


Mapping the Value of a Sample in a 2D DFT to Cycles/Pixel If I have an image and its 2-D DFT of that image, what is the mapping between the value of the DFT at (u,v), and the frequency in the spatial domain in the x and y components, in cycles/pixel? I want to be able to say "There is a peak in the DFT at (u,v), which corresponds to this frequency". I don't think the answer is u cycles/pixel, since the values shift with different image sizes, even if the image content remains the same.

Well the Nyquist frequency is 0.5 cycles / pixel and the range of frequencies can be considered to range from -0.5 to +0.5 cycles / pixel. So if the dimension in a given axis is N pixels then each pixel is 1 / N cycles per pixel.

Depending on your particular DFT/FFT implementation you may have zero frequency in the middle of the image (M/2, N/2) or in the corner (0, 0). For the former case the spatial frequency at u, v will be:

(M/2 - u) / M, (N/2 - v) / N


Frequency prewarping of a bilinear transform (Tustin transform) sadly I'm not too well aquinted with discretazation methods. At the moment I struggle to reproduce a bilinear transform with frequency prewarping.

I have the following transfer function in the s-domain:

$$\begin{align}
\\
G(s)=\frac{-ms^2}{s^2+2abs+b^2}
\\
\end{align}$$

There m,a,b are just some constants. Now I would like to transform this transfer function into the z-domain. In Matlab I am using the c2d-function. I compared all the available methods with a sampling frequency of 512Hz.

One can clearly see, that with a pole being around 100Hz that is already too close to the Nyquist frequency, so that the Tustin method has to be wrapped in order to give a more accurate result.

Now I tried to reproduce the results from the matlab c2d-function by doing the bilinear transform myself.

I first tried it without prewarping using:

$$\begin{align}
\\
s=A\cdot\frac{z-1}{z+1},  A=\frac{2}{T}
\\
\end{align}$$

Comparing my result with Matlab does look really nice:

Then I tried it with prewarping, using:

$$\begin{align}
\\
s=A\cdot\frac{z-1}{z+1}, A=\frac{b}{\tan(b\cdot\frac{T}{2})},  b=2\pi\cdot100
\\
\end{align}$$

This results in:

With the warping the result doesn't look the same as with the matlab c2d-function. So where is the difference??

In matlab I use:

G(z)=c2d(S,1/512,['Method','tustin', 'PrewarpFrequency', b]);

And I have derived:

$$\begin{align}
\\
G(z)=\frac{-mA^2\cdot z^2+2mA^2\cdot z-mA^2}{(A^2+2abA+b^2)\cdot z^2+(2b^2-2A^2)\cdot z+(A^2-2abA+b^2)}
\\
\end{align}$$

If I put in $A=2/T$ I achieve to get the same result as the Matlab c2d-function. But with $A=b/\tan(\frac{b\cdot T}{2})$ I don't get the same result as the c2d-function in Matlab.

Putting in some values for the constants I get he following:

Tustin wo prewarp using c2d yields: $\frac{-6.781 z^2+13.56 z - 6.781}{z^2-0.8456 z + 0.8669}$

Tustin with prewarp using c2d yields: $\frac{-7.997 z^2+15.99 z-7.997}{z^2-0.615 z+0.8217}$

My self derived Trustin-transform with prewarp yields: $\frac{-6.216 z^2+12.43 z-6.216}{z^2-0.6266 z+0.8599}$

So apparently I am doing something different than the matlab function. Maybe someone can enlighten me by pointing out where I am at odds. I also failed to recreate the discretized transfer function using the matched-method from c2d. Basically my problem is that I need to build nice discrete transfer functions while being unable to use matlab. So I would like to understand how the algorithm works. Thank you in advance for any help.

There might be some undocumented features for the c2d function. Namely, if I follow the documented way of specifying pre-warping (using c2dOptions) one would get:

  opts = c2dOptions('Method', 'tustin', 'PrewarpFrequency', b);
  G = c2d(S, T, opts);

When I used this I do get the same result as when I would manually do pre-warping using:

  z = tf('z', h);
  sp = b / tan(b * T / 2) * (z - 1) / (z + 1);
  G = -m * sp^2 / (sp^2 + 2 * a * b * sp + b^2);

However, the notation that you use is not specified in the documentation for as far I could find. It does seem to give a better approximation of the continuous transfer function in this case. However, if I use the following transfer function

$$
P(s) = \frac{s}{s^2 + 0.1\,s + 1},
$$

with a sample time of $T = 1$ your method behaves in an unexpected way as shown in the figure below. In that figure D1 uses pre-warping using c2dOptions, D2 uses your method of pre-warping and D3 uses manual pre-warping (this completely overlaps with D1 in the figure).

It can be noted that D2 does follow the overal shape of $P(s)$ better, but it does not match it at the specified pre-warp frequency and the phase is much lower (when zoomed out to a bigger frequency spectrum it looks like a delay). On the other hand D1/D3 exactly matches $P(s)$ at the pre-warp frequency as expected, but the overall shape does deviate more the further away you get from the pre-warp frequency.


Tone generation in Matlab (scientific paper) I'm trying to generate sounds as described in a scientific paper that I found on the Internet.
What I need to do is to generate simple tones in a certain pattern, however I have some issues understanding how to do this correctly.
The attached image contains the full description of the method. I'm stuck at  "Thirteen tones of this sound .." .I'm confused about how this thirteen tones are selected and how to arrange them correctly.
If anyone could provide some help, it would be much appreciated.

My code so far:

function testaudio()

sampleRate = 44100;
duration  = 0.5;
freqs = [100 200 300 400 500 600 700 800 900 1000 1100 1200];

[signal t] = generateTone(freqs, sampleRate, duration);
sound(signal, sampleRate);

%generates a tone (pure or complex)
function [y t] = generateTone(freqs, sampleRate, duration)
    timeSamples = duration * sampleRate;

    t = linspace(0, duration, timeSamples);
    y = zeros(size(t));
    env = generateEnvelope(t);

    for i=1:length(freqs)
        temp = env.*sin(2*pi*freqs(i)*t);
        y = y + temp;
    end

    % amplitude normalisation
    y = y' / max(abs(y));

end

function y = generateEnvelope(t, alfa)
    alfa = 1;
    y = (-t.^alfa) .* log(t./max(t));
    y(1) = 0;
end

end

The section on 12 added sidetones is referring to something different to the section on thirteen notes.

In regards to 13 notes: By diatonic I think the mean a diatonic musical scale, with the lowest note at 148.83Hz.

From this site:  http://www.phy.mtu.edu/~suits/notefreq446.html   148.83Hz is note $D_3$

So I would choose 13 consecutive notes in the D-scale starting with $D_3$. The notes are played according to Figure 1. Each note lasts 8.40s/72 = 7/60 s.

However the notes look equally spaced in frequency in Figure 1. So I think they might mean a chromatic scale with equal temperament. That would make more sense as the 13th note would be one octave above the first note, and each note would be $^{12}\sqrt{2}$ apart. Code then:

low_note = 148.83;
notes = low_note .* nthroot(2,12).^(0:1:12);

signal = [];

% First section
for k = 1:12
    signal = [signal, generateTone(notes(13-k+1), 44100, 7/60); % High note
    signal = [signal, generateTone(notes(k), 44100, 7/60); % Low note
end
% Second section
for k = 13:-1:2    % Reverse
    signal = [signal, generateTone(notes(13-k+1), 44100, 7/60); % High note
    signal = [signal, generateTone(notes(k), 44100, 7/60); % Low note
end
% Third section
for k = 1:12
    signal = [signal, generateTone(notes(13-k+1), 44100, 7/60); % High note
    signal = [signal, generateTone(notes(k), 44100, 7/60); % Low note
end


Convolution for non-signal-processing background I am a civil engineer and am analyzing traffic data recorded by capturing vehicle movements over a highway for a specified time period. The database I am dealing with contains observations at every 0.1 seconds (total about 4 million rows) and has some 'noise' in position (lateral and longitudinal) columns. A research paper suggests to 'smooth' the velocity and acceleration data (derived from position) by using following formula:

Someone told me it is a 'convolution'. Since I have no background in electronics/ DSP, I had no idea what convolution is. After studying some online resources I can understand a little about the concept but I need a thorough understanding of this topic from a resource which is targeted towards people who have no background in DSP. Could anyone please guide me to the required resources? I am using R for analysis.

I disagree with the assertion in the answer by user7090 that the calculation in question is not a convolution; it is indeed a convolution.  Let us consider the case $D=2$. Ignoring $Z$ for now,
the five terms in that sum
$\displaystyle \sum_{k = i-2}^{i+2} x(t_k)e^{-|i-k|/\Delta}$ on the left in the OP's question are
$$\begin{align}
x(t_{i-2})e^{-|i - (i-2))|/\Delta} &= x(t_{i-2})e^{-2/\Delta}\\
x(t_{i-1})e^{-|i - (i-1))|/\Delta} &= x(t_{i-1})e^{-1/\Delta}\\
x(t_{i-0})e^{-|i - (i-0))|/\Delta} &= x(t_i)e^{-0/\Delta} = x(t_i)\\
x(t_{i+1})e^{-|i - (i+1))|/\Delta} &= x(t_{i+1})e^{-1/\Delta}\\
x(t_{i+2})e^{-|i - (i-2))|/\Delta} &= x(t_{i+2})e^{-2/\Delta}
\end{align}$$
and so the sum is very much a convolution of the  sequence $x$
whose $i$-th term is $x(t_i)$
with a length-$5$ sequence $h$ whose $n$-th term is $e^{-|n|/\Delta}$
for $|n| \leq 2$ in this case (and a length-$(2D+1)$ sequence
 with $|n| \leq D$ more generally).

What you, the OP, need to understand is the following (and don't
get bogged down in words like convolution).  You are doing a
smoothing of the data that you have collected.  Each $x(t_i)$ of the
original data is being replaced by a $\tilde{x}(t_i)$ which is formed
as a weighted sum of $x(t_i)$ and its immediate neighbors.  The
weights assigned to $x(t_i)$ is the largest, the weight assigned to
its nearest neighbors $x(t_{i+1})$ and $x(t_{i-1})$ are smaller (but equal),
the weights assigned to slightly more distant neighbors $x(t_{i+2})$ and $x(t_{i-2})$ is yet smaller, and so on, till the we get to $x(t_{i+D})$ and $x(t_{i-D})$ which get the smallest weights of all. Anything more
than $D$ away gets a weight of $0$.  The weights in question are
normalized in the sense that the total weight is $1$, but let
us look at the raw weights which are alleged to be those $e^{-|i-k|/\Delta}$.
What you really have (before normalization of the weights) is
that sum (let us use $D = 3$ for simplicity) which is

$x(t_i) + \left.\left.e^{-1/\Delta}\right[x(t_{i+1}) + x(t_{i-1)}\right]
+ \left.\left.e^{-2/\Delta}\right[x(t_{i+2}) + x(t_{i-2})\right]
+ \left.\left.e^{-3/\Delta}\right[x(t_{i+3}) + x(t_{i-3})\right]$

$= e^{-3/\Delta}x(t_{i-3}) + e^{-2/\Delta}x(t_{i-2}) +
e^{-1/\Delta}x(t_{i-1}) + x(t_i) + e^{-1/\Delta}x(t_{i+1}) + e^{-2/\Delta}x(t_{i+2})
+ e^{-3/\Delta}x(t_{i+3})\quad (1)$
where on the right side we have the elements occurring in the
natural order, flanking the central element $x(t_i)$ on both sides.
Now, the sum of the raw weights exhibited in $(1)$ is
$$e^{-3/\Delta} + e^{-2/\Delta} +
e^{-1/\Delta} + 1 + e^{-1/\Delta} + e^{-2/\Delta}
+ e^{-3/\Delta} \tag{2}$$
which is just the element $Z$ in the question, and so the
normalized weights are just the ones exhibited in $(2)$
divided by $Z$. A better way of writing $Z$ might be
$$Z = \sum_{k=-D}^D e^{-|k|/\Delta} = 1 + 2\sum_{k=1}^De^{-k/\Delta}
\tag{3}$$
and perhaps a better way of expressing the first sum might be
$$\tilde{x}(t_i) = \frac{1}{Z}\sum_{k=-D}^D x(t_{i-k})e^{-|k|/\Delta}$$
which expresses more clearly that we are forming a weighted sum
of $x(t_i)$ and its nearest neighbors and that the weights assigned
to neighbors at distance $k$ decrease exponentially as a function
of $k$ (here $k \geq 0$ of course since we are not distinguishing
between the directions of separation).


Can two modulations be used in the same frame Can two modulations be used in the same frame?
For qam/qpsk/psk/ modulation scheme for the data bits can i have a fixed preamble with psk modulation?
does this help in achieving faster carrier freq lock and timing lock?

Can two modulations be used in the same frame.

Yes. Nobody's stopping you....

  For qam/qpsk/psk/ modulation scheme for the data bits can i have a fixed preamble with psk modulation?

If you want, why not? It's actually pretty common that systems with a flexible modulation scheme (e.g. IEEE 802.11, "wlan") have a header symbol that's in a fixed modulation that tells the receiver what the rest of the frame is.


CMV: COLA constraint is not necessary In the literature on the short-time Fourier transform (STFT) it is often stated or implied that the constant overlap-add (COLA) constraint must be satisfied in order for the STFT to be invertible.

Example: https://ccrma.stanford.edu/~jos/sasp/FBS_Perfect_Reconstruction.html

People seem go to a lot of trouble to choose a window/overlap combination that satisfies COLA.

I'm challenging the community here to "change my view" that the COLA constraint is sufficient but not necessary. I contend that you can achieve perfect invertibility under much weaker conditions and so have access to a very broad range of windows and overlaps.

I've written this blog post to examine this question. I would love to know if I am missing something.

This script provides an example of a window that is not COLA compliant, but which creates an invertible STFT.

You just triggered one  old itch, so this is a partial (both meanings) answer only. For discrete-time, uniformly sampled STFT, one often uses the usual fixed-window with $3/4$, $1/2$, $1/4$ overlaps, and somehow uses  direct closed form inverses. However, as long as there is redundancy (and completeness), there exists a infinite option for inverses.

In a   collaborative work (with Jean-Christophe Pesquet and Jérôme Gauthier), following other works on vector-frames, we used the formalism of complex oversampled multirate filter banks

to model analysis STFT, with the aims of

checking conditions under which the STFT were invertible, given a hop between time-frames,
building "better" inverses (better localized in time or frequency) with optimization.

Our aim at that time was double:

for data in higher dimensions $d$, reduce the STFT redundancy $r^d$ as far as we could,
-for 1D data, reduce the number of channels and increase the hop as we could as similar quality for real-time event detection.

All was based on polyphase matices, as we did consider windowed transform only as a special case;  so the main outcomes were the following:

given any redundant analysis filter bank (complex or not, windowed-DFT like or not), check whether it is invertible
in most cases, anything redundant is almost surely inversible,
windows with negative values can be used, liked Kaiser windows, see On the Optimization of Oversampled DFT Filter Banks, 2007

you can reduce the degrees of freedom in the oversampled synthesis to perform better concentrated inverses after selection/denoising for (close-to) real-time computations,
ensure the Hermitian character in the inverse for scalar threshold.

There were remaining issues, among which: the difficulty to go back to the orthogonal setting when the redundancy tends to one. Weirds behaviors were observed, possibly related to Farey sequences, but I have no definitive clue. And
References can be found in Optimization of Synthesis Oversampled Complex Filter Banks, 2009. Some related codes is available at SURE-LET Optimal oversampled Complex Filter Banks synthesis toolbox. We do plan a better release, but it has been remaining upcoming for months.

However, we did not investigate the non-COLA effect, and the necessary condition is not stated. My belief is that it is not necessary (due to redundancy), but this should be proven.


"Sinusoidal" in DSP terminology This is probably a bad question but I'll ask anyway because I'm often confused by it. Many textbooks explain how the movement of $z$ around the unit circle can be broken up into a sine wave along the horizontal axis axis and a cosine wave on the vertical axis. Both of these waves are termed sinusoidal.

So, if both of these are termed sinusoidal, then is the movement around the unit circle also termed "sinusoidal" or is it a combination of two sinusoids?

in my opinion, either of these two functions of time can be rightly labeled "sinusoidal":

$$ x(t) = A \sin(\omega t + \theta) $$

$$ y(t) = A e^{j(\omega t + \theta)} $$

where $A$, $\omega$ can be any non-zero real number (i wouldn't call DC a "sinusoid") and $\theta$ can be any real number.


Normalise data which has peaks Some methods

#1

YData = YData / max(abs(YData));

#2

maxYdataAbs = abs(max(YData));
minYdataAbs = abs(min(YData));
YData = YData / (0.5*(maxYdataAbs + minYdataAbs));

I am not happy with neither one.

Assume that you have the following data.

Only one or two sharp peaks.
I am not happy with #1 or #2.

I am not sure anymore if the best way is to normalise to the range $[-1,1]$.

How should you normalise such a data?

You could normalize the average amplitude, i.e.

YData = YData / mean(abs(YData));

Or you could normalize the signal power to one, i.e.

YData = YData / sqrt(mean(abs(YData).^2));

If just the peaks are bothering you, you could use dynamic range compression, but that would introduce nonlinear distortions. As Phonon hinted, please tell us why you are not happy with the results and why you want to normalize your data in the first place. That would allow for a more precise answer.


How to shift the center frequency of an IIR filter I have a pre-designed in-series 14th-order IIR filter in second-order sections. Problem is, when I simulate and measure the behavior in my codec, there is a frequency shift downward (as compared to the MATLAB prediction). Is there a simple way to alter the coefficients in such a way to shift the center frequencies of each band? Thank you.

The problem was an incorrect sampling rate issue, as a lower sampling rate than the designed rate caused the center frequency shift downwards. Today I learned something about the relationship between sampling rate and IIR coefficients (there is a proportional relationship between sampling rate and center frequency).


Step response given impulse response $h = e^{-|t|}$ I was given the above homework problem and I know how the step response is achieved, but I'm not sure how to set up the integrals in order to evaluate the problem. I've attached my work, but it is wrong because I keep getting the wrong answer. The solution is given as

$$S(t) = u(-t)e^t + u(t)(2-e^{-t})$$

*Note: We have not yet covered laplace transform. In my class, we are expressing the step response as running sums/integrals of the impulse response.

HINT: In your work you implicitly assume that $t\ge 0$, so your solution is (probably) correct for $t\ge 0$, but you must also consider the case $t<0$.


Proof for the solution of homogenous difference equation Suppose the output $y_h[n]$ of a linear, time-invariant system is described by by the following equation for input $x[n]=0$,

$$\sum_{k=0}^{N}a_k y_h[n-k] = 0$$

My book states the $y_h[n]$ is in fact a member of a family of solution of the form:

$$y_h[n]=\sum_{m=1}^{N}A_mz_m^n$$

where the complex numbers $z$ are the zeros of $A(z)=\sum_{k=0}^Na_kz^{-k}$.

What is the proof that $y_h$ is in fact of the form given?

A simple substitution would result in:

$$\sum_{k=0}^N (a_k \sum_{m=1}^NA_mz_m^n) = 0$$

but I don't see to proof showing a set $\{z_1,z_2,...,z_N\}$ and $\{A_1, A_2,..., A_N\}$ indeed exist for all n. Sorry if this seems like a poor question.

If you substitute your second equation back into the first equation, you obtain

$$\begin{align}\sum_{k=0}^{N}a_k\sum_{m=1}^NA_mz_m^{n-k}&=\sum_{m=1}^NA_mz_m^n\sum_{k=0}^Na_kz_m^{-k}\end{align}\tag{1}$$

If $z_m$ are zeros of the polynomial

$$A(z)=\sum_{k=0}^Na_kz^{-k}\tag{2}$$

then Eq. $(1)$ equals zero, as required for the homogenous solution.


texture matching between patches of an image I have an image A which I have divided into 4 x 4 subband images. For a given patch P1 in image A, I need to find another patch P2  in the same image A which is most similar: have same texture.

Till now my appraoch is to create feature vectors Fi for all i patches and using Euclid's distance formula find out which patch's feature vector is closest to given patch P1's feature vector.

Currently I have added following features:

  Mean & Standard Deviation of brightness (using L channel of LAB colorspace)
  Mean & Standard Deviation of color values in A channel of LAB colorspace
  Mean & Standard Deviation of color values in B channel of LAB colorspace

Although I am getting similar patches but I still think the matching can be improved if I incorporate more prominent features (which I'm unaware of).

Following are my queries:

  I have a doubt whether this is a good way to compare color between two patches.
  Please suggest some more good features which can be helpful to get a proper differentiation.
  Although most of my images are proper but are there any restrictions/drawbacks with the approach if image is dark and/or noisy? If yes, is there a good alternative approach or
  feature(s)?

Please help!

The mean and standard deviation are two measurements of a distribution. Others you could also use are higher order moments like 'skewness' (how skewed the distribution is) and 'kurtosis' (how 'peaky' the distribution is).

However, what I would try is a histogram of the values for each of the channels. For example, if you used a histogram with 16 bins, you would get a feature vector with 48 values. The histogram bin values can be thought of as a descriptor of the distribution.

To account for changes in brightness, you need to introduce some illumination invariance. Already the LAB colour-space splits the image into luminance (L) and color (AB) channels. To achieve illumination invariance I would normalise the L channel histogram values using the norm of whatever measure you are using (in this case, euclidean distance, it would be the l2 norm). Normalising the A or B channels is probably not necessary. The normalisation accounts for any linear changes in luminance.

Depending on the images, you could also add vectors like:

histogram of oriented gradients (HOG)
histogram of connected components (make B/W then count how many pixels are white in the 8 connected neighbours)
k-means clustering of LAB values to get k centre LAB values which you can compare to other images.

and so on.


Extract not blurry images from video I need a tool or code for extract not blurry images from video (it can be used for panorama or SFM building from video)

something like
http://research.microsoft.com/en-us/um/redmond/groups/ivm/sharppanoramas/
(It is worth noting that software for building panorama from video is already exist, but I don't know why they don't use it for SFM)

I want to get images from video and put them to vsfm

For example virtualdub gives blurry images when I try to convert video to images.

One straightforward approach is to do blur detection (see here and here) and prune the frames that are detected.


Two generated tones creates a side tone - how to avoid the side tone I'm trying to create two tones at different frequencies. So I have one tone at 15000 Hz and one at 15400 Hz using Audacity. The problem i'm running into is there is a third tone created at a much lower frequency I can tell exists from both my ears and a spectrogram.

How can I get rid of the lower tone? Where does it come from? I only have two tracks

I have tried the following:

Launch Audacity.
Generate a 15000 Hz tone in the track created by default.
Add a new track.
Generate a 15400 Hz tone in the new track.

A lower frequency tone appears during playback.

The reason is that both tracks have high levels, so their sum exceeds 1.0; and Audacity applies clipping or limiting. This non-linear operation is what creates the overtone.

Simply decrease the level of both tracks to -6dB (or use an amplitude of 0.5 during signal generation) to make the overtone disappear.


adjust mean of signal using exponential I have discrete signals whose values are between 0 and 1.

I wish to post-process such a signal such that the mean of its values equals 0.5, yet keeping maximum and minimum values 0/1 intact.
Intuitively, I think an 'element-wise' exponential would be the right approach, yet I can't seem to come up with the method to determine the value of the appropriate exponential.

Fictive matlab example to illustrate:

>> s=[0.2 0.7 1 0.5 0.4 0.2 0 0.3];
>> mu=mean(s)
mu =  0.41250
>> s_adjusted=s.^0.69998;
>> mu_adjusted=mean(s_adjusted)
mu_adjusted =  0.50000

What is the best method to determine the 0.69998 exponential ?

You need to do it numerically, as a related toy problem $2^a + 3^a = 4$ appears to have no symbolic solution. Bisection method (binary search) is probably the easiest solver to implement:

minexponent = 0;
maxexponent = 2;
targetsum = 10;
precision = 32;
loop precision times {
  exponent = (minexponent + maxexponent)*0.5;
  if (sum(s.^exponent) < targetsum) {
    maxexponent = exponent;
  } else {
    minexponent = exponent;
  }
}
return exponent;

The smallest possible exponent value is 0 for your problem. If you are unsure about the largest possible exponent value, you can test the easy-to-check exponents 1, 2, 4, 8, 16 etc. (as many as needed before the sum goes below the target) before the above code, exponentiating into an auxiliary array by multiplying each element by itself.

Note that if all of the elements are zero, or if all of the elements are one, then you cannot get any other mean than that by exponentiation.

The plot below shows 1000 runs of the algorithm, and the value of exponent at each iteration.


Direct Digital Synthesis - Super Sampling - Unknown Harmonics I'm using a Xilinx FPGA (Virtex) with 4 DDS cores (each supplied a 250MHz clock) used in parallel to provide a samples to a DAC38J82IAAV from TI, 16 bit DAC running at 1 Gsps. The four cores super sample (interleave their samples) to provide samples at 1 Gsps (250MHz clk * 4 cores = 1 Gsps) for the DAC. I'm trying to output frequencies between 20MHz and 500MHz with the DAC.

The DDS cores are configured in rasterized mode with modulus of 10000. I set the Phase Increment (PINC) to be (desired_freq * 10000)/250 and the Phase offset for each core to be 0, PINC/4, PINC/2, 3*PINC/4.

This works great for generating frequencies above 250MHz, but for frequencies below 250MHz, I see harmonics at Fout, 2*Fout, 3*Fout and so on. I would expect to see harmonics at Fout, Fsample-Fout, Fsample+Fout and so on, which I do see for frequencies above 250MHz.

Why am I seeing multiple harmonics at Fout interval for frequencies below 250MHz?

I think that with a 250 MHz sampling clock, you should be able to generate a signal with a 125 MHz max frequency.
Also, what type of DDS do you use? Is it based on a table of sine/cosine or is it a CORDIC algorithm? It has an impact on precision.
Is the signal path from DDS to DAC ok? Do they work at the same sampling frequency?
If not, a resampling module should be used or correctly configured.


Convolution and Cross Correlation on 2D Image I have read that convolution and cross-correlation are the same thing, but convolution flips 180 degrees (images), or time reverses (sequences) the kernel, before performing the elementwise multiplication of the input and the kernel.

I understand how this works on 1D signals, however I do not know how to perform these two operations for 2D images, and it's been really hard to find any good example on the web.

I would very much appreciate an example, showing both Convolution and Cross-correlation of an image. Clearly specifying what the indexes stand for and how the output image is calculated.

For example (feel free to use any other example if you think it would be better):

If as you said you understand well the 1-D convolution/cross-correlation functioning (the Wikipedia first graph explains it in a clear way), the 2-D version is very similar!

This website explains 2-D convolution in a simple way with clear indicies and examples. In a nutshell, the kernel has to be flipped for convolution and this means your kernel example would become
[M L; K J] (sorry for the formatting). Then you multiply elementwise and sum the results, as you would do in the 1-D case.

Cross-correlation functioning is the same without flipping the kernel.


Is FIR decimation the fastest possible method for integer downsampling? For a digital music synthesizer I am developing, I have a stream of incoming samples $x_i$ and would like to downsample by $M$ to obtain a stream of samples $y_{Mi}$. I am currently bandpassing the signal with an FIR which is a brick-wall filter windowed by Blackman-Harris, and it sounds great to my ears. But I have to compute $y_{Mi} = \sum_{j=0}^{N-1} x_{Mi-j} h_j$ for every outgoing sample. $N$ is usually 64 to 256, and $M$ is usually 8 or 16, so an FFT convolution doesn't seem worth it for performance.

If there a better method than naive FIR for integer decimation of audio at these scales? Can I make some reasonable tradeoffs in quality to achieve better performance? For example, polyphase decimation is a nice trick, but since RAM/cache is not an issue, I believe they should run in the same speed since they involve an equal number of multiplications.

Another way to design a FIR filter involves matching the impulse response to an analog prototype's impulse response.  The idea being that  an IIR can often be truncated to FIR. The FIR will not be symmetric  and not have linear phase but  the group delay can be much smaller.
I believe that Friedlander has a paper  on how to do the approximation. There was a Matlab Signal Processing Toolkit function that did this but I haven't had much cause to use it.
Comet also can do nonsymetric FIR designs.

I also suggest that you try using Parks Mclelan if it can give you a shorter filter than using the window method


Which Machine or deep learning algorithm is appropriate of this issue? Suppose I have $n$ features as $Y$ = ($y_1 , y_2 ...., y_n$), and a matrix of $J$ of dimension $M$x$N$, one feature of $Y$ is selected randomly to be convolved with one random column of $J$ resulting a new vector called, for example, $X$.

The question, if I have the matrix $J$ and the vector $X$, can I detect the selected feature from $Y$ and selected column from $J$ using any deep learning algorithm? such that CNN or DNN ...... etc.   What's about if I don't have the matrix $J$ am I still able to detect the selected feature $y$ and column of $J$ based on $X$ ?

For example, let's $Y = [1,2,3,4]$ and matrix $J$ is any random matrix of dimension $4$x$4$. Hence, one random number of $Y$ either 1 or 2 or 3 or 4 is going to be convolved  with one column of matrix $J$, let's say element 2  of vector $Y$ is convolved with with second column of matrix $J$, so the resulted vector is $X$. So the question can we estimate the element selected from $Y$ based on matrix $J$ and resultant vector $X$ using machine or deep learning algorithm? I think yes we can, but what's the most appropriate algorithm to do that?

thank you

Your question can be interpreted in many ways. First, convolution with a number is not a common operation. Convolution is defined by:

$$(f*g)=\int_{-\infty}^{\infty}f(\tau)g(t-\tau)d\tau$$

Assuming a constant signal $g(t)=a\in \mathbb{R}$ this is reduced to:
$$(f*g)=\int_{-\infty}^{\infty}af(\tau)d\tau=a\int_{-\infty}^{\infty}f(\tau)d\tau$$

I am assuming, therefore, that you meant multiplication and not convolution. Correct me otherwise and I will edit the answer accordingly.

Secondly, assume for your example the given parameters:

$$Y=[1,2,3,4],J=\left[
 \begin{matrix}
  1 & 2 & 3 & 4\\
  1 & 2 & 3 & 4\\
  0 & 0 & 0 & 0\\
  0 & 0 & 0 & 0
 \end{matrix}
\right], X=\left[\begin{matrix}4 \\ 4\\0\\0\end{matrix}\right]$$

It is easy to see that a couple of use-cases could give us this result. For example, it could be that the chosen $y$ was $y_4=1$ and the chosen column was the most right column $j_3$. On the other hand, the combinations $\{y_1,j_1\},\{y_3,j_0\}$ are also plausible. I have purposedly gave a very bad, linearly dependent, matrix $J$ to illustrate that the choice of $J$ in the design process or formulation of the problem is of major interest and is fundamental for answering you question. The following example illustrates the complete opposite:

$$Y=[1,2,3,4],J=\left[
 \begin{matrix}
  1 & 0 & 0 & 0\\
  0 & 1 & 0 & 0\\
  0 & 0 & 1 & 0\\
  0 & 0 & 0 & 1
 \end{matrix}
\right], X=\left[\begin{matrix}x_0 \\ x_1\\x_2\\x_3\end{matrix}\right]$$

Obviously, $x_0=y*j_0$,$x_1=y*j_1$, and so on. Here, it is easy to deduce $y$ and only a small number of $X$ are even possible.

Generally, if we call $J$ a dictionary, you are looking at something called Orthogonal Matching Pursuit (OMP) problem, which is not a machine or deep learning problem. Otherwise, you can try and use a deep learning model to try and learn the dictionary, though this would be unrecommended as such a solution is not compatible with such a problem in my opinion.

Hope this answers your question. If not, or I misinterpreted your question, feel free to comment and I will try to supply more information if I can.


Getting stuck into an open source image processing/computer vision project I have a degree in Electrical Engineering (focusing on VLSI). However, I am more mathematically minded, and only recently discovered that image processing seems to contain lots of interesting maths. I am looking to break into this field in a about a year or so in spite of not having any formal experience in this area.

I was wondering: would contributing to code development on one of the open source projects help my cause? Would this enable me to overcome my lack of formal education or experience in this field?

I would also appreciate it if someone could point me to good open source resources. I know of GitHub and Nasa's Vision workbench. However, I have not contributed to any open source projects and I don't know what to expect.

I believe it would if you intend to practice your knowledge, rather than go into academia.

I suggest you look into the code base of opencv, simplecv, and scikit-image, as suggested. Once you become comfortable with using them, you can submit patches for bugs and new algorithms.

Nevertheless, you will also need to understand what you are doing at a theoretical level, so get some books too. There are also courses online, such as this and that.

Good luck!