Slow sampling of signal and averaging over n periods Can one in theory achieve as high accuracy as needed by sampling a sinusoidal signal over n Periods? Is there a formula that connects the signal frequency (f), sampling frequency (fs), periods (P) and accuracy (Acc) or can someone explain how accuracy rise with number of periods? Let's say I want to know the signals amplitude with an error of less than 0.1%.

Lets say I want to sample a signal with a sampling frequency exactly four times that of the signal frequency. Let's say the first sample hit the sinusoidal signal at it's peak, the second at a zero-crossing, the third at it's negative peak and the fourth at a zero-crossing. This average will yield (1 + 0 + abs(-1) + 0)/4 = 1/2. Averaging over n periods would not change this average value in theory. Let's say the first value hits 45 degrees off the peak and the next 3 samples follows with 90 degree phase difference. The average would now be: (0,71 + abs(-0,71) + abs(-0,71) + 0,71)/4 = 0,71. Averaging over n periods would not change this average value either.

But I have been thinking.

If one introduces a small phase shift of let's say 1 degree per period based on the difference between the signal frequency and the sampling frequency, one could over 360 periods have sampled along the "whole" sinus with samples at 360 different points on a sinusoidal signal. This of course assumes the signal is periodically (and it is in this case). This would give a pretty good estimate of the sinusoidal signals amplitude, but how much is the error here? Is this something that can be calculated easily with some sort of formula or do I have to use MATLAB?

And what happens in "real" situations when the sampling frequency can't be exactly 1 degree off per period, but maybe 0,9384283437 and we have to sample e.g., 324986234 (just a big number to prove the point that an unrealistic number of periods have to be sampled) periods to be able to sample an integer number of periods. I understand that here there will be a trade-off between needed accuracy and number of periods needed to be sampled, but how can I calculate this?

If what I have written above makes sense and is correct, it seems doable to obtain the wanted accuracy by adjusting the phase shift and number of periods to sample over.

Other: I can understand how some would think of this long text as more of an explanation or help to calculate different things. And this is mostly true as I have no experience or academic background in digital signal processing and therefor need to get confirmation on my idea, thought and understanding of the subject. Things I ask here are things I have not found the answers to by any text book or by googling, even though they for sure are available a lot of places, none of them have been able to completely show me or make me understand this topic and I am therefor asking by writing this text above. Hopefully others may benefit from this questions and the answers provided. Thanks! =)

Accuracy is independent of the sampling frequency and frequency, as long as your signals frequency is smaller than half of the sampling frequency. In the noiseless case its also independent of the number of periods.

Your calculation of the mean of the sample values has no connection to calculating the amplitude of the signal. To calculate the amplitude one would normaly use a DFT.

Here is a sample Matlab code using your examples, and showing that both give the same amplitude:

x1 = abs(fft(cos([0:pi/2:3/2*pi])));
x2 = abs(fft(cos([0:pi/2:3/2*pi]+pi/4)));

Both signals show the same amplitude (which must be scaled properly to get the exact value).

In the case of noise, using multiple periods is kind of averaging, where the following rule holds (for uncorrelated noise): doubling the number of averages (periods) gives 3 dB more SNR (signal to noise ratio).


downsampling DFT with aliasing (Schuam's DSP Outline, 2nd edition, problem 6.11(c), page 241).

Is there a DFT down-sampling property that looks something like this:

Given $x[\![n]\!]_M$ we want to downsample from M to N to obtain $y[\![n]\!]_N$...

$$
\Bigg(y[\![n]\!]_N = x\Big[\!\Big[\  [\![n]\!]_M\  \Big]\!\Big]_N \Bigg) \xrightarrow{DFT} \left( Y[\![k[\!]_N = X\left[\!\left[ \frac{M}{N}k  \right]\!\right]_M \right)
$$

$$
\text{for } (M > N)\ \&\ (M \ modulus\  N = 0)
$$

For Example, if I downsample $x[\![n]\!]_6$ which was modulo 6 and convert it to modulo 3 and assign it to $y[\![n]\!]_3$ which is modulo 3, then it looks like this in frequency domain:

$$
\Bigg(y[\![n]\!]_3 = x\Big[\!\Big[\ [\![n]\!]_6\ \Big]\!\Big]_3\Bigg) \xrightarrow{DFT} (Y[\![ k]\!]_3 = X[\![ 2k]\!]_6)
$$

Is this relationship true, and is there an easier way to write this relation as a transform?

Maybe I can invent a math notation to keep my DFT table from getting to crazy:

$$
\Bigg(y[\![n]\!]_N = x[\![n]\!]_{M\ to\ N} \Bigg) \xrightarrow{DFT} \left( Y[\![k[\!]_N = X\left[\!\left[ \frac{M}{N}k  \right]\!\right]_M \right)
$$

First of all very similar questions have been asked here and even I have answered some of them, but now I will provide an answer using yet another approach.

Our aim is to show the effect of down-sampling on the DFT.
Let $x[n]$ be a sequence of length $N$. And let's define $y[n]$ by downsampling $x[n]$ by an integer factor $M$ as shown:

$$ x[n] \longrightarrow \boxed{ \downarrow M } \longrightarrow  y[n] = x[Mn]  \tag{1}$$

$y[n]$ is of length $K = N/M$. At the moment assume $K$ is integer.

I will make use of the relation between DTFT $X(\omega)$ and DFT $X[k]$ for finite length sequence which states:

$$ X[k] = X(\omega)|_{w = \frac{2\pi}{N}k } = X( \frac{2\pi}{N}k  ) ~~~,~~~ k = 0,1,2,...,N-1 \tag{2}$$

From the theoretical DTFT definition, it can be shown that the DTFT $Y(\omega)$ of $y[n]$ is given as:

$$ Y(\omega) = \frac{1}{M} \sum_{m=0}^{M-1} X( \frac{ \omega + 2\pi k}{M} ) ~~~,~~~ -\pi \leq \omega < \pi \tag{3}$$

And we define the $K$-point DFT of $Y[k]$ as:

$$ Y[k] =  Y(\omega)|_{\omega = \frac{2\pi}{K} k } = Y(\frac{2\pi}{K} k) ~~~,~~~ k = 0,1,...,K-1  \tag{4}$$

Note the range of DFT index $k$ for $Y[k]$. Since $y[n]$ is a $K$-point sequence we have defined a $K$-point DFT of it.

Finally plug Eq(4) into Eq(3)

$$ Y[k] =  Y(\frac{2\pi}{K} k) = \frac{1}{M} \sum_{m=0}^{M-1} X( \frac{ \frac{2\pi}{K} k + 2\pi m}{M} ) $$

$$ Y[k] =  Y(\frac{2\pi}{K} k) = \frac{1}{M} \sum_{m=0}^{M-1} X( \frac{2\pi}{KM} k + \frac{2\pi}{M} m ) $$

Now $KM = N$ and we replace $\frac{2\pi}{M}$ with $\frac{2\pi}{N}(N/M) = \frac{2\pi}{N}K $ to get

$$ Y[k] =  Y(\frac{2\pi}{K} k) = \frac{1}{M} \sum_{m=0}^{M-1} X( \frac{2\pi}{N} k + \frac{2\pi}{N} Km ) $$

$$ Y[k] =  Y(\frac{2\pi}{K} k) = \frac{1}{M} \sum_{m=0}^{M-1} X( \frac{2\pi}{N} ( k + Km) ) $$

And finally we recognize that $X( \frac{2\pi}{N} ( k + Km) )$ is actually the $N$-point DFT $X[r]$ of $x[n]$, evaluated at $r = k + Km$ and insert it into the equation:

$$ \boxed{ Y[k] = \frac{1}{M} \sum_{m=0}^{M-1} X[ k + Km] ~~~,~~~k = 0,1,...,K-1 } \tag{5}$$

Note that DFT index k for $Y[k]$ is of modulo-K where as DFT index $r = k + Km$ for $X[r]$ is of modulo-N. Hence Eq(5) can alo be written like:

$$ \boxed{ Y[(k)_K] = \frac{1}{M} \sum_{m=0}^{M-1} X[ (k + Km)_N] } \tag{6}$$

The following is a MATLAB/ OCTAVE code to demonstrate the result:

N = 32;        % Length of x[n]
M = 4;         % Downsampling factor (integer)
K = N/M;       % length of y[n] (assuming N divisble by M)

x = randn(1,N);     % x[n]
y = x(1:M:end);     % y[n] = x[Mn]

X = fft(x,N);       % N-point DFT X[r] of x[n]

Y = zeros(1,length(y));  % K-point DFT Y[k] of y[n]
k = 0:length(Y)-1;       % DFT index range k

for m = 0:M-1            % implements Eq.6 above
    Y = Y +  X( mod(k+K*m,N)+1 );
end
Y = Y/M;

Y2 = fft(y,K);          % Comput Y[k] directly from y[n] for         checking.

figure,stem(abs(Y))     % display that they are the same..
hold on
stem(abs(Y2),'r+');


Determining height of object from images Given two images of an object taken from the same distance but different vertical angles, is there a way to determine the physical height of the object?

If you are only interested in the measurement of a single dimension (such as the height or the width or radius) of an object then you do not need multiple images. A single image would have worked if you could determine some scale factors which would depend on the distance of the camera to the object and orientation angles (which would be defined in a scaling matrices) etc.

If the measured dimension was inclined in 3D wrt the camera, then you should also apply projective geometry corrections which would also be defined in respective matrices.

Two images, on the other hand, would let you create the complete 3D environment through which any measurement is easily accomplished by their respective 3D point coordinates inside the world coordinate system. However it's considerably more complex and sensitive to create a precise definition of 3D points from just two images due to several reasons. You need very strict control on the imaging setup.


Is there an equivalent time-domain operation to the spectrum averaging in Welch's method? In Welch's method a time-series is segmented and then the subsequent FFTs of the segments are averaged. Is there an equivalent time-domain operation?

I guess doing the the Welch estimate and producing an impulse response using an FIR filter might do what I want, but I'm looking for a way to do it in the time-domain directly (getting rid of noise, keeping signal properties)...

Yes, you can average in time lag (autocorrelation) as well as in frequency. You can also combine the two domains, as long as it is reasonable to treat the data as being stationary.

A. H. Nuttall and G. C. Carter, "Spectral estimation using combined time and lag weighting," in Proceedings of the IEEE, vol. 70, no. 9, pp. 1115-1125, Sept. 1982.


About voice filtering I have a waveform containing

voice
some children screaming
some times when there is no noise at all
some times there are cars driving by when talking

I want to detect the start and ending times of the loudest voice or the areas in which the same speaker is talking to cut them out.

What I did was a FFT on the waveform to realize that voice means there is(in a small time unit of a few ms) always a high amplitude at the base frequency to be followed by equidistant frequencies, being apart in the integral multiple of the base frequency.

But this also applied to some noises.

The goal is to detect voice and cut the parts containing voice out(ie. keep them) with all background sounds. It is not about only having voice but finding out where voice is and cutting it out with all cars children or whatever was surrounding the microphone while recording. The problem is that I get many false positives when music(no voice) is being played or some squealing sound can be heard.

What is a better approach to filter out the times voice is at in a wav file with these given features?

I don't think filtering is the way to go. The spectra of the different signals (the desired voice, the screamy children and the cars) are probably overlapping, so there is no real way to get rid of all that noise without spoiling your signal.

Actually, the approach would depend in how loud the noise is in comparison to the signal. If the speaker is near the microphone, and therefore the relative noise level is ''low'', this sounds like a VAD kind of problem.

If noise level is quite low, an energy-based VAD may work: after windowing the input signal (let's say to 20ms windows), each window's energy can be computed. Then, a threshold can be applied in order to decide if that window is noise or voice. You would simply discard the windows that aren't load enough to meet the threshold. This is unarguably the simplest approach.

A simple improvement to this kind of VAD, assuming you only have one speaker, is to reduce the bandwidth in which you apply the stated algorithm. This way, if the kids voice spectrum begins at about e.g. 500Hz, and the cars induce some low frequencies up to e.g. 80Hz, you could use a band pass (e.g. 100Hz - 300Hz, assuming the speaker fundamental frequency is about 200Hz) filtered signal to decide wether each window is or isn't noise, and select those who aren't in the original signal. I know this is not ideal, as too many coincidences are needed to have a reliable system.

Unfortunately, this won't work when noise gets louder. Furthermore, some phonems may be discarded due to their low energy, such as /s/ or /z/. This just makes speech sound weird. There are some more advanced VAD techniques such as Sohn's A Statistical Model-Based Voice Activity Detection. You can find some Matlab implementations in the net. I have used this in addition to spectral substraction noise reduction (which doesn't seem to be needed in your case) in a noisy environment (babble + WGN) and results were susprisingly good.

Matlab code for using Voicebox's Sohn VAD.

in_file='yourfile.wav';
[s,fs,wmode,fidx]=readwav(in_file,'p',-1,3500); %3500 is the number of samples to discard at the beginning of your file.
[y1,zo]=vadsohn(s,fs); %y1 is 1 when there is voice, 0 when there isn't
y=s(1:length(y1)).*y1; %apply mask
y(1)=0;

for i=2:length(y)
    if y(i-1)<y(i)
        StartTime=(i-2)/fs
    elseif y(i-1)>y(i)
        EndTime=(i-2)/fs
    end
end


Relaxation time in a stationary stochastic process Foreword: I'll give some background here because I'm not sure all my premises are correct. I encourage you to IGNORE THE BACKGROUND, unless you're really interested / bored :)

The Question

I want to calculate the eigenvalues of this system.

$\frac{dM_0}{dt}=k+t_0 M_1 - M_0 (t_1+d)$

$\frac{dM_1}{dt}=t_1 M_0 - M_1 (t_0+d)$

In matrix form, I want the system to look like $m^{'}=Ax$, and I'd compute the eigenvalues from the square matrix A. But here I have $m^{'}=Ax+b$ -an inhomogeneous set of ODEs.

My intuition is that eigenvalues don't depend on $k$ ($b$) at all, but I have been unable to confirm this by any means (consulted many ODE/eigenanalysis texts, and the internet). So the question is: can I ignore $k$ when computing the eigenvalues?

Background

In order to construct the steady state probability distribution function (PDF), I run a simulation many times, record the "final" state of the system, and repeat. At the end of N simulations, I have a vector of final states $[X_1 X_2 ... X_N]$, which as $N\rightarrow \infty$ the histogram of which converges to the probability distribution function.

It seems to me that the optimal time $T^*$ to simulate each individual run is when the stochastic system has "relaxed" enough that it no longer "remembers" its starting condition. A good metric for this is autocorrelation: the lag time for which autocorrelation goes to zero (or some arbitrarily small value) seems like the time I should choose for $T^*$. If I choose a time $T<T^*$, every simulation is correlated with the start condition, and not independent observations. Thus, the PDF will never converge on the correct distribution. If I choose a time $T>T^*$ then I do "extra" computation. Since it turns out the simulation is VERY inefficient, even simulating a little past $T^*$ results in significant computational overhead.

The question is, how to choose $T^*$? Suffice to say, numerically computing it is out of the question. One idea suggested to me is to calculate the eigenvalues of the system. The inverse of an eigenvalue corresponds to the relaxation time of a system, which is equivalent to the lag-time tau at which autocorrelation = $\frac{1}{e}$. Thus, the time $T^*$ could be chosen instantaneously by selecting the longest relaxation time of the system (smallest eigenvalue)--or some multiple of this time to get the correlation sufficiently small. Bonus points if you weigh in on whether this is a good idea or bad, and why--but not necessary to answer the question.

Considering the complexity of the question, and the rating of your reputation, i assume you are still alive.

In the following, i am assuming $k$ is a constant, since it is not a variable -fixed or random- expressed in the question.

If you are still interested, which is a numerically strong assumption, the system can be rewritten trivially as:

$$\frac{d}{dt}
\left(\begin{array}{cc}
x_1  \\
x_2
\end{array}\right)=
\left(\begin{array}{cc}
-(t_1+d) & t_0 \\
t_1 & -(t_0+d)
\end{array}\right)
\left(\begin{array}{cc}
x_1  \\
x_2
\end{array}\right)
+
\left(\begin{array}{cc}
k  \\
0
\end{array}\right)
=
\left(-d
\left(\begin{array}{cc}
1 & 0  \\
0 & 1
\end{array}\right)
+
\left(\begin{array}{cc}
-1 & 1  \\
1 & -1
\end{array}\right)
\left(\begin{array}{cc}
t_1  \\
t_0
\end{array}\right)\right)
\left(\begin{array}{cc}
x_1  \\
x_2
\end{array}\right)
+
\left(\begin{array}{cc}
k  \\
0
\end{array}\right)$$

Which has the standard form:
$$Dx=Ax+x_0$$

The constant vector $x_0=[k;0]$ or equivalently the constant vector $A^{-1}x_0$ is a state reference.

Through the variable change:
$$x+x_0 \to x$$
you can safely discard the $x_0$ vector.

This is an autonomous system which, unless you forget to include perturbations or inputs, is a very boring system. I am sure you did.


Doppler Shift and the Pulse Repetition Interval Suppose that I have a signal $x(t)$ consisting of $N$ pulses at a given frequency $\omega_{c}$. I space them $T$ seconds apart, ie the Pulse Repetition Interval (PRI) for the waveform is $T$.

Now this signal bounces off a target and experiences a large Doppler shift, such that the compression of the pulses are evident. My question is -- does the PRI change as well?

In other words, are each of the waveforms still $T$ seconds apart, or does the PRI scale with the Doppler shift?

I guess the question really is, should the model of the second pulse be:

$$ x(t-T) \underset{v}{\rightarrow} x(\frac{t -T - \tau}{v}) $$

or

$$ x(t-T) \underset{v}{\rightarrow} x(\frac{t -T}{v} - \tau) $$

or

$$ x(t-T) \underset{v}{\rightarrow} x(\frac{t}{v} - T - \tau) $$

Where $\tau$ is the target delay and $\frac{1}{v}$ is the Doppler shift.

Note: A similar question is asked in the DSP Stackexchange question
Doppler Shift of Phase Coded Pulse Compression Waveform

Actually the PRI also changes. There is a complete explanation about this topic in "Cobbold - Foundations of Biomedical Ultrasound" chapter 10 (Pulsed Methods for Flow Velocity Estimation and Imaging).

Consider a case where PRI of emitted pulse is 1 s and the wave speed is 9 m/s and the object is 10 m away and moving with speed of 1 m/s toward the antenna. The first echo will be received 2 seconds after the emission (after 1 seconds, the object moved 1 meter toward the antenna and placed 9 meter away from the antenna and the transmitted wave traveled 9 meter away from the antenna). For the next pulse the echo will be received 1.8 seconds after the emission (considering the object is moved closer to the antenna during the PRI between first and second pulse). So the first echo received at t=2 and the second echo received at t=2.8 and the PRI would be 0.8 s (the same is true for the next pulses until the object passes over the antenna).

Also the transformation you have here is something like x(t - (x0-vt)/c ).


How to make a signal conjugate symmetric? Take the simple frequency-domain band-pass filtering operation below . . .

NFFT = 128;
x = randn(NFFT,1);
H = zeros(NFFT,1);
H(10:20) = 1;
y = ifft(H.*fft(x), 'symmetric');

This gives a real output because I'm using the conjugate symmetric flag to the ifft operation.

I want a function that returns the conjugate symmetric version of H, so I do not have to rely on the builtin symmetric option in Matlab's ifft. NFFT can be any positive integer. This could be called something like this . .

H(10:20) = 1;
H = MakeConjSym(H);

Conjugate symmetric means

$$f(-x) = f^{\ast}(x)$$

i.e. the sign of the imaginary part is opposite when $x<0$

The FFT of a real signal is conjugate symmetric. One half of the spectrum is the positive frequencies, and the other half is the negative. The negative coefficients are conjugate of the positive.

So if you do filtering, your envelope must do both the positive frequencies and their corresponding negative frequencies, so that the imaginary bits cancel out.

In your example, H does only one half. That is why the output has imaginary bits in it. What you want is

NFFT = 128;
x = randn(NFFT,1);
H = zeros(NFFT,1);
H(10:20) = 1;
H(end-20+2:end-10+2) = 1;    % Other half
y = ifft(H.*fft(x));


Convert low pass continuous time filter design to bandpass, discrete time I am trying to convert the continuous time transfer function of a second order lowpass Butterworth filter is given by:

To a bandpass fourth order bandpass digital filter, I first apply the mapping to bandpass:

We therefore obtain the bandpass continuous time transfer function

In order to convert this to a digital filter with with band edge frequencies of 2000Hz and 2819.3 Hz, and a sampling rate of 8 kHz. We first normalise and apply bilinear prewarping

approximating slightly we get

Then applying the bilinear transform

We get the fourth order discrete time filter

However
Checking the results at the normalised angular frequencies 0 and at the bandpass center frequency 3, I find that the frequency response is not as desired. is this an arithmetic mistake or have I missed something, your help would be very much appreciated.

There is an error in the analog bandpass filter transfer function. If you apply the LP-BP transformation as you stated it you should get the following transfer function:

$$H(s)=\frac{s^2(\Omega_2-\Omega_1)^2}{s^4+\sqrt{2}(\Omega_2-\Omega_1)s^3+(\Omega_1^2+\Omega_2^2)s^2+\sqrt{2}\Omega_1\Omega_2(\Omega_2-\Omega_1)s+\Omega_1^2\Omega_2^2}$$

With the given band edge frequencies (note that the lower band edge of the discrete-time filter must be $2000\,\text{Hz}$, not $4000\,\text{Hz}$) you correctly got $\Omega_1=2$ and $\Omega_2=4$. Applying the bilinear transform to the analog bandpass transfer function with these band edges finally gives the following discrete-time transfer function:

$$H(z)=\frac{z^4-2z^2+1}{(10+3\sqrt{2})z^4+(12+2\sqrt{2})z^3+20z^2+(12-2\sqrt{2})z+10-3\sqrt{2}}$$


Distortion of a signal consisting of two periodic components We have a signal
$$x(t)=U_1\cos\omega_1t+U_2\cos\omega_2t$$
whose frequencies are $f_1=100\,\text{Hz}$ and $f_2=600\,\text{Hz}$. This signal is being discretisated in time using unipolar rectangular pulse with 1 V  amplitude, length $\tau$ and period $T=\frac{1}{f_0}$. Can we reconstruct this signal without distortion if we use ideal low-pass filter if sampling frequency is:
a) $f_0=1\,\text{kHz}$
b) $f_0=1.3\,\text{kHz}$

My attempt:
when we want to tell whether we can or cannot reconstruct original signal when sampling with specific frequency we compare sampling frequency with the frequency of the given signal, sampling frequency must be at least as twice big as signal frequency (NYQUIST criterion). However i cannot determine frequency of this signal since it has two components.
Since we are doing time discretisation with rectangular pulse we want to find fourier transform of $f(t)=x(t)s(t)$ where $s(t)=\sum_{n=-\infty}^{\infty} \frac{\tau}{T}sinc(n\omega_0\tau/2)e^{jn\omega_0t}$.
I know that multiplication in time is convolution in frequency but i don't know how to implement it in this case since
$X(j\omega)=\frac{U_1}{2}2\pi\delta(\omega -\omega_1)+\frac{U_1}{2}2\pi\delta(\omega +\omega_1) + \frac{U_2}{2}2\pi\delta(\omega -\omega_2) + \frac{U_2}{2}2\pi\delta(\omega +\omega_2)$
How can i do that and how can i tell, when i find fourier transform, if signal can be reconstructed or not?
Any help appreciated!

The question is trying to make you understand Aliasing

The low-pass filter is introduced (as a hint) to imply that the output is band-limited to half the Sampling frequency (Why ? See below for aliases).
This would allow you to understand how out of band signals can sneak into your sampled signals if left unchecked. Look at p5 here

You should read in detail about sampling.
Skipping the real math, as rule of thumb you can do as below to answer your questions quick

  a) for $f_0=1\,\text{kHz}$ the signal cannot be reconstructed

Once you sample your signal with a rate you create aliases at frequencies
$\pm f_{signal}\pm n*f_{sampling}$
where n can take any positive value

If you calculate all combinations for all your frequencies you end up with
[100 400 600 900 1100 1600 1900 2100 2900 ...]

Since you'd sampled with 1000Hz, and are using an (ideal) LPF at 500Hz, you'd take all components from the set above between [0-500Hz] in which case you end up with [100 400].

Do you see the 400Hz odd man out here?This is termed as the alias of the 600Hz signal.

So, for an input of [100 600] and you end up with [100 400].

  b) for $f_0=1.3\,\text{kHz}$ reconstruction works

Repeat steps above.

You'll find that the only components between [0-650Hz] are [100 600].

There you go !


Does resampling cause phase shifts? I have a signal sampled at 40MHz, I would like to resample it to 37MHz. The signal is not periodic, I did resampling with Matalb resample function and it doesn't cause phase shift (as far as I understood). Matlab applies an anti-aliasing FIR filter and compensate for the delay introduced by the process.

I want to do the same process in Python. I know there is a resample function in scipy.signal but the documentation is not clear for me. Does the resample function introduces phase shifts?

if so, there is also a decimate function in scipy as well which applies similar FIR filter. Should I use decimate function over resample?

Upsample first by 12 and then downsample by 13, you'll get very close to 37Mhz, 36.92MHz to be precise, should be good enough.

If you want exactly 37Mhz, then upsample by 37 and then downsample by 40, but this is more calculation intensive with not much gain, so better opt for the first strategy. Also antialsiaing demands are less for the first case.


What is the Fourier transform of $\delta(t-a)$ As I know the Fourier transform of $\delta(t)$ is equal to $1$, on the other hand, shift with $a$ in time domain is equal to multiplying with $e^{-jaw}$ so the Fourier transform of $\delta(t-a)$ is equal to $e^{-jaw}$, but on the other hand the Fourier transform of Dirac comb ($\sum_{n=-\infty}^{+\infty}\delta(t-nT)$which is equal to sum of Dirac functions) is also another Dirac comb, This doesn't make sense to me. The Fourier transform of Dirac comb should be sum of the $e^{-jnTw}$.

Your argumentation is perfectly right. The only thing you are missing is that the sum of complex exponentials becomes a comb function in the limit. See also this link and this matlab code:

N = 2048;
n = (0:(N-1))';

k = 4*(0:20);

signals = zeros(N,length(k));
for kk = 1:length(k)
    signals(:,kk) = exp(-1j*2*pi*k(kk).*n/N);
end

plot(abs(sum(signals, 2)));


Removing periodic noise from image using Fourier transform I'm trying to get rid of some periodic flicker noise through post-processing of the recorded images. The reason for these artifacts is that the electronic rolling shutter of the camera reads each line sequentially, and this mechanism results in brightness variations in the image, when the picture is taken under fluorescent lighting (due to 50 or 60 Hz AC power frequency).

I tried calculating the Fourier transform of the image and suppressing the components that might cause flickering. I inspected the magnitude spectrum (as you can see below for an example image) and set the vertical components around the DC component to a very small value (I preserved the DC component). In the end, I managed to eliminate the flicker because I suppressed the flicker frequency, but also the image was sort of deformed since I, probably, also suppressed some other frequencies that contributed to the signal's energy.

I would like to know what would be a better way to detect the flicker frequency and eliminate/suppress only that frequency? In my current implementation, I could be doing something fundamentally wrong, so I would appreciate any guidance and hints. Below, I show the resulting image that I got for the kind of filtering that I explained above. You can also find my Python code below.

import cv2
import numpy as np
from matplotlib import pyplot as plt

img = cv2.imread('images/flicker2.jpg',0)
f = np.fft.fft2(img)
fshift = np.fft.fftshift(f)
# calculate amplitude spectrum
mag_spec = 20*np.log(np.abs(fshift))

r = f.shape[0]/2        # number of rows/2
c = f.shape[1]/2        # number of columns/2
p = 3
n = 1                   # to suppress all except for the DC component
fshift2 = np.copy(fshift)

# suppress upper part
fshift2[0:r-n , c-p:c+p] = 0.001
# suppress lower part
fshift2[r+n:r+r, c-p:c+p] = 0.001
# calculate new amplitude spectrum
mag_spec2 = 20*np.log(np.abs(fshift2))
inv_fshift = np.fft.ifftshift(fshift2)
# reconstruct image
img_recon = np.real(np.fft.ifft2(inv_fshift))

plt.subplot(131),plt.imshow(img, cmap = 'gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(132),plt.imshow(mag_spec, cmap = 'gray')
plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([])
plt.subplot(133),plt.imshow(mag_spec2, cmap = 'gray')
plt.title('Magnitude Spectrum after suppression'), plt.xticks([]), plt.yticks([])
plt.show()

plt.figure()
plt.subplot(121),plt.imshow(img, cmap = 'gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(122),plt.imshow(img_recon, cmap = 'gray')
plt.title('Output Image'), plt.xticks([]), plt.yticks([])
plt.show()

Just eyeballing the first image tells me that you have a very simple band pattern. They're in fact mostly horizontal (in the spatial domain). You see a clear cross in the Fourier domain. That's because most of your spatial features in that hallway are oriented horizontally and vertically. Your horizontal band similarly show up as a pair of dots above and below the origin. The bands are slightly rotated to the right, and so will the pair of dots.

You'll need to erase a very small area around these dots. But how to locate them? The easiest solution may be to take your recorded image, open it in Paint or a similar program, and paint pure black and white bars over the actual bands. Now take the FFT of this bar pattern, and you get the exact location of the two points you're looking for. It in general won't be exactly a single FFT bin, so you could fit a gaussian curve to the neighborhood of those two points, and use the opposite 1-gauss(x-xc,y-yc) as the frequency domain filter.


Convolution of discrete signals using convolution sum I want to perform the convolution of the following discrete signals: $$h[n]=u[n-2] $$ and $$x[n] = (0.5)^nu[n+2]$$.
That's what I've done so far: $$\sum_{k=-\infty}^{\infty} (0.5)^ku[k+2]*u[n-2-k]$$
So: $$\sum_{k=0}^{\infty} (0.5)^ku[k+2]*u[n-2-k]$$
I know that the upper limit does not go to infinity, however, I cannot determine the correct upper limit. Any thoughts?

Well, you know that $u[n]=0$ for all $n<0$. So from the two terms in your sum, all terms will be zero where either $k+2 < 0$ (i.e., $k < -2$) or $n-2-k<0$ (i.e., $k>n-2$). From this your sum will need to run from $-2$ to $n-2$. Can you take it from here?


EEG signal processing I have an EEG data set downloaded from physionet. I want to play that data into a brain monitoring device like for example BIS monitor or Narcotrend or any other similar brain monitoring device.I use Matlab GUI to visualize the data. Can I send it from a Matlab GUI to any one of the brain monitoring devices with a suitable hardware  in between.Any ideas please I am new to EEG signals.

I'll assume that you want to do online monitoring.

You can always store the coming data from the GUI in a vector in MATLABA.
Then you can run simple code to plot them.

This also depends on the type of measurements you have. If the set is capable on detecting power ranges as well you should have the following readings:

High Alpha, Low Alpha, High Beta, Low Beta, High Gamma, Low Gamma, Delta, Theta, Attention, and Meditation.


Is a Kalman filter suitable to filter projected points positions, given Euler angles of the capturing device? My system is the following. I use the camera of a mobile device to track an object. From this tracking, I get four 3D points that I project on the screen, to get four 2D points. These 8 values are kinda noisy, due to the detection, so I want to filter them to make the movement smoother and more realistic. As a second measurement, I use the device's gyroscope output, which provides three Euler angles (i.e. the device attitude). Those are more precise and at greater frequency (up to 100 Hz) than 2D positions (around 20 Hz).

My first attempt was with a simple low-pass filter, but the lag was important, so I now try to use a Kalman filter, hoping it will be able to smooth the positions with little delay. As seen in a previous question, one key point in a Kalman filter is the relation between the measurements and the internal state variables. Here the measurements are both my 8 2D point coordinates and the 3 Euler angles, but I'm not sure about what I should use as internal state variables and how I should connect the Euler angles to the 2D points. Hence the primary question, is a Kalman filter even suitable for this problem? And if yes, how?

Low Pass Filtering

It would be good to know what you mean by "simple low pass filter".

For example, if your measurements at time $k$ are

$$
p_{k} = \left[ \begin{array}{c} x_k \\ y_k \end{array} \right]
$$

and your low pass filtered estimates are:

$$
p^{\tt LPF}_k = \alpha p^{\tt LPF}_{k-1} + (1-\alpha)p_k
$$

then you will have quite a large group delay in the filter of about $1/(1-\alpha)$ (for alpha close to 1).

Modeling the Signal: Simplistic Approach

To use the Kalman filter (or any similar approach), you need to have a model for how your measurements are acquired and updated.

Usually this looks like:

$$
p^{\tt TRUE}_{k+1} = {\mathbf A} p^{\tt TRUE}_k + {\mathbf B} \epsilon_k
$$
where $\epsilon_k$ is the process (driving) noise, ${\mathbf A}$ is the state transition matrix, and ${\mathbf B}$ is your input matrix.

And then your measured $p_k$ are:
$$
p_k = {\mathbf C} p^{\tt TRUE}_k + {\mathbf D} \nu_k
$$
where $\nu_k$ is the output (measurement) noise,  ${\mathbf C}$ is the output matrix, and ${\mathbf D}$ is your measurement noise matrix.

Here, the "state" of the model is chosen as the true positions, and the things you measure are the output.

You can then apply the Kalman filter equations to this to get state estimates $\hat{p^{\tt TRUE}_k }$ of the true position.

However, this approach is simplistic because it does not use any knowledge of how the points might move (nor does it use your 4 points and any knowledge you might have about how they move together).

Modeling the Signal: Starting a Better Approach

This page shows how to set up the problem involving the positions and euler angles.
It's doing something different from what you need, but the state is:

$$
p^{\tt TRUE}_{k} = \left[
x_k\  y_k\  z_k\
\dot{x}_k\  \dot{y}_k\  \dot{z}_k\
\ddot{x}_k\  \ddot{y}_k\  \ddot{z}_k\
\phi\ \psi\ \theta\
\dot{\phi}\ \dot{\psi}\ \dot{\theta}\
\ddot{\phi}\ \ddot{\psi}\ \ddot{\theta}\
 \right]^T
$$

and the measurements (output) is

$$
p_k = \left[ x_k\  y_k\  z_k\  \phi\ \psi\ \theta\ \right ]^T
$$

All the model on that page is really doing is saying:
$$
x^{\tt TRUE}_k = \sum_{n=0}^k \dot{x}^{\tt TRUE}_n n\Delta t + \frac{1}{2}  \sum_{n=0}^k \ddot{x}^{\tt TRUE}_n (n\Delta t)^2
$$
(but for each of $x,y,$ and $z$).

This is just the classic "equations of motion". See equation (3) here.


Efficient FFT computation of a zero-padded vector I don't think this question has a good answer but will ask nevertheless since it has been bothering me for a few days.

I am interested in computing as efficiently as possible the N-point FFT of a N/2 vector that has been zero-padded with N/2 zeros.  (I am actually interested in doing the same for vectors zero-padded with many more zeros but let's discuss the simplest case simple here.)

A simple obvious optimization consists in skipping half of the computations in the first FFT stage which correspond to multiplying a zero input value by a twiddle factor.  One could further skip a quarter of the computations that have zeros in the second stage, and so on, but the savings quickly dwindle down.

I was hoping we could do better.  After all, an FFT of a zero-padded vector is essentially a sinc interpolation of a smaller FFT taken on the non-zero smaller input vector.

In fact, writing down the math shows exactly that.  If I call the FFT of the non-zero N/2 vector y, the FFT of the zero-padded vector ypad, @ represents the circular convolution, and SINC the N/2 FFT of the first N/2 twiddles of an N-point FFT, we have:

ypad[2k]   = y[k]
ypad[2k+1] = ( y  @ SINC )[k]

Problem is: computing the circular convolution of y by the SINC function for each odd point is very costly and doesn't lead to a faster implementation that computing the N-point FFT.  My last hope is that once we have computed the circular convolution in one point, the circular convolution for other points can be derived with fewer computations, but I currently don't think so.

Has anyone an idea on how to compute efficiently the N-point FFT from the N/2-point FFT of the non-zero vector or is that a dead end?

You can decompose the FFT that you describe into some smaller transforms. Take a look at what you're calculating when you zero-pad a vector $x[n]$ of length $N \over 2$ to length $N$ and calculate a DFT:

$$
X[k] = \sum_{n=0}^{\frac{N}{2}-1} x[n] e^{-j2\pi nk \over N}
$$

To split this up into two transforms, first look at just the even-indexed terms in the $N$-point sequence $X[k]$:

$$
\begin{align}
X[2k] &= \sum_{n=0}^{\frac{N}{2}-1} x[n] e^{-j2\pi n2k \over N} \\
&= \sum_{n=0}^{\frac{N}{2}-1} x[n] e^{-j2\pi nk \over \frac{N}{2}}
\end{align}
$$

Via inspection of the above, we can conclude that the even-indexed terms of the zero-padded DFT are just equal to the $\frac{N}{2}$-point DFT of the original input sequence (with no zero-padding). What about the odd-indexed terms?

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

This is a similar expression to the even-indexed terms $X[2k]$. We take the original $N\over 2$-length sequence, multiply it by a complex exponential sequence $e^{-j2\pi n \over N}$, then calculate the $N \over 2$-length DFT of the result. This makes intuitive sense if we write the exponential as $e^{-j2\pi 0.5n \over \frac{N}{2}}$: multiplying by such a complex exponential yields a shift in the frequency domain of half of a $N \over 2$-length DFT bin width. This gives the interpolation-by-2 effect that the zero-padding is known to yield.

In summary: To decompose an $N$-point DFT that consists of an $N \over 2$-point signal $x[n]$ followed by $N \over 2$ zeros, do the following:

Calculate an $N \over 2$-point DFT of the original signal $x[n]$. This make up the even-indexed values in the zero-padded DFT result.
Multiply the original signal by the complex exponential function $e^{-j2\pi n \over N}$. Calculate an $N \over 2$-point DFT of the product. This makes up the odd-indexed values in the zero-padded DFT result.

The caveat: While it depends on your DFT size $N$, this may not effect any tangible reduction in complexity beyond just calculating the zero-padded DFT directly. I have an application where I've tried to make a similar optimization myself in the past, although it is on a very different platform from what you described (mine runs on high-performance server machines using highly-optimized FFT libraries like Intel Math Kernel Library). While you may be able to achieve a theoretical reduction in the total number of arithmetic operations that you need, due to other effects like memory accesses, caching, and so on, it's common that trying to do 2 $N \over 2$-point FFTs is slower than just doing one $N$-point FFT directly. In my case, I've never been able to beat the existing library's $N$-point transform.

To estimate the complexity reduction by doing this, we can estimate the number of arithmetic operations that each method needs. Using the rule of thumb that an $N$-point FFT requires approximately $5 N \log N$ arithmetic operations, we can estimate the savings in operations:

$$
\text{savings} = 5 N \log N - \left(2\left(5 \frac{N}{2} \log \frac{N}{2} \right)+5 \frac{N}{2}\right)
$$

$$
\text{savings} = 5 N \log N - \left(5 N \log \frac{N}{2} +5 \frac{N}{2}\right)
$$

The first term in the parentheses counts the approximate number of operations required by the 2 $\frac{N}{2}$-point DFTs, while the second term counts the number of operations needed for the $\frac{N}{2}$-point complex vector multiplication (there are 5 scalar operations per complex multiply).

The problem is, as $N$ increases, the amount of savings diminishes because of the logarithmic scaling of the DFT. For typical DFT sizes, you will likely find that this approach yields little complexity reduction.


Uncorrelated model error A dynamical system evolves according to $x_{k+1} = M_k x_k + w_k$

$w_k$ denotes the model error.

In a textbook, it is specified that $w_k$ is temporally uncorrelated and

$E[w_kw_j^T]$ = $Q_k$ if j=k and zero otherwise; $Q_k$ is n x n symmetric positive definite

Could someone explain why $w_k$ is temporally uncorrelated and how $Q_k$ is symmetric positive definite?

If $\mathbf{w_k} $ was temporally correlated, it would have memory and would be equivalent to an augmented state.  In other words, it can be written as a state variable with its own independent noise term. You would be able to rewrite the equation to the same form you originally showed.  It is the most basic form of the model.

The symmetry of $Q_k$ is because $E\{ ab\}=E\{ba\}$ for any terms $a$ and $b$ of the vector $\mathbf{w_k}$

Edit to address comments.

Let $\mathbf{x}=\mid x_0,x_1, \dots x_{N-1} \mid^T$ be a random Gaussian distributed random vector.
$$
E\left\{ \mathbf{x} \mathbf{x}^T \right\}=\mathbf{R} \quad \text{Is typically Full rank and also typically} = \sigma^2 \mathbf{I}
$$
when the elements of $\mathbf{x}$ are independent.
There are some exceptions to full rank, like like a linear dependence of some of the elements of $\mathbf{x }$ such as $x_i=x_j$ of more generally $x_i=\sum_{k=0, k\ne i}^{N-1} \alpha_k x_k$ .

The easiest example is to consider $E\{ x_i x_i \} = E\{ x_i^2\}$  .  A Gaussian random variable can have zero mean but a variance $\sigma^2$. The matrix $\mathbf{R}$ would have nonzero diagonal terms, and is Full Rank.

The Matrix $\mathbf{R}$ isn't required to be diagonal and has nonzero off diagonal terms when $E\{x_i x_j\} \ne 0$ for some $i$ and $j$.  For your problem, your book will likely show the covariance of the evolution of the the state vector and that will most likely have off diagonal nonzero elements.
$$
E\left\{ \mathbf{x}\right\} E\left\{ \mathbf{x}^T \right\} \quad \text{Is rank one}
$$


Determining the fundamental frequency/pitch of a note What I am trying to achieve is getting the fundamental frequency of a note played by an instrument. What I have already done is performing an FFT on a samples of audio file, and here's what I get:

After that I just tried to find the maximum value in that dataset and I get "4138". Unfortunately, the note being played in the audio file is an E and from internet I know that the E note should have 5274.04 Hz or 2637.02 Hz.
I continued reading some additional info on the internet and found interesting info that in some cases fundamental frequencies could be missing from the sound, so I have tried doing some autocorrelation on the samples I got from the audio file to find the fundamental frequency and I got something like this:

It looked like a frequency domain for me so I performed FFT on that dataset but got similar results.
Also I have read that some people suggest using HPS (Harmonic Product Spectrum) for finding the fundamental frequency, I have tried that but obviously the result was the same, as the FFT itself didn't change.
P.S. Of course I use Hamming windowing function for samples. The library I use for FFT is JTransforms.

Ok, everything was working quite fine, the only thing that I have missed is another formula that is used to get the "real" frequency from the FFT power spectrum. And the formula is as follows:

  frequency = indexMax * Fs / L;
  where:
  indexMax - index of the maximum in power spectrum;
  Fs - sampling rate of an audio file;
  L - length of the power spectrum array.

UPDATE
This approach is working fine for such trivial tasks as making a tuner or something like that, but as

  robert bristow-johnson

pointed out, this implementation will, in some cases, ignore the fundamental frequency and potentially show you frequency of one of the harmonics. So a better approach would be

Autocorrelation
AMDF
ASDF
YIN
Mcleod Pitch method


$l_2/l_2$ guarantee on sparse Fourier transform I am starting my studies now on signal processing, and really didn't find nothing on "$l_2/l_2$ guarantee" of a certain function, in my case:

$$||\hat{x} - \hat{x}'||_2 \leq C\text{ min }_{\text{k-sparse}\; y} ||\hat{x}-y||_2 $$

My reference is "Simple and Practical Algorithm for Sparse Fourier Transform" by Hassanieh, Indyk, Katabi, Price.

What does it means to say that the function is $l_2/l_2$ guarantee or $l_{\infty}/l_2$ guarantee? A reference would be of great help!

A signal $s$ approximation from $x$ is often quantified by   the closeness of the approximant $s'$ and the desired behavior. Traditionally, one uses norms ($\ell_p$, $p\ge 1$),  pseudo-norms ($0<p < 1$) or the so-called $\ell_0$ count index. A very traditional formula uses a linear transformation $\phi$, and the approximation operator $A$, $s'=A(\phi x))$:

$$ \ell_q(s -A(\phi x))\,.$$

Minimization is often looked for in specific signal classes.
Looking for a $k$-sparse approximation is notoriously complicated: satisfying that $s'$ has at most $k$ nonzero components is generally NP-hard (B. K. Natarajan, Sparse Approximate Solutions to Linear Systems, 1995). So one may want a faster  approximant to be bounded by a measure based on the error with the best $k$-sparse  solution:

$$ \ell_q(s -s')\ \le C\min_{s_k} \ell_p(s -s_k)\,,$$

where $C$ is a constant as close to $1$ as possible. Sometimes, this error is guaranteed with some probability $p$. One often finds the terms "mixed norm error guarantee".

Many references can be found in the compressed sensing literature. On of the early source is A. Cohen, W. Dahmen, and R. DeVore. Compressed Sensing and best k-term approximation, 2006. Other references include $\ell_2/\ell_2$-for each sparse recovery with low risk, A. C. Gilbert et al.  You can find an $\ell_2/\ell_1$ in
  One sketch for all: Fast algorithms for Compressed Sensing, A. C. Gilbert, M. J. Strauss, J. A. Tropp.

A nice survey could be Sparse Recovery Using Sparse Matrices, Anna Gilbert, Piotr Indyk.


limiter algorithm wanted: mixing two channels without clipping I have two audio tracks (instrumental recordings) named track A, and track B. The tracks will be played in an application and the user of this application has the ability to change the balance of these tracks. The user controls a slider ranging from $x=0$ to $x=1$ with which the balance can be changed from hearing track A only, to hearing track A and B mixed equally, to hearing track B only. Mathematically, I have volume functions for track A and B, $V_A(x)$ and $V_B(x)$ respectively, that should satisfy:

$V_A(0) = 1$
$V_A(1) = 0$
$V_A(0.5)=V_B(0.5)$
$V_B(0) = 0$
$V_B(1) = 1.$

The resulting track R is calculated as follows:

$R(t) = A(t) V_A(x) + B(t) V_B(x)$ where $0 \le x \le 1.$

A naive approach would be the to let $V_A(x) = 1 - x$ and $V_B(x) = x$ but I don't want the volume of track A to go down drastically as the user moves to the equal mix at $x = 0.5$. A better approach is $V_A(x) = \sqrt{1 - x}$ and $V_B(x) = \sqrt{x}$ (see How to mix two signals without changing the overall loudness?). But ideally, I would have the following function (with the non-constant parts possibly some other power function):

$V_A(x) = \left\{
 \begin{array}{ll}
  1  & \mbox{if } x \le 0.5 \\
  2  (1 - x) & \mbox{if } x > 0
 \end{array}
\right.$

$V_B(x) = \left\{
 \begin{array}{ll}
  2x  & \mbox{if } x \le 0.5 \\
     1 & \mbox{if } x > 0.5
 \end{array}
\right.$

The problem is that both track A and track B are mastered at -0.1 dBFS (they will be played on a device and need to be as loud as possible), so obviously the last two functions will (most likely in my case) cause clipping.

I have two ideas on how to handle this:

Use my ideal function and apply some type of limiter (or would it be better to limit track A and B first, and then mix?)
Find the highest value for $V_A(0.5) = V_B(0.5)$ such that $R(t)$ never clips and then fill in the volume curve somehow. I can even limit some peaks so that this value could probably be quite high.

What are people's thoughts on this? My specific question is about idea 1. How does one go about creating a limiter? Does anyone have a good reference to some algorithms?

You can keep the mixing itself linear but add a mapping to the mixing variable $x$ so that you can control its dynamics. The closest this is to reality is the use of logarithmic potentiometers.

So, given $C_1, C_2$ being your mixing channels and $x$ being your mixing variable, the mixed track $M$ at some time instance $n$ is:

$$ M[n] = C_1[n] \cdot x + C_2[n] \cdot (1-x)$$

At $x=0.5$, both tracks are contributing equally, provided that they are mixed to similar levels. If the tracks are not mixed to similar levels, then this inserts a bias on $x$ on either side.

Anyway, so now another "problem" arises which is how to control $x$ and even how to limit it.

To achieve this, we setup a mapping via a function of the form $y = f(x)$. Ideally, $x,y \in \left[0 \ldots 1 \right]$ So now our mixing process looks like this:

$$ M[n] = C_1[n] \cdot f(x) + C_2[n] \cdot (1-f(x))$$

(The right part could also have been $f(1-x)$)

If you want $f$ to limit your values in a subset of $x$, then you could use a sigmoid function. The sigmoid function has a linear part (where the mapping between $y,x$ is almost identical) and two "non-linear" parts called the shoulder (at the top) and the knee (at the bottom). These are the points where the function starts to plateau and no matter what the value of $x$ is, the value of $y$ is pinned down to a specific limit.

Now, obviously, for some sigmoid function $f$, its output $y \in \left[-1 \ldots 1\right]$. But that is not much of a problem because you can simply $y = \frac{f(x)+1}{2}$ and now your $y$ goes from $0$ to $1$. If you don't want it to hit $1$, then you can now multiply the output of your $f$ with another scaling factor. So $y = q \cdot \frac{f(x)+1}{2}$, where $q$ could be $0.8$ in your case.

Another thing you can do with this scheme is control the "dynamics" of the mixing. So, the user might be moving the slider linearly but their effect taking place non-linearly or with some other rate. This could account for the fact that perception is not linear (for example).

For this purpose, you can set up a logarithmic function $f$ that goes from 0 to 1. The logarithmic also has a "limiting" region towards its maximum value.

The general idea is to keep the mixing linear but not necessarily the way that the mixing variable varies.

The other thing you might want to consider if your tracks are "difficult" might be the use of a compressor before the mixing, so that the levels of both audio channels are kept within reasonable limits before they get mixed down. But, inserting a compressor introduces more controls (more variables) that someone would have to fiddle with depending on the tracks. Maybe you don't want this (?)

Hope this helps.


How we can represent a bandpass signal equation in matlab The general form of bandpass signal is given by
$$s(t) = \Re\{S_b(t)\exp(j2 \pi f_ct)\}$$
How this can be represented in matlab.

Usually (almost always), simulation is done using the discrete time baseband signal. Lets say you want to simulate 32 BPSK symbols (either +1 or -1). The MATLAB code for this is:

 symbols = sign( rand(32, 1) - 0.5 );

Now that you have the symbols, the next step is to upsample. You can either use the MATLAB function or do it yourself:

 symbols_up = upsample(symbols, samplesPerSymbol);

Then create and use a pulse shaping filter (just use rcosdesign):

 pulseFilter = rcosdesign(rolloff, span, samplesPerSymbol, 'sqrt');
  txSig = filter(pulseFilter, 1, symbols_up);
  figure, plot(txSig)


Programmatically (e.g. in Python) adding clipping to an audio file? I am trying to replicate the audio artifact exemplified here to a set of 1000 audio files. In particular, I would like to be able to 'clip' the louder parts of each sound clip so they become distorted as in this video. Is there a way I can do this with a script?

Thanks!

Simple clipping :
threshold = 0.5
If x > threshold
  x = threshold
elseif x < -threshold
  x = -threshold
end

Real-world clipping can be significantly more complex than this, involving various time-constants, asymmetry, heating-effects,....


low pass filter and FFT for beginners with Python I am new to signal processing and especially to FFT, hence I am not sure if I am doing the correct thing here and I am a bit confused with the result.

I have a discrete real function (measurement data) and want to set up a low pass filter on that. The tool of choice is Python with the numpy package. I follow this procedure:

compute the fft of my function
cut off high frequencies
perform the inverse fft

Here is the code that I am using:

import numpy as np
sampling_length = 15.0*60.0 # measured every 15 minutes
Fs = 1.0/sampling_length
ls = range(len(data)) # data contains the function
freq = np.fft.fftfreq(len(data), d = sampling_length)
fft = np.fft.fft(data)
x = freq[:len(data)/2]
for i in range(len(x)):
if x[i] > 0.005: # cut off all frequencies higher than 0.005
    fft[i] = 0.0
    fft[len(data)/2 + i] = 0.0
inverse = np.fft.ifft(fft)

Is this the correct procedure? The result inverse contains complex values, that confuses me.

The fact that the result is complex is to be expected. I want to point out a couple things:

You are applying a brick-wall frequency-domain filter to the data, attempting to zero out all FFT outputs that correspond to a frequency greater than 0.005 Hz, then inverse-transforming to get a time-domain signal again. In order for the result to be real, then the input to the inverse FFT must be conjugate symmetric. This means that for a length-$N$ FFT,

$$
X[k] = X^*[N-k], k = 1, 2, \ldots , \frac{N}{2} - 1\;\;\;\;\;\;\;(N\;\; even)
$$

$$
X[k] = X^*[N-k], k = 1, 2, \ldots , \lfloor \frac{N}{2} \rfloor\;\;\;\;\;\;\;\;(N\;\; odd)
$$

Note that for $N$ even, $X[0]$ and $X[\frac{N}{2}]$ are not equal in general, but they are both real. For odd $N$, $X[0]$ must be real.

I see that you attempted to do something like this in your code above, but it is not quite correct. If you enforce the above condition on the signal that you pass to the inverse FFT, then you should get a real signal out.

My second point is more of a philosophical one: what you're doing will work, in that it will suppress the frequency-domain content that you don't want. However, this is not typically the way a lowpass filter would be implemented in practice. As I mentioned before, what you're doing is essentially applying a filter that has a brick-wall (i.e. perfectly rectangular) magnitude response. The impulse response of such a filter has a $sinc(x)$ shape. Since multiplication in the frequency domain is equivalent to (in the case of using the DFT, circular) convolution in the time domain, this operation is equivalent to convolving the time domain signal with a $sinc$ function.

Why is this a problem? Recall what the $sinc$ function looks like in the time domain (below image shamelessly borrowed from Wikipedia):

The $sinc$ function has very broad support in the time domain; it decays very slowly as you move in time away from its main lobe. For many applications, this is not a desirable property; when you convolve a signal with a $sinc$, the effects of the slowly-decaying sidelobes will often be apparent in the time-domain form of the filtered output signal. This sort of effect is often referred to as ringing. If you know what you're doing, there are some instances where this type of filtering might be appropriate, but in the general case, it's not what you want.

There are more practical means of applying lowpass filters, both in the time and frequency domains. Finite impulse response and infinite impulse response filters can be applied directly using their difference equation representation. Or, if your filter has a sufficiently-long impulse response, you can often obtain performance benefits using fast convolution techniques based on the FFT (applying the filter by multiplying in the frequency domain instead of convolution in the time domain), like the overlap-save and overlap-add methods.


FFT question as it relates to Neural Networks / Supervised Learning Models So I've been researching wavelet transforms and FFT.  I want to feed wavelet transforms of a 1D signal in time into a neural net and train against a target variable at each time step.  The idea being the WT/FFT can help the net pick up on data structures.

Only issue being I have no idea how to feed the coefficients of the FFT into the model because my approximation signal is length 87 while the original series had 44,000 samples.

I looked into zero-padding to restore the number of samples but that is apparently a technique used on the time-series itself before decomposing with FFT.

Any help appreciated, hard to find material on exactly this subject.

paper im recreating: https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0180944

I stumbled across this.
http://ataspinar.com/2018/12/21/a-guide-for-using-the-wavelet-transform-in-machine-learning/

hope it helps
john


Energy per symbol calculation I am having trouble calculating my Energy per symbol and Es/No. Using these parameters and equation
$N_s = $Number Symbols
$R_s = $Symbol Rate
$T_s = $Symbol Time
$sps = $samples per symbol
x = signal = QAM-16 with no noise added
Then using these, two equations:
$E_s = E_{total} / N_s = sum(abs(x)^2) / N_s $
$E_s=P_{avg} / (sps/R_s)  = P_{avg} sps T_s$
This second equation comes from Energy = Power x Time. So Energy per symbol = Power x Time for Symbol, but my symbol spans more samples because its upsampled, one symbol spans sps samples. So Energy per symbol = Power x time for symbol x sps.
Where
$P_{avg} = sum(abs(x)^2) / N $
N = Number Samples
I have tried both equations after the RRC transmit and after the RRC receive filters, both times these equations give two different answers.
First equation gives a very large number and the second is incredibly small number in comparison. Ive tried removing sps from the second equation but does not get anywhere close to first equation.
Ultimately I need to work out Es/No and I cant get Es right.  Any thoughts please? :/

The first equation, $ \frac{1}{N_s}\sum |x[n]|^2 $, is correct. Applying the RRC filter kind of smears the energy of the symbol over several samples so that when you add them all back up again you get $E_s$, and of course divide by $N_s$ for multiple symbols. The RRC filter should be normalized as someone in the comments noted.
The second equation I can't say much as it is not clear how you compute $P_{avg}$. Secondly, dimensional analysis is a fancy way of saying writing out the units and seeing what cancels out. For example, in the second equation:
$$\text{sps }T_s\rightarrow \frac{\text{samples}}{\text{symbol}}\frac{\text{seconds}}{\text{symbol}}=\frac{\text{(samples)(seconds)}}{(\text{symbol})^2}$$
I've found this method useful as a "sanity check" for myself many times. The above units should suggest to you that something is amiss (also why do you need two equations for the same quantity?).
Edit
The question was updated to include definition of $P_{avg}$. The OP is computing average power per sample and trying to convert it to average power per symbol. Again, write out the units it becomes clear:
$$\frac{\text{energy}}{\text{sample}}\frac{\text{samples}}{\text{symbol}}=\frac{\text{energy}}{\text{symbol}} $$
sps = 4;
Ns = 100;
Es = 10;
rrc = rcosdesign(0.25, 6, sps, 'sqrt');
s = sqrt(Es)*qammod(randi([0, 15], Ns, 1), 16, 'UnitAveragePower', true);
x = upfirdn(s, rrc, sps);

equation1 = sum(abs(x).^2)/Ns;
equation2 = sum(abs(x).^2)/length(x)*sps;


what is the relationship between a spectrogram and the uncertainty principle heisenberg? what is the relationship between these two things
Perhaps more resolution in a spectrogram is equivalent to knowing more the position of the electron and less resolution is knowing the velocity of the electron.

In Quantum Mechanics, state of system specified by a vector (wave function which is a vector in function space), and you could use different basis to represent this vector (imagine one vector in 2 different coordinate system which lead to different components in those coordinate system but they both represent the same vector). For systems composed from moving particles, one basis is the position basis and the second one is momentum basis. The transform between these two basis is the Fourier transform. Considering the wave function give us a probabilistic interpretation, to find the momentum or position of a particle you have to find the average and after that to have a sense about error, you have to find the standard deviation.
The Heisenberg Uncertainty principle, tell us the multiplication of position error and momentum error could not be smaller than some value, for any possible state of system.
Considering this is a property of Fourier transform you could extend that to signals where the signal amplitude and it's spectrum are different representation of same thing in different basis, and say multiplication of effective bandwidth of pulse (as standard deviation of power spectrum around central frequency) and the effective pulse width (as standard deviation of signal's power around it's center in time) could not be smaller than some value for every possible signal.


How can we prove that Spectral density of In-Phase component and Quad-Phase component of a sample function of narrow-band noise have same definition? Suppose we have a narrow-band noise $N(t)$. We analyse this to In-Phase and Quad phase components by the equation:

$$N(t) = N_I(t)  \cos(2\pi f_c t) - N_Q(t) \sin(2\pi f_c t)$$

Now a book (by Simon Haykin) states that:

$$ S_{N_I}(f) = S_{N_Q}(f) =
\begin{cases}
 S_N(f-f_c) + S_N(f+f_c) & -B\le f \le B \\
 0 & \text{otherwise}
\end{cases}
$$

How do we prove this?

This is the Image of Power Spectral Density of Narrowband Noise.

Synthetic Method:

If $\{\hat X(t)\}$ and $\{\hat Y(t)\}$ are zero-mean uncorrelated  low-pass WSS processes with identical autocorrelation function $R(\tau)$ and identical power spectral density $S(f)$ enjoying the property that $S(f) = 0$ for $|f|>B$, then
$$\hat{N}(t) = \hat X(t)\cos(2\pi f_ct) - \hat Y(t)\sin(2\pi f)ct$$ is a band-pass process whose autocorrelation function is $R(\tau)\cos(2\pi f_c\tau)$ and power spectral density is
$${S}_{\hat{N}}(f) = S(f-f_c)+S(f+f_c) = S_N(f).$$
There is no pretense that $\hat{N}(t)$ is the same as $N(t)$ for any $t$ or that $\hat X(t)$ and $\hat Y(t)$ are the complex baseband representation of $N(t)$ but the point is that by looking only at the power spectral densities, one cannot tell whether the synthesized process $\{\hat N(t)\}$ is the same as, or different from, the given bandpass process $\{N(t)\}$.

Analytic Method: (Suspend disbelief and take with a large grain of NaCl)

Suppose that $\{N(t)\}$ is a zero-mean bandpass WSS process whose power spectral density $S_N(f)$ has value $0$ except when $|f| \in [f_c-B, f_c+B]$ and consider the random variables
$N(t)\cos(2\pi f_ct)$ and $-N(t)\sin(2\pi f_ct)$. They obviously have zero mean, and if $N(t)$ has variance $\sigma^2 = R_N(0)$, then these variables have variances \begin{align}\sigma^2\cos^2(2\pi f_ct) &= \frac 12\big(\sigma^2+\cos(2\pi (2f_c)t)\big)\tag{1}\\
\sigma^2\sin^2(2\pi f_ct) &= \frac 12\big(\sigma^2-\cos(2\pi (2f_c)t)\big)\tag{2}.
\end{align}
Their covariance is $$-\sigma^2\cos(2\pi f_c t)\sin(2\pi f_c t)
= -\frac 12\sigma^2\sin(2\pi (2f_c)t).\tag{3}$$

We take the viewpoint that each of the random processes $\{N(t)\cos(2\pi f_ct)\}$ and $\{-N(t)\sin(2\pi f_ct)\}$, neither of which is a WSS process, is the sum of a low-pass WSS process $\{X(t)\}$ (respectively $\{Y(t)\}$) and a "double-frequency" process which not WSS, and which we can eliminate by (ideal) low-pass filtering these two processes with bandwidth $B$. Put another way, the result of (separately) low-pass filtering the processes $\{N(t)\cos(2\pi f_ct)\}$ and $\{-N(t)\sin(2\pi f_ct)\}$ leaves us with WSS processes with variance $\frac 12\sigma^2$ each (the double-frequency stuff in $(1)$ and $(2)$ has been eliminated by the low-pass filtering), and these processes are uncorrelated (the nonzero covariance in $(3)$ is double-frequency stuff that has been eliminated by the low-pass filtering).

Multiplying by $\cos(2\pi f_ct)$ or $\sin(2\pi f_ct)$ is modulation (frequency content at frequency $f$ is moved to frequencies $f+f_c$ and $f-f_c$).  Applying this notion to the frequency band $[f_c-B, f_c+B]$, we see that any content there is moved down to be centered at $0$ and also moved to the band $[2f_c-B, 2f_c+B]$. Similarly, any content in the band $[-f_c-B, -f_c+B]$ is moved down to be centered at $0$ and also moved to the frequency band $[-f_c-B, -f_c+B]$. The latter is the stuff that gets eliminated by the low-pass filtering.  What is left is the stuff moved down to DC and this "explains" why Haykin claims that $$ S_{N_I}(f) = S_{N_Q}(f) =
\begin{cases}
 S_N(f-f_c) + S_N(f+f_c) & -B\le f \le B \\
 0 & \text{otherwise}
\end{cases}
$$


Highpass filter square wave Generated a 1kHz square wave in Audacity and normalized it to 0dB (fs = 44100Hz):

Added a 24dB/oct high pass filter (std Audacity plug-in), cutoff frequency set to 100Hz. Result:

Zoomed version:

Why "clipping" so badly?
EDIT: When normalized to -3dB before filtering there's no "clipping" present.
Filtering same signal with low pass filter does not have this issue.

That's perfectly normal and expected. You are removing harmonics and as a result you change the time domain wave form. There is no guarantee or expectations that the maximum amplitude will be maintained.
In fact, a square wave has the lowest possible Crest factor (ratio of peak to RMS) and ANY change will increase the Crest factor. In order to not clip you need to remove a substantial amount of power (which your lowpass probably does).


How to apply the CA-CFAR methods in practice? Applying constant false alarm rate is a critical step for estimation of performance of all detectors. Has anybody experience on practical implementation of this estimator for energy detector? Can anybody explain and introduce a code to me?

Constant false alarm rate detection is more of a general technique that underlies a class of approaches; you can't really say that there's "an" implementation of CFAR processing. CFAR techniques are useful when you need to detect the presence of a signal where there is some uncertainty as to what the magnitude of the signal will be, especially relative to any background noise or interference. The general principle of constant false alarm rate detection is that the detector makes some effort to estimate the background level and set a detection threshold adaptively relative to the observed ambient environment. This approach makes the performance of the detector (typically measured quantitatively using the probability of detection $P_d$ and probability of false alarm $P_{fa}$) insensitive to the bulk magnitude of the signal under observation; only relative changes in the statistics of the signal over time (such as a transition from noise only to signal-of-interest plus noise) are important.


Relating bin numbers to frequency non uniform fft I understand how to relate bin numbers of FFT output to frequency when dealing with a uniform FFT. How does this process work with a non-uniform FFT?

The magic of NUFFT is that you create a (oversampled) grid of equally spaced points by performing the interpolation and convolution of your unevenly spaced data. This is mostly done by using Gaussians. After that you apply the Fast Fourier Transform and minimise the effects of the filter you applied previously. Although this method was described by Dutt & Rokhlin in 1993. Currently you can find other, like ones based on Multiresolutional Analysis and Least Square Interpolation.

Nevertheless, you always end up with evenly spaced frequency bins - there are no unevenly spaced frequency bins due to jitter in the time domain.

If you wish to dig into math, then please refer to (among many others):

Particle Simulation Based on Nonequispaced FFT

NUFFT and MRI


How to get velocity from PSD graph Hello everyone!
I have a graph, which is Spectra inside cavity. PSD vs. frequency.
I need to get velocity [m/s] from PSD [dB/Hz].
Does anyone know how to do that?

You can get the sound pressure level from PSD, and the velocity is related to the sound pressure. According to the equation of motion,
$$
\rho \frac{\mathrm{d}\vec{v}}{\mathrm{d}t} = -\nabla p \tag{1}
$$
where $\rho$ is the density of the medium, $\vec{v}$ is the particle velocity, $\nabla=\frac{\partial}{\partial x}\vec{i}+\frac{\partial}{\partial y}\vec{j}+\frac{\partial}{\partial z}\vec{k}$ is the gradient operator, and $p$ is the sound pressure. When the sound is not so loud, the equation of motion can be linearized and the total derivative $\mathrm{d}\vec{v}/\mathrm{d}t$ becomes partial derivative $\partial{\vec{v}}/\partial{t}$.
Furthermore, if you make a plane wave assumption, then the particle velocity amplitude is
$$
v = \pm \frac{p}{\rho_0 c_0} \tag{2}
$$
for plane waves of forward and backward propagation, respectively. $\rho_0$ and $c_0$ are respective the medium density and sound speed without acoustic disturbance. Their product $\rho_0 c_0$ is the characteristic specific acoustic impedance of the medium.
However, your data is acquired in a cavity with an extremely high level of PSD, thus the linearization is obviously inappropriate. The total derivative in Eq. (1) should be
$$
\frac{\mathrm{d}\vec{v}}{\mathrm{d}t} = \frac{\partial \vec{v}}{\partial t} + \frac{\partial\vec{v}}{\partial x} \frac{\partial x}{\partial t}+ \frac{\partial\vec{v}}{\partial y} \frac{\partial y}{\partial t}+ \frac{\partial\vec{v}}{\partial z} \frac{\partial z}{\partial t} \\
= \frac{\partial \vec{v}}{\partial t} + v_x\frac{\partial\vec{v}}{\partial x}+ v_y\frac{\partial\vec{v}}{\partial y}+ v_z\frac{\partial\vec{v}}{\partial z}
\tag{3}
$$
In addition, velocity is also related to frequency, but it seems that you want an overall velocity. You can obtain the time-domain sound pressure according to the Parseval's theorem, and then calculate the particle velocity.


Intuitive explanation of coherence I am a biologist so apologies for the basic question.

I am trying to get an intuitive understanding of what the coherence between 2 signals actually means. I have read a couple of introductory texts and an idea that keeps coming up is that (roughly) two signals are coherent if you can accurately predict the phase of one given the phase of the other (i.e. if there is a constant phase difference).

But I am getting confused when it comes to thinking about coherence with regards to a particular frequency. It feels like any wave at a particular frequency will have a constant phase difference relative to all waves with the same frequency, but from the explanation in the paragraph above that would mean that coherence should be 1 across every frequency for every signal. I am definitely fundamentally misunderstanding this, if someone could please explain where I'm going wrong that would be very much appreciated!

Sorry again for asking something so simple, I am very new to this!

It feels like any wave at a particular frequency will have a constant phase difference relative to all waves with the same frequency, but from the explanation in the paragraph above that would mean that coherence should be 1 across every frequency for every signal.

I think you're falling into the frequency-domain trap.  The frequency domain is nice, but it was invented to use when thinking in the time domain is harder.  If thinking in the time domain is easier -- don't try to wedge things into the frequency domain, unless you have to for some other part of the problem.

So if you have two signals that are predominantly sinusoidal, and if they're both frequency modulated, then you'll have something like $$s_n(t) = \cos\left( \int \omega_n(t) dt\right)$$ where $\omega_n(t)$ is a random process.  If $\omega_1(t)$ exactly equals $\omega_2(t)$ then your signals will be coherent -- even if $\omega_1$ and $\omega_2$ are time-varying.  If $\omega_1(t)$ and $\omega_2(t)$ have equal long-term averages, but small differences at any given time, then your signals will show some coherence (depending on how you assign numbers to "small differences" and "some coherence").  If $\omega_1(t)$ and $\omega_2(t)$ are unrelated, then your signals will have no coherence.

But note: this is really flipping hard to analyze in the frequency domain and not too bad to think about in the time domain.  So thank M. Fourier for all his hard work, do this analysis in the time domain, and then use the frequency domain the next time you really need it.


Sampling a low-frequency sound doesn't work at some sampling rates When sampling a low-frequency sound (like a 55Hz sine wave) I hear an additional high frequency sound at certain sampling rates, but not at others. This happens in Audacity and in my own program. For example, I hear both a low and a high pitch when sampling at 44,100 Hz, but only a low pitch when sampling at 48,000 Hz. What could this be?
More information:
The high frequencies are present for other input frequencies and it's not dependent on amplitude.
I'm sampling the sine wave using the standard sin() function in C++:
for( int i = 0; i < N; i++ )
{
    float time = float(i) / float(samples_per_second);
    output[i] = sin( 2*pi * 55 * time ) * 0.5;
}

Not sure how Audacity is doing it. With Audacity I generate a sine tone, and then changing just the project's sampling rate will make the high pitch appear/disappear.
The same thing happens playing back a bass drum sound (in Audacity and my program). I sample it just by indexing into a table loaded from a .wav file.
Even more information:
I noticed that the high frequency artifact sounds the same as sampling at a different rate with no interpolation, like this:
float in_idx = 0.0f;
for( int out_idx = 0; out_idx < N; out_idx++ )
{
    output[out_idx] = input[ int(in_idx) ];
    in_idx += 44100.0f / 48000.0f;
}

So maybe Windows or the driver or the hardware are deciding not to use interpolation for some reason.

Could be a few things

You may be borderline clipping. Your sine wave has an amplitude of 1, which is just at the edge of clipping (depending on how its rendered). Try it with an amplitude of 0.5
Your hardware is sloppy. For example, cheap laptop sound cards often cut corners in the anti aliasing filters and or clipping management
Your operating system is poorly configured, particular if you use Windows. If you change sample rate in the project, it may just trigger the operating system  to do sample rate conversion. If the sound card only supports 48 kHz natively, the OS will up-convert 44.1 kHz to 48 kHz and the quality of the sample rate converter can be fairly bad. On Windows all audio goes through the kernel mixer (unless you specifically bypass it) and that often messes with the audio as well.
Driver problem. The interface between the OS and the sound card could be sloppy.


Is there such thing as harmony independent timbre characteristic of a sound? I was wondering if any of you know about any technique to extract a harmony (e.g. harmony in a shape of a chromagram) independent timbre characteristic of a sound frame.

What I'm after is a version of spectrum that will make all the piano notes look almost the same despite the fact they occupy different parts of the actual frequency spectrum.

Googling that didn't give me much yet.

Thanks a lot in advance!

EDIT: Great folks in comments pointed out some ambiguity in my question. I'm interested in a metric that is independent of the notes themselves, but captures the whole set of timbral characteristics that could help describe the instrument independently of the note content. It could be a psychoacoustical metric. I'm really curious if anyone knows anything that could fit this task. Thanks a lot!

well, since timbre is a multi-dimensional perceptual phenomenon, there is no single metric.  so lotsa different metrics.

possibly the most common metric is obtained by high-pass filtering the note, measuring the average power, and normalizing that against the power of the unfiltered note and also normalizing against the pitch.  and it yields a "brightness" measure.  the HPF is a digital differentiator.

let your input note be $x[n]$.

differentiator output:
$$ x_1[n] = x[n] - x[n-1] $$

"sliding" average-power of HPF

$$ x_2[n] = \sum\limits_{m=0}^{M-1} (x_1[n-m])^2 \, w[m]$$

$x_2[n] \ge 0$ and $w[m]$ is a window function.

Normalized against the note power (so that louder notes having the same timbre do not get a larger measure of brightness:

$$ x_3[n] = \frac{x_2[n]}{\sum\limits_{m=0}^{M-1} (x[n-m])^2 \, w[m] + \epsilon} $$

$x_3[n] \ge 0$ and you can multiply $x[n]$ by nearly any value you want and get virtually the same $x_3[n]$.  $\epsilon>0$ is necessary to add to the denominator to preclude division by zero.  but $\epsilon \ll |x[n]|^2$ for most non-zero $x[n]$.

now this is normalized against different note loudness, but is not normalized against the pitch of the note.  higher pitched notes will result in higher values of $x_3[n]$, even if the waveform is otherwise identical.

to normalize against the pitch of the note, you will need a pitch detector or a fundamental frequency detector.  that's another whole problem and i have discussed it here on SE.


How should I pre-process a real valued signal in order to use Kay's estimator? I have 100,000 samples of a signal $x[n]$ that was sampled at 20kHz. The data is vibration data  from a rotating machine, and contains a significant spectral component related to the speed of the machine's rotation.

Because the speed of the machine varies over the duration of the sample, using the peak of the FFT does not yield the result I am looking for.

So I want to use estimators such as Kay's estimator that allow short-term estimates, but assume a signal model of:

$x[n] = A \exp(j \omega n + \theta) + z[n]$

where  $n$ = 0... 99,999, $A$ is the amplitude, $\omega$ is the frequency to be estimated, $\theta$ is the initial offset, and $z[n]$ is the complex noise.

However, my signal is real-valued and looks more like:

$x[n] = A \cos(\omega n + \theta) + z_r[n]$

where $z_r$ and $A$ are now real-valued.

How do I transform my real-valued signal into a complex-valued signal, so that I can use Kay's estimator?

The tool for converting real signals into their analytic representation is the Hilbert transform.
Suppose your signal were a projection of some helical rotation with variable amplitude onto the real-time plane, like in the image below.

Source
The Hilbert transform produces such complex signal given its real part. It is a linear transform and it's very easy to do in frequency domain. Without going into its mathematics and derivations too much in depth, the Fourier Transform imaginary part of your signal is the same as the one of your real signal multiplied by $-j$ (rotated by 90 degrees). By symmetrical properties of real signal, you get the following relationship:

All your negative ferquency components become 0.
Your DC component stays the same.
All your positive frequency components double

In Matlab, for example, you would do the following:
a = rand(1,201);

hilbert_a = ifft( [ 1, 2*ones(1,100), zeros(1,100)] .* fft(a) );

or simply use built-in the hilbert function.


Difference between real and complex impulse responses? I want to know the characteristic differences between real and complex impulse responses.

Thanks, JK

Your question might be the result of a misunderstanding. There is only one type of impulse response, and this is a system's response to a Dirac impulse $\delta(t)$ (in continuous time), or - in discrete time - to a discrete impulse $\delta(n)$. For a real-valued system, i.e. a system with real-valued parameters, the impulse response is always real-valued. For a complex-valued system, i.e. a system with complex-valued parameters (e.g. complex-valued filter coefficients), the impulse response is complex-valued. Note, however, that complex-valued systems are just a convenient way to describe certain real systems, just like complex numbers are often more convenient to establish certain mathematical results.


Beat detection and calculation of start time Beat detection
I'm trying to create a beat detection function of a machine producing  in the language R. I hosted the audio file on dropbox.
It was recorded with 16 bit and 16kHz. I started with some filtering. The frequencies can be seen here. Most backgroundnoise is around 100Hz and I don't really have information above 4000Hz, so i used a bandpassfilter from 300Hz to 4000Hz. The filtered audiofile is also hosted on dropbox and the frequencies now look like this.
I then downsampled to 4000Hz to save data and computation time. Finally I used the hilbert function in combination with autocorrelation to get something like this.
With the peak I can calculate that the machine cycle is 70/min, which is correct.
This is the first time I'm doing something like this, so how can I improve the preprocessing? Should I downsample even more to save recources?
This method only looks at the loudness of the signal, disregards the frequency completly. Do you think thats good enough for something like this or can you recommend a method that also includes the frequency of the signal?
Calculation of start time
I also have to detect where the machine begins its production. I cut the signals into 1.5 second packets and labeled them. I used some simple metrics like mean and max of the loundness to cluster them like this. My idea was to check wheter its running or not, then use the packets of the running machine to calculate the cycle time. Then I can loop back the amount of the cycle time and check, if the machine is running at that point or not.
I don't really like that approach though. Do you have a better idea? Are there some better features than mean and max?

Interesting. I would be very tempted to use correlation and given the OP wants to specifically find the occurrence of a first pulse implies that the correlation interval is restrained to be one pulse cycle wide. (Unless we know further information, and specifically how many pulses occur after a start position- if this is known then further processing can be done on the entire "burst"). The first question to determine the applicability of an optimized correlation approach is how correlated each of the pulses are to each other, within a given sequence and from sequence to sequence (repeatability over time). This can be evaluated by cross-correlating various captures to a representative pulse and looking at the statistics on the peak magnitude of the cross correlation function, and how much margin this allows relative to the "noise" correlation everywhere else.  If that were to suggest correlation could be used (by providing a sufficient probability of detection vs probability of false alarm, which such metrics would be clear from this evaluation), I would then determine the mean reference pulse and use that to detect a pulse present or not present condition using a correlation function. This is basically and optimally a matched filter; optimally if the background noise is white. If the noise is not white, then targeted noise rejection filters could also provide significant advantages - this can be determined by capturing the background with no pulses present and evaluating the time and frequency statistics of those captures to see if there are particular details that can be targeted and eliminated (for example 60 Hz spurs from power lines, etc...).
Where a correlation / matched filter approach wouldn't work well is if the content of each pulse in time and frequency is mostly random from pulse to pulse, in which case an optimized receiver would instead consist of a power detector with a pre-selection filter around the dominant frequency content of the pulse prior to detection.  The first cross-correlation testing described above would reveal the optimum strategy based on those statistics. In either case, detecting the first (or any) occurrence limits the observation time for decision to be one pulse interval, which then limits the available SNR to make any such decision for any approach (unless other intelligence is known that make a first pulse unique from any other).


How to use function spectrum in R to perform spectral analysis I intend to perform spectral analysis with R stat:spectrum and then plot power spectrum.

My code is listed here:

  t <- seq(0,1024,by=0.1)
  x <- cos(2*pi*t) + 0.75*sin(2*pi*4*t) + 2*sin(2*pi*6*t)
  x.spec <- spectrum(x,log="no",span=5,plot=FALSE)
  del <- 0.1
  spx <- x.spec$freq/del
  spy <- 2*x.spec$spec
  plot(spy~spx,xlab="frequency",ylab="spectral density",type="l")

the signal has 3 periodic components with 1 Hz, 4 Hz, 6 Hz, but the final spectrum only shows the first 2 frequencies, without 6 Hz. And I don't know why the range of x.spec$freq is [0,0.5] before scaling.

Please tell me my possible concept misunderstanding and guide me how to plot the spectrum correctly. Thank you very much.

Hi: Your sampling frequency is $f_{s} = 10 $ but the frequency of the third component is $f_{3} = 6$.. Therefore, the nyquist criterion for  the maximum frequency is not being met because $2 * 6 = 12$ which is not less than $f_{s} = 10$. I'm new at this but hopefully this is the reason why. Someone will surely correct me if not.


Fourier Series representation of a signal Use the defining equation for the Fourier Series coefficients to evaluate the Fourier Series representation of the following signal:
  $$x(t)=\sum_{m=-\infty}^{+\infty}=(\delta(t-m/3)+\delta(t-2m/3))$$

I calculated $T=2/3$ and $w_0=3\pi$, however, I'm not sure whether $X[k]$ will be $$X[k]=3/2 * \int_{0}^{2/3}(2\delta(t)+\delta(t-1/3)+2\delta(t-2/3))*e^{-jk3\pi t}dt$$or$$X[k]=3/2 * \int_{0}^{2/3}(\delta(t)+\delta(t-1/3)+\delta(t-2/3))*e^{-jk3\pi t}dt$$
And if $X[k]=3/2 * \int_{0}^{2/3}(2\delta(t)+\delta(t-1/3)+2\delta(t-2/3))*e^{-jk3\pi t}dt$, is this following calculation correct (since I never did integral with $\delta$ before)?
$$X[k]=3/2*(2e^{-jk3\pi 0} + e^{-jk3\pi 1/3} + 2e^{-jk3\pi 2/3})=3+3/2*e^{-jk\pi} + 3e^{-jk2\pi}$$

You're getting yourself into unnecessary trouble by choosing the integration limits exactly at those values of $t$ where you have Dirac impulses. Note that you can choose any integration limit as long as you integrate over one period of the given function.

E.g., if you choose some positive $\epsilon$ satisfying $0<\epsilon\le\frac16$ and you integrate from $-\epsilon$ to $\frac13+\epsilon$ then you have the relevant portion of the signal during one period inside the integral, and there's no question as to the scaling factor of the Dirac impulses. The result of the integration will be a real-valued constant for even $k$, and another real-valued constant for odd $k$. I'm sure you can derive the final result yourself.


How to analyze reverb in existing sound recording? I'm wondering if there is a technique by which to analyze an existing sound recording and extract some details of the reverb and echo that are present (either due to natural environment, or added during production).

Are there software packages or libraries that implement that sort of analysis? (Is there even a chance of recovering the audio minus the reverb?)

You may be able to reverse engineer the reverb if you have a segment of relatively clean reverb 'tail'. This will allow you to play around and get some metrics such as decay time and frequency response.

Convolution reverb works by recording an impulse response in a particular location/machine after a burst of noise. This is then convolved with the input signal to apply the reverb.


Discrete Time Convolution Example How did we get this final answer in example in fig attached.
I understand that delta(n) is a DT unit impulse where it is 1 at n = 0 and zero elsewhere.
I understand that u[n] is a unit-step where it is 1 for n>=0 and 0 elsewhere.

How do these two functions equate in final answer below and how did we get each of the final form of the three terms in the final answer? and why is the 3rd term still in u[n-2] form?

Draw the three discrete time signals (shown on the left side) on top of one another. You will see the summation leading to $\delta[n]$ and $\delta[n-1]$ term as follows. For $n=0$, only the first term survives and hence $3\delta[n]$. For the second term, put $n=1$ that makes $3/2+2=7/2$ at index $1$.

Then, obviously third term extends in $u[n-2]$ form as a summation from all three terms at each subsequent index. This is because $n=2$ brings $3/4$ from the first term, $1$ from the second term and $1$ from the third term, totaling to $11/4$ onwards.


What are the measurable factors of image sharpness? I need a way to present data on an essay written about image sharpening algorithm. So I've thought that if I measure a factor of sharpness from image or relative sharpness (between test img and reference img) and another factor like PSR for noise, then I could analyze resultant data with maths and build the essay around that. Is it the right approach to measure image properties?
I'm novice when it comes to photography and graphics so any kind of help is much appreciated.

Some references on image sharpness metrics:

Encoding Visual Sensitivity by MaxPol Convolution Filters for Image Sharpness Assessment, IEEE Transactions on Image Processing, 2019
A Fast Approach for No-Reference Image Sharpness Assessment Based on Maximum Local Variation, IEEE Signal Processing Letters, 2014
Image Sharpness Assessment Based on Local Phase Coherence, IEEE Transactions on Image Processing, 2013
$S_3$: A Spectral and Spatial Measure of Local Perceived Sharpness in Natural Images,  IEEE Transactions on Image Processing, 2012


How to classify accelerometer data? I am trying to detect if a car did accelerate or did brake by using the accelerometer of the iPhone.

In the figure below I plotted the collected data. (To collect the data the phone was laying flat in the car facing the top of the phone to the front).

The issue I am having now is that visually I can clearly see where the car accelerated and where it did brake (see below). But can't figure out how to do that programatically.

What I basically want is to know how often a car did brake or accelerate during a certain time window. (Window size around 10-30 seconds).

Any help would be highly appreciated.

Thanks in advance!

You should consider doing an STFT on the accelerometer signal. This will allow you to visualize the frequency of the signal. If you can determine the frequency of the brakes, you should be able to set a threshold for that specific braking frequency and count each time it is exceeded. You would also be able to measure how long the braking occurs.

How are you processing your signal? Matlab and python have the specgram function, which will allow you to see the actual frequency braking. From there you will be able to determine the brake frequency as well as a threshold.


Dedicated DSP for delayed replay I'm looking for a dedicated DSP chip that will allow sampling and delayed replay of a signal, given the following parameters:

buffer should be integrated
configuration of the DSP chip (if needed) should be possible via host CPU (SPI or serial)
starting and stopping the sampling/replay should be possible via a dedicated pin
reading the samples via the host CPU is NOT required!
sampling should be done with a frequency of at least 2MHz (since a temporal resolution of 1 µs is required)
since this is for replaying a digital signal, a low resolution (even 1 bit per sample) is possible
buffering should allow for 250 ms period

I did not find any stand-alone DSP with this characteristics.

Hard to believe, that you actually looked. You probably don't even need a DSP for that, a normal microcontroller should be sufficient.

For example, the STM32F4-Discovery board comes with an ARM Cortex-M3 that has several ADC channels that can do 2.4MSps each and up to 7.2MSps in interleaved mode. There's a DAC for signal output as well. But since your signal is digital anyways, you can probably sample with just a normal digital input pin. If 1bit per sample is enough, 250ms of a 2MHz sampled signal require 62500bytes which fits into the memory of the ARM quite well.


Looking for an arcsin algorithm Does anyone have a simple algorithm for computing a reasonably accurate arcsine? By "simple" I mean some sort of polynomial that requires <= 5 multiplies per output sample. And by "reasonably accurate" I mean an algo whose error is no more than 10% when the input argument is close to plus or minus one. I searched the web for a while but found nothing immediately useful.

Here's just a polynomial version:

$$
\arcsin(x) = x + \frac{1}{2} \frac{x^3}{3} + \frac{1 \cdot 3}{2 \cdot 4} \frac{x^5}{5} + \frac{1\cdot 3 \cdot 5}{2 \cdot 4 \cdot 6} \frac{x^7}{7}
$$

function y = arcsin_test3(x)
    y = x.*(1+x.*x.*(1/6+ x.*x.*(3/(2*4*5) + x.*x.*((1*3*5)/(2*4*6*7)))))
endfunction

which seems to have five multiplies (assuming you can save the result of x.*x) and three additions.

And the scilab plot is:

Top is scilab's asin vs this one, bottom is the error between the two.

Original Answer

The square root here might be a hassle, but I thought I'd write it up because it looks like fun. :-)

This page  suggests:

  from page 81 of the
  Handbook of Mathematical Functions, by Milton Abramowitz and Irene
  Stegun:

$$
   \arcsin(x) = \pi/2 - \sqrt{1 - x}(a_0 + a_1*x + a_2*x^2 + a_3*x^3),
$$
where
$$
   a_0 = 1.5707288\\
   a_1 = -0.2121144\\
   a_2 = 0.0742610\\
   a_3 = -0.0187293
$$

I've implemented this in scilab and it works OK, except around $x= -1$.  Just reflecting the $0 \le x \le 1$ over to $-1 \le x \le 0$ makes for a much better approximation.

The top plot shows scilab's asin function against the above approximation (in dashed red) against my change in green.

The bottom plot shows the error for my change (plotting that and the original on the same axes means the green looks zero everywhere).

// 25770
function y = arcsin_test(x)
    a0 = 1.5707288
    a1 = -0.2121144
    a2 = 0.0742610
    a3 = -0.0187293

    xx = abs(x)

    y = %pi/2 - sqrt(1-x).*(a0 + a1*x + a2.*x.*x + a3.*x.*x.*x)

endfunction

function y = arcsin_test2(x)
    a0 = 1.5707288
    a1 = -0.2121144
    a2 = 0.0742610
    a3 = -0.0187293

    xx = abs(x)

    y = %pi/2 - sqrt(1-xx).*(a0 + a1*xx + a2.*xx.*xx + a3.*xx.*xx.*xx)

    y = y.*sign(x);
endfunction

x = [-1: .0100001 : 1];

clf
subplot(211)
plot(x,arcsin_test2(x),'g.');
plot(x,arcsin_test(x),'r:');
plot(x,asin(x))
subplot(212)
//plot(x,(arcsin_test(x) - asin(x)),'r:')
plot(x,(arcsin_test2(x) - asin(x)),'g.')