Single Side band Phase noise I want to know why the single-side band phase noise is more customary rather than dual-side band phase noise?(even both side band are not the same)

In addition, Could you explain this sentence according to formula for phase noise measurement problems?

This is where things get ugly:

Can you define what "one-sideband noise power in 1 Hz bandwidth" is?

In my world, it'd just be the PSD of that signal multiplied by 1 Hz. But then, there'd be no need to call it "one-sideband noise power"; it'd just be the noise power density in 1 Hz.

So, "customary" might very much be a "relative" statement. I'd expect that term to only show up in literature that tries to explain signal properties based on a real-valued, analog technological background (maybe trying to bring the complex signal's view to an audience of people that learned the real-valued signal theory first).
 Frankly, I find it an elegant move that much of the modern literature starts right off with complex signals, and just shows that a simple multiplication with a complex sinusoid can move it to baseband, and then introduces the "normal" (as in: non-quadrature) mixer based on that understanding.

This is probably where that "upper and lower sidebands overlap" comes from:

If you take a noisy bandwidth $B$ of passband around a carrier at $f_0$ (you observe $]f_0-\frac B2; f_0 + \frac BS]$), and mix it down using a quadrature mixer (giving you I and Q signals, which we commonly interpret as real and imaginary part of a complex-valued signal), you get signal in $]-\frac B2; \frac B2]$.
The PSD of noise in that range then is the same as the PSD of noise in the passband. That happens because, mathematically, the quadrature mixer just multiplies with $e^{-j2\pi f_0}$, leading to a "clean" shift of the positive frequencies down to baseband (by $-f_0$).
If you, however, just multiply with a $\cos(2\pi f_0)$, you'll mathematically shift the positive to baseband, but also you shift the negative frequencies, which are a mirror of the positive frequencies, up (by $+f_0$).
Now you have the upper sideband being overlaid with the mirror of the lower sideband.
So, looking at this differently, this is an effect of using real-valued signals to represent something that would become complex baseband when mixed with a quadrature mixer.

Note that the statement is obsolete (or plain wrong), because frankly, unless you really want this for some reason, measurement devices either

mix onto an intermediate frequency (IF), so that none of the signal "crosses" 0 Hz when being downconverted, or
employ quadrature mixers, where that overlap doesn't happen.

Update I'll have to revise my above statement about the obsolence of the statement, because your book,

  Phase Noise and Frequency Stability in Oscillators by Enrico Rubiola

is a book on mixer theory, where you of course analyze the workings of the elementary mixers (which can be combined to form e.g. quadrature mixers), so that the author inherently works with the more "physical" real-valued signals, and then he's obviously right:

if you're characterizing the mixing element itself, you can't tell noise above and below the LO apart.


Mathematical proof of relationship between derivative and convolution What I'm looking for is the mathematical proof of why we can calculate derivatives as a simple convolution of a mask, and how do we get that mask?
I know it has something to do with how the derivative is approximated

Can anyone direct me to some book or anything?

There's no mathematical proof per se; the fact is that finite differences are one method of approximating derivatives. Remember the definition of a derivative:

$$
f'(x) = \lim_{\Delta x \to 0} \frac{f(x+\Delta x) - f(x)}{\Delta x}
$$

The above equation states that to calculate the derivative of a function $f$ at position $x$, one takes the difference between the function evaluated at $x$ and the function evaluated at some small distance $\Delta x$ away. The ratio of the difference in the function value between the two points and the distance traveled along the $x$-axis is defined to be the derivative, as the distance between the two points becomes infinitesimally small.

A finite-difference derivative method, then, is just a way of approximating the above relationship using some finite spacing instead of the infinitely-small spacing implied by the limit. This type of approach makes sense for real-world applications where you might not analytically know the function $f(x)$  (and therefore can't evaluate its behavior across arbitrarily-small intervals.

With that said, there are a number of ways to approximate the derivative of a signal in this way; one example is the causal first-order difference:

$$
x'(t) \approx \frac{x(t) - x(t-T)}{T}
$$

This is simple to implement, but performs best if $T$ is small relative to the inverse of the signal's bandwidth. A slightly more complicated method can provide better results in some cases (the second-order difference):

$$
x'(t) \approx \frac{x(t) - 2x(t-T) +x(t-2T)}{T^2}
$$

I've used the "backward second-order difference" form above to show how they could be applied using causal filters, but note that they will have a resulting delay commensurate with the order of the approximation.


Variance of the Discrete Gabor Transform for practical implementation When we design a discrete Gabor filter for a 1D signal, how do we determine the variance of the filter, depending on the time step of the signal? Is there any empirical theorem/result that we should know about the Discrete Gabor Transform?

What I look for in such parameter choice problems is empirical experiments. Whatever works best for my application is the correct value to use.

You can develop a parameter search algorithm looking for the optimal values in some training data set. Keep in mind that what works for one application may not perform well in others.

I'm not aware of mathematical arguments estimating  gabor's $\theta$. These estimates are typically conservative anyway.


Error correcting code? - using RGB space with a lookup table This problem could be thought about as how to transmit a signal with 3x8bit RGB in a way that is robust to the distortion of antialiasing blends colors.
The context: I need to rewrite a short utility library, to get it working with the Brave browser (My actual question isn't about brave per se.)
canvas-color-tracker - example of it being used and src/index.js  is short
The HTML canvas is used for mouse over and clicks to determine what object is under the cursor.
To use it with circles, you register your circle to get a color, then draw the circle in that color. Then to determine what is under the cursor, you get the color with getImageData.
The trouble is what to do about anti-aliasing, the browser doesn't allow that to be turned off, and I want to be able to return nothing rather than collide with another uninvolved circle when it's edges blend into another color)

Assuming the anti-aliasing blends linearly How to create an ordering of points in 8bit RGB space, such that the minimum number of points are colinear without knowing the final total number of points? What is a reasonable heuristic for this?

Addendum, I am doing this because Brave will choose a color channel at random and randomly vary the lowest bit of it - bit fiddling as a privacy feature prevent fingerprinting.

This is a bit of a "I hope I understand your question right" answer, but seeing I have to leave for the day:
Randomizing the LSBs of all the three color channels would obviously reduce the amount of information a single pixel could carry would reduce the amount of information carried in that color from 24 to 21 bits.
Picking one randomly and randomly varying that has an entropy $1+ \log_2 3$ bits, so you reduce the information contained in that pixel's color from 24 to ca 21.415 bit. This doesn't look like it seriously introduces any privacy there.
It will be impossible to find a 24-bit code that has a rate higher than 21.415/24 that can correct 3 bit errors, so, honestly, you usually simply wouldn't need to use error correction at all.
Still, I think they're doing that for a reason – probably it breaks an existing system of identification. There's no "provable" privacy to be won here, just an "against a single implementation" privacy.

Assuming the anti-aliasing blends linearly

Not the case, usually; you use low-pass filters that are more than a linear interpolation to get acceptably-looking edges.

How to create an ordering of points in 8bit RGB space, such that the minimum number of points are colinear without knowing the final total number of points?

Oh, that's a very interesting problem!
However, we do know from linear algebra how to find linearly independent vectors in a vectors space (e.g. simply choosing canonical base vectors that are all zero but for one entry, or taking a set of random vectors and doing Gram-Schmidt on them).
In this case, our vector space is 3-dimensional. In a three-dimensional space, you get exactly 3 bases in each orthogonal base system, soooo, that's three code points that are orthogonal only! That's very few.
Now, if you add two of these code points, you get another code point, which isn't colinear with any of the existing code points, and since 3-choose-2 is 3, you get three new code points, for a total of 6. You can also add up all three original bases and get one new. OK, so 7 code points.
You can iteratively repeat that, but you have to be careful: You will find multiple ways to calculate the same code point later on, or even multiples of that, and at some point, your elements' coordinates can get larger than your 8-bit integers can hold. But since essentially computers are fast, and the number of things you'll get will be still quite low, you can probably really just brute force your way out of this reasonably by going through all possible pairs of points you have, adding them, and checking that the result is not a multiple of any existing point.


Verb for subtracting the mean? Is there a verb for removing the mean or DC component of a signal? I am trying to name a flag for doing just that.

"Normalize" can mean many things, so I don't think that would be a good choice.

I've heard this referred to as demeaning a signal, or applying a DC notch (or, in some more hardware-related domains, a DC block). More generally, it might be called just highpass filtering the signal to remove the DC component.


Why does this FFT based signal deliver a sinusoidal result? I have noticed that there is an identical pattern between two signals:

one is a signal of n identically spaced samples over a single period of a sinusoid

the other is a signal constructed of the real component of the second (i.e. first periodic) Fourier bin of n circular shifts of the signal, so pseudocode would be:

allocate a final result, size n
loop from 1 to n
    newSignal = circularly shift signal one step to the left
    newValue = take the FFT, second (first periodic) bin, real component
    finalResult(n) = newValue

The only difference between these two is the scaling. If one's input is known to be a single cycle of a sinusoid, this is a nice way to remove all noise without knowing the phase of the input. Obviously it works in part because the second Fourier bin corresponds to the frequency of the input. But that's as far as my intuition gets me. In terms of theory, I don't understand why a signal composed of the real component of the second Fourier bins of the circular shift of each step of the signal looks the same as the original signal. Can anyone shed insight?

Let me see if I understood you correctly. You have one period of a sinusoidal signal:

$$x[n]=\sin\left (\frac{2\pi n}{N}\right ),\quad n=0,1,\ldots N-1$$

Its DFT is given by

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

For $k=1$ we have

$$X[1]=\sum_{n=0}^{N-1}x[n]e^{-j2\pi n/N}$$

For the given signal $x[n]$ we get

$$X[1]=\frac{N}{2j}$$

For the left-shifted versions of $x[n]$ we get

$$X_l[1]=\sum_{n=0}^{N-1}x[n+l]e^{-j2\pi n/N}=e^{j2\pi l/N}\sum_{n=0}^{N-1}x[n]e^{-j2\pi n/N}=e^{j2\pi l/N}\frac{N}{2j}$$

So we have

$$\Re\{X_l[1]\}=\Re\left\{e^{j2\pi l/N}\frac{N}{2j}\right\}=\frac{N}{2}\sin\left (\frac{2\pi l}{N}\right )=\frac{N}{2}x[l],\quad l=0,1,\ldots, N-1$$

Q.E.D.


Logic for bit growth in radix-2 fft In radix-2 FFT, why is the bit growth equal to $1+\sqrt{2}$ per stage?

I want to know the logic as to how that value has been obtained. Is it just empirical or has it been derived?

FFT error analysis is a well-studied field with a significant literature.  For a radix-2 Decimation-in-Time (DIT) butterfly, bit growth per stage can be 2.4142..., and for a Decimation-in-Frequency (DIF) stage it can be 2.8284... .  Other factors are related to: matter of rounding (eg: round up, round down, round to nearest, truncate, etc.), hardware considerations (number of guard bits, overflow/underflow detection/correction, etc.), algorithm type, number representation, etc., etc.

Although the maximum numbers from a butterfly can be determined by inspection, significant analytical work has also been done in expressing average errors, error variance, etc., subject to the conditions mentioned above.

There's other info available, such as: https://groups.google.com/forum/#!msg/comp.dsp/c5c9DKYE1Qw/8SNY4gRqqBcJ


How can I show DTFT result in MATLAB? I wanna show DTFT result and convolution result in f-domain are same.

Additionally, I wanna show sampled function and inverse fourier transform result are same too.

How can I show this? If I use just FFT or DFT, they don't show picture-like.
If I have to use other function (e.g diric, freqz, etc), then please let me know how can I use, and input value mean.

The DTFT is continuous, so you have to approximate the DTFT on a computer.

Zero padding “interpolates” the DFT, so  a lot of zero padding will approximate the DTFT


Detecting Bass in a wav-file In this code, i am doing a STFT on my audio file. I want to find the rhythm of file. Bascially the bass, which is located in the lower frequencies. I think i need to find the BPM of my file. The tempo of a song is measured by BPM, which is i need.
There is a method,that if there's strong beat is to take the magnitude of the STFT of the waveform, and then auto-correlate it in only the time dimension. The peak of the auto-correlation function will be the beat, or a submultiple of it.

This is equivalent to breaking up the signal into a lot of different frequency bands, finding the amplitude envelope of each, autocorrelating each envelope, and then summing them. The noise and other parts of music are averaged out by the cross-correlation operation.

I tried to find this method, how excatly this can be performed. I was searching for an open-source-code or any kind of instructions, how this can be done. But i wasn't able to find anything. Does anyone can help me to reach the steps explained in the methods?

%% MATLAB
%% STFT = SHORT TIME FOURIER-TRANSFORMATION
%___________________________________________
clear all;
[y,fs]=audioread('UnchainMyHeart.wav');
% audioread = Reading the WAV-File
% y = A vector, which contains our audio signal
% fs = Sampling frequency
% 'UnchainMyHeart' = Name of the WAV-File

%PARAMETER FOR THE CODE
%_____________________________________________
t_seg=0.09; %Length of the segment, on which we use the STFT

fftlen = 4096;
% Length of the fft

% L,M,H bands
% Start & Stop frequency defenition of bands
fL = [  1,  300 ]; % Low frequencies
fM = [  301,  5000 ]; % Medium frequencies
fH = [ 5001, 22050 ]; % High frequencies
%_____________________________________________

segl =floor(t_seg*fs);
% Length of the segment(50ms), which is being multiplied by the
% sampling frequency
% The result is rounded off with the function "floor"

windowshift=segl/3;
% Defining the size of the window, which goes to next segment and so
% on,until to the end of the audio signal

window=hann(segl);
% Hanning Function, which is being stored in a vraiable
% Matlab usually works better with variables than actual numbers. Hence
% the workspace. This way, you can avoid errors.

window=window.';
% transpose vector from a line vector to a row vector
% In the workspace, you can see, what kind of vector it is

si=1;
% Start index

ei=segl;
% End index

AOS=floor( length(y)/windowshift - 1);
% Determining the numbers of segments in my audio signal

f1=figure;
% New figure is being openend
n=0:1:fftlen-1;
f=n/(fftlen-1)*fs;
% Definig frequency vector

Ya=zeros(1,fftlen);
% The values between "1" and "fftlen" are being filled with zeros

for m= 1:1:AOS

y_a = y(si:ei);
y_a= y_a.*window;
Ya=fft(y_a, fftlen);
YA=abs(Ya);

drawnow;
% Updates the graphical objects

% indices of bands
LI = find((f>=fL(1))&(f<=fL(2)));
MI = find((f>=fM(1))&(f<=fM(2)));
HI = find((f>=fH(1))&(f<=fH(2)));
% mean values of bands
ML = mean(YA(LI));
MM = mean(YA(MI));
MH = mean(YA(HI));

figure(f1);
plot(f, 20*log10(YA));
ylim([-90 50]);
xlim([0 fs/2]);

% Limiting the length of my y-axis
title('Spectrum of audio signal');
% Title
xlabel('f(Hz)');
% Name of x-axis
ylabel('abs(Ya)/dB');
% Name of y-axis
grid on;
% Generating grid raster

si=si+windowshift;
%Start index is being updated
ei=ei+windowshift;
%End index is being updated

end

You should be more clear about Detecting Bass.

Bass instrument detection (i.e. polyphonic pitch tracking) and Bass Drum detection (to compute "rhythmic features" or tempo/BPM) are two different tasks. I'll assume you're talking about the latter.

I recommend looking at the open source music-analysis library called librosa, you can do tempo analysis in one line, on audio you've read into memory. You can also use the computed beats to compute various other beat-synchronous features:

import librosa
y, sr = librosa.load(librosa.util.example_audio_file(),sr=None)
tempo, beats = librosa.beat.beat_track(y=y, sr=sr)
print(tempo)
print(beats)

Links to digest:

https://github.com/librosa/librosa/blob/master/examples/LibROSA%20demo.ipynb
http://librosa.github.io/librosa_gallery/


Implementing Frequency Domain Convolution in MATLAB (Convolution Theorem) I'm trying to prove convolution in time domain is same as multiplication in frequency domain but I'm not getting the same answer in matlab.

Here is the code:

You need to set the length as well in your fft command.
The two signals are of length $5$ and their convolution is of length $5+5-1=9$. So use this instead:

a=[1 2 3 4 5];
b=[6 7 8 9 10];
x=fft(a,9);
y=fft(b,9);
ifft(x.*y)

ans =

Columns 1 through 8

6.0000   19.0000   40.0000   70.0000  110.0000  114.0000  106.0000   85.0000

Column 9

50.0000

conv(a,b)
ans =

 6    19    40    70   110   114   106    85    50


period detection for a series of photons I have a dataset which records the arrival time for a series of photons.
I am using lomb-scargle method of scipy.signal to try to find some kind of period.

Since the arrival time is random,how to express the intensity of the signal?
If I bin the dataset and let the photon counts in every bins be the intensity , the dataset changes to a evenly spaced new dataset and information is lost.

lomb-scargle is good at finding period for unevenly spaced data.
Is it necessary to bin the data? If I bin the data, part of the data info will be lost(?) and I can use fft directly.

As discussed in the comments to your question the recorded data that consists of arrival times of photons must first be converted to photon counts $n$ for this is the light intensity. The result of this conversion will be a vector $n(k)$ (with discrete time $k$) that represents the evolution of light intensity over time and that can be analyzed in more detail with an FFT, for example. The Lomb-Scargle method won't help here as you'd like to analyse light intensity, not photon arrivals. It could be of use if your data was a sequence of unequally space light intensity values. But as I understand your question this is not the case.

In order to obtain $n(k)$ count the number of photon arrivals in every time interval $t_0 + kT$, where $t_0$ is some constant time offset (e.g. first arrival time) and $T$ is the length of the equidistant time intervals. $T$ determines the sampling frequency $f_\mathrm{s} = 1/T$ of $n(k)$ and consequently the range of any spectral analysis, that is $0 \ldots f_\mathrm{s}/2$. On the other hand $T$ should be sufficiently large to obtain a reliable photon count in each time interval.

To analyse the periodicity of $n(k)$ calculate its FFT and look for distinct peaks. Functions like periodogram also come handy sometimes.


How to calculate phase for discrete time complex valued signal I am studying Digital Signal Processing course for my Computer Science and Engineering major.

I have to calculate magnitude and phase of following
$$x(n) = e^{(-0.1+0.3j)n}$$
where
$$-10\le n \le 10$$

Here the "e" is Euler's constant. I can calculate the magnitude for an n, by putting a value of n from range and calculating value of imaginary and real part. But how do i calculate the phase? My concept is not clear on this. Kindly explain in detail and as easily as possible.

For a complex number expressed in the polar form, $x = r e^{j\phi}$ , the magnitude is $$|x| = |r e^{j\phi}| = |r| |e^{j\phi}| = |r|$$ where $r$ is a real number and the phase is $\phi$ in radians.

When the complex number is expressed in rectangular form $$x = a + j b =r \cos(\phi) + j r \sin(\phi)$$ then the angle is expressed as $\phi = \arctan(b/a)$

For a sequence of complex numbers $$x(n) = e^{(-0.1+0.3j)n}$$ for each n,
one can find its magnitude and phase once writing it in the polar or rectangular form where the polar form is $x[n] = ( e^{-0.1 n}) (e^{j 0.3 n})$.

Then its magnitude is $|x[n]| = |e^{-0.1 n}| |e^{j 0.3 n}| = e^{-0.1 n}$ and its phase is $\phi = 0.3 n$ radians for each n.

Note that principal value of the phase is the angle beetween $-\pi$ and $\pi$, any value other than this range can be equivalently mapped into it.


How do I make sense of the cosine wave having Fourier Transform coefficients which have infinite magnitude? To illustrate my question better, consider the Fourier Transform of an aperiodic (as a periodic cosine wave has a Fourier Transform not Fourier Series) cosine wave

$$f(x) =
\begin{cases}
\cos(2\pi f_0x),  & \text{$ 0 \le x \le 2\pi$} \\
0, & \text{else}
\end{cases}$$

$$F[f(x)] = \frac{1}{2}(\delta(f - f_0) + \delta(f + f_0)) $$

This above implies an aperiodic cosine wave consists of two Fourier Transform coefficients which have infinite magnitude (as a delta function has infinite magnitude). Therefore, does this imply that a perfect cosine wave can not be produced in real life?

Your function is both integrable, and square integrable, because it is continuous on the finite support $[0\;2\pi]$. In this context ($L_1 \cap L_2$), its Fourier transform is pretty well defined as a classical function (not a distribution like the Dirac), and continuous. So, your formula is incorrect.

So there is no infinite magnitude. Moreover, one cannot simply say that Dirac have "infinite magnitude".

Yet, this would not imply that a perfect cosine wave can not be produced in real life; produced by whom? However, there is little chance that somebody could observe, or capture, a perfect cosine wave in real life (caveat: if real life can be described by a continuous time variable, and amplitude could be continuous as well).


Identifying the magnitude and impulse response from pole zero plot quickly I have an exam next week and it's verty certain that a task of this kind will be there.

Are there some good tips how to match the right pole zero plot to the right responses? No proof is needed in the task.

At the moment I try to find two of the same kind and kinda guess which of the two same looking responses could be which.

What I know so far for the impulse response:

Amount of zeros = amount of poles -> dirac impulse at 0
Complex conjugated poles cause oscillation

For the magnitude response:

Zeros in the left half plane have low pass characteristics/highpass charactresistics in the right half plane
Poles outside the unit circle means unstable system

Are there quick ways without computing $\left\{\left|H(e^{j \Omega } )\right|\right\}$ ?

For A:
zero at z=1, and complex poles, so it has to be one of the oscillating impulse responses: 1 or 5. D has the poles closer together so it oscialltes less, so D is probably 1 and A is 5.

For B
Zeros are symmetric so magnitude response is probably symmetric.

For C
Amount of zeros = amount of poles, so dirac at 0 (2 or 3)

For D
See A

For E
No idea, should have a dirac at 0, but no oscillation?

Really unsure about the magnitude responses.

Some of the cues you want to look for are the following:

A zero on the unit circle causes a zero in the magnitude response at the corresponding frequency ($z=1$ is DC, and $z=-1$ is Nyquist)
For a zero at DC ($z=1$), the sum of all impulse response coefficients must be zero.
A zero near the unit circle causes a dip in the magnitude response at the corresponding frequency. The closer the zero to the circle the more pronounced the dip.
A pole close to the unit circle causes a peak in the magnitude response at the corresponding frequency. The closer to the circle the sharper the peak.
Complex conjugate poles cause oscillation. The closer to the unit circle they are, the slower does the impulse response decay.
If all poles are at the origin, the impulse response has a finite length (FIR).
If there are fewer (finite) zeros than poles then the first few samples of the impulse response are zero. The number of zero coefficients equals the difference between the number of poles and (finite) zeros.

For PZ-plot A, the magnitude response must be zero at DC and it must peak at $\pi/2$ (i.e., at $0.5$ in the magnitude plot). The impulse response (IR) must be oscillating and its values must add up to zero (which you can of course check only approximately). Furthermore, since there are two poles and only one zero, the first IR coefficient must be zero. So we are left with IRs 1 and 5, and with magnitudes b and c. Since we have a PZ-plot which is qualitatively the same as A (namely D), just with a smaller pole radius, we can easily guess that the peakier magnitude and the more slowly decaying IR correspond to A, whereas the others correspond to D. So we get

  A-5-b

  D-1-c

PZ-plot B must be FIR (all poles in the origin), and the first coefficients must be zero (3 poles, only 2 zeros). The magnitude must be symmetrical with respect to $\pi/2$ ($0.5$ in the magnitude plot), and there must be (equal) dips at DC and at Nyquist. Only magnitue plot "a" satisfies these criteria. The only matching finite IR is "4". So we have

  B-4-a

PZ-plot C is also FIR, and the first coefficient is not zero (because the number of poles equals the number of zeros). The magnitude must have a slight dip at $0.5$ and it must be symmetrical (because the PZ-plot is). This gives

  C-2-e

This would leave us with

  E-3-d

but the latter is non-sense, as I've mentioned in my comments. The PZ-plot E corresponds to a finite IR, but that IR must be complex valued. From the plot, let's estimate the zero at $0.8j$, which gives us

$$H(z)=c\cdot (1-0.8jz^{-1})\tag{1}$$

with an arbitrary scaling constant $c\neq 0$. This corresponds to an impulse response

$$h[n]=c(\delta[n]-0.8j\delta[n-1])\tag{2}$$

which is obviously complex-valued. The magnitude response must have a dip at 0.5, but no dip at -0.5 (note that for real-valued impulse responses the magnitude response is symmetrical). IR 3 is just an impulse at zero, so the corresponding PZ-plot must be empty, because there are neither poles nor zeros. IR 3 just shows the real part of the actual impulse response. The magnitude response "d" can be the response corresponding to E, but only for positive frequencies, which are not sufficient to characterize the system. So the PZ-plot E is misleading, to say the least. The same is true for the single impulse shown in 3, because it corresponds to none of the PZ-plots and to none of the magnitude plots shown in the figure.


Why do envelope extraction before pitch detection? I am referring to this report on pitch estimation.

The author does envelope extraction using Hilbert transform before applying AMDF.

What is the advantage of doing that?

Is speech a narrowband or wideband signal? I read somewhere that hilbert envelope extraction works better for narrowband signals.

When I extracted envelope of a speech signal using Hilbert transform, in some places the envelope was almost the same as the original waveform(Figure 1), whereas in other places it tracked the upper portions of waveforms (Figure 2).

Red curve is the Hilbert envelope, blue curve is the speech waveform.

Fig 1

Fig 2

My guess is that the Hilbert Transform is a pre-processing step which attempts to make the actual pitch detection stage (AMDF, in this case) more robust.

The basic reasoning is this: pitch detection of speech is made more difficult by the time-varying spectral characteristics of the signal.  On the other hand, pitch detection of an impulse train is easy: just count the number of samples between each pulse (or look for the first peak of the auto-correlation function/the first dip of the AMDF etc). Pre-processing the speech signal to make it more pulse-like therefore improves the results.

In your diagrams, the Hilbert transform does make the speech pulses more impulsive, though I think that other pre-processing steps (e.g. inverse-filtering by the estimated spectral envelope) is likely to yield better results.


FFT starting at zero signal processing is not my domain of expertise , but I need to make use of it, to do an HRV analysis. These are steps I do starting from my raw ECG signal.

ECG signal -> RR peak detection -> find RR intervals -> cubic spline interpolation -> IIR high pass filter -> FFT

The peaks I get is a mirrored image, I guess its because its a pure signal from a digital patient simulator like all at 60 bpm or all at 120 bpm etc and follows the rule of Hermetian symmetry.

However my FFT output at pure 60 bpm should have indicated a peak at 1 Hz and a 120bpm at 2 Hz, is what I understood. But my all peaks start from 0.

Please help

PS : Implementing this all in an android application.

PS : This is what my output looks on a 60 bpm pure signal.

PS : Added log data

FFT output : http://www.yourfilelink.com/get.php?fid=1337132

The FFT output will be bins evenly spaced from DC (bin 0) to 1 bin less than your sampling frequency. Thus your sampling frequency is at N for bins 0 to N-1, where N is the number of samples in your time domain sequence.

Within the first Nyquist zone of $-f_s/2$ to $+f_s/2$ where $f_s$ is the sampling rate, the bins from 0 to N/2-1 are the positive frequency components, while the bins from N/2 to N-1 are the negative components. Specifically bin N-1 is -1, bin N-2 is -2, etc so you could view the domain as a circle from 0 to N-1 and repeating again.

The reason your peaks are a mirror image is because your time domain signal is real, and for real signals the positive and negative frequencies have a conjugate symmetric relationship (same magnitude, opposite phase).

The reason you see a large response in bin 0 is because your time domain signal has a non-zero mean. Bin 0, being the DC term, is the mean of your signal (scaled by a factor of N). You likely don't care about DC, and would therefore benefit by removing the mean of your signal before taking the FFT.

The FFT is an algorithm that implements the DFT efficiently, which is as the name implies a discrete version of the Fourier Transform. I think it is helpful to see that the Fourier Transform is simply a correlation as we sweep frequency. (Where correlation is complex conjugate multiply and integrate, or equivalently discrete as complex conjugate multiply and accumulate). Thus for the DFT, each bin is a correlation of a single rotating phasor (your traditional sine wave is made of two equal and opposite rotating phasors as shown by Euler's Identity: $\frac{1}{T}(e^{j\omega t}+e^{-j\omega t}$).  So for the DFT, given as: $$X[k]=\sum_{n=0}^{N-1}x[n]e^{-j(k\omega_0) n}$$

So we see that for the first bin (k=0), we just have N times the mean of x. If x has a DC value, we would get a strong correlation for that calculation. Also observe when k =1, which corresponds to the next bin, the exponential is spinnning at its lowest possible rate $\omega_0$ corresponding to one cycle over your time domain sequence. So in that calculation we are computing the correlation of the lowest possible frequency that is spinning in the opposite direction ($e^{+j\omega t}$): such a signal in our x(t), after this computation would have the same result as bin 0: we would get N times the mean after having "derotated" the signal that was present there. Thus every bin goes through the same "averaging filter" which is a Sinc function in frequency (not a great filter), so our output at any given bin will be sensitive to some degree of frequencies everywhere in our spectrum (and as we see with the Sinc function there will be discrete frequencies that are at null positions and therefore will be locations where there is no "leakage" into other bins).

This should give you enough insight into some important fundamental points to dig further into the DFT/FFT. In addition here are some other related posts I think will help you:

Frequency measurement from acoustic emission sensor data(VDC RMS)

What proportion of a padded FFT should be actual values

Signal Processing/FFT gives very high magnitudes for low frequencies

You will want to read up a little more about windowing too; the links above and other posts (by search) include more information on that. Ultimately fred harris' paper has been a great resource for me on that topic.


How to apply the gain factor G to SOS coefficients when designing a IIR filter on sptool (MATLAB)? After designing an IIR filter on sptool, I export the coefficients (SOS) and gain (G) into a .mat file. When filtering, we can use sosfilt(SOS, x), which results in the correct output but with the wrong scaling. I guess this is what the gain variable G is for, but the documentation of sosfilt() never mentions this.
The impulse response of the filter has values of several powers of ten (~10^16), whereas in sptool it had very small coefficients (< 1).
How is one supposed to use the SOS and G variables when exported from sptool, to correctly scale the weights of the exported IIR filter?

The easiest way would be
y = prod(G)*sosfilt(SOS,x);

An alternative would be put to the cumulative gain into the first section.
sosScaled = SOS;
sosScaled(1,1:3) = sosScaled(1,1:3)*prod(G);
y = sosfilt(sosScaled,x);


How does one interpret an element of the "transfer matrix" used to calculate frequency domain granger causality (via VAR models)? I am attempting to gain a better mathematical understanding for how autoregressive models can be used to infer frequency-domain granger causality. All freq. domain measures of causality that utilize autoregressive models convert the time-domain parameters to a frequency-based representation (i.e. take the Fourier transform of the original parameters), and then take the inverse of that representation to get the so-called "transfer matrix". This transfer matrix is referenced in many causality-related papers, and is the backbone of all freq. domain measures. Here is the derivation for the transfer matrix:

$a_{i,j}(k)$ - Autoregressive coefficient describing influence of node $j$ on node $i$ at lag=$k$

$p$ - Model-order (i.e. how many past time-points to use to predict the present)

$\textbf{A}(k)$ - $Q$-by-$Q$ autoregressive parameter matrix at lag=$k$. This matrix is estimated via model fitting.

$U(t)$ - Innovation noise (i.e. the difference between the model's predictions and observed data at time $t$.)

\begin{align}
   {Y(t)} = \sum_{k=1}^{p}\textbf{A}(k)Y(t-k) + U(t)
\end{align}

$\bar{\textbf{A}}(f)$ - Frequency representation of $\textbf{A}(k)$

$\bar{{a}}_{i,j}(f)$ - Frequency representation of autoregression describing influence of node $j$ on node $i$ at freqeuncy $f$

$\bar{a}_{i,j}(f) =
\begin{cases}
                                   1-\sum_{k=1}^{p} a_{i,j}(k) e^{-\textbf{j}2\pi f k} & \text{if $i=j$} \\
                                   -\sum_{k=1}^{p} a_{i,j}(k) e^{-\textbf{j}2\pi f k)} & \text{otherwise}
  \end{cases}
$

\begin{align}
    \text{where } \textbf{j}=\sqrt{-1}
\end{align}

Note that above piecewise equation stems from taking the FT of model coefs. At least as far as I understand it. I'm not actually positive as to why this equation is piecewise and not identical for all $i$/$j$ combinations. I haven't worked with FTs very much.

$\textbf{H}(f)$ - Transfer matrix

\begin{align}
    \textbf{H}(f)=\bar{\textbf{A}}^{-1}(f)
\end{align}

Now for my question--does anyone have an intuitive definition for what the transfer matrix, or a given element of the transfer matrix ($h_{i,j}(f)$), actually represents? I can easily and intuitively understand what $a_{i,j}(k)$ represents, but I've yet to fully grasp what $h_{i,j}(f)$ represents. How does one interpret such a value in layman's terms?

Hi: At a quick glance,  the Geweke, 1982 and 1984 references on the last page of this link are probably necessary reading for understanding it.

I realize that this doesn't answer your question but I'm trying to see if I'm able to send a link without actually typing the link.
and I don't know how to do that in a comment. Hopefully it works here.


Converting PLL Discriminator into Doppler error I'm writing a GPS receiver SDR software, and I have calculated the PLL discriminator ouptut as atan(Q_Prompt,I_Prompt)
I would like to know the formula that converts the discrminator output to doppler error

Frequency is the derivative of phase versus time, so you can estimate the frequency by calculating the difference between the phase of two successive samples and dividing by the sample time which is the radian frequency in radians/sec- so divide by $2\pi$ to get Hz.

For efficient receiver implementation instead of computing the arctangent and performing the computations above consider using the cross product frequency discriminator given by the imaginary term of the complex conjugate multiplication of the current I,Q symbol with the previous, which is beautiful as it simplifies to:

$$Imag((I_1+jQ_1)(I_2-jQ_2))= j(I_2Q_1+Q_2I_1)$$

Which is proportional to the phase difference between the two symbols and therefore the proportional frequency error term. This would provide the frequency error for a frequency lock loop.  Similarly by comparing it to a sample on the real axis we see that the phase error for a phase lock loop is simply proportional to Q for each sample.

The I,Q samples above would be the prompt correlation samples and the complex vector I+jQ should be normalized to constant magnitude for the computations above.


How is the matrix $\mathbf R_x$ not Toeplitz in case of a signal missing one term? I am solving a question that says if we have sequence $x(n)$ of a signal missing one term then we have to find autocorrelation matrix $R_x$ as follows:

$$R_x = E\{\mathbf {xx^H}\}$$
Now if I take $x(n) = [x(0) , x(2),x(3)]$ that is missing term $x(1)$
I found following value of $R_x$

$R_x =  $$\left(
        \begin{matrix}
        r_x(0) & r_x^*(2) & r_x^*(3) \\
        r_x(2) & r_x(0) & r_x^*(1) \\
        r_x(3) & r_x(1) & r_x(0) \\
        \end{matrix} \right)
$

Now I check that this matrix is Toeplitz or not, by the definition of Toeplitz (constant diagonal term) $R_x$ is Toeplitz. But when I checked the solution it says the above matrix is not Toeplitz and there is no explanation given. I wonder what I am missing in definition of a Toeplitz matrix that does not fit in the above matrix. Can anyone explain this?

For $\boldsymbol{R}_x$ to be toeplitz it should be:

$$\boldsymbol{R}_x = \begin{bmatrix}r_{x}(0) & r_{x}(2) & r_{x}(3) \\ r_{x}(2) & r_{x}(0) & r_{x}(2) \\
r_{x}(3) & r_{x}(2) & r_{x}(0) \end{bmatrix} $$

Since it is not, it is not a toeplitz matrix.


Information i can take from coeficient intensity map in DCT transform What useful information i can take from the coefficient intensity maping after performing DCT transform?

The DCT is used in many places for its "energy compaction" properties. This is mostly useful in compression. For example, if I have a signal and I take its wavelet transform, then I can "compress" the signal by throwing away certain coefficients. These coefficients are not known before-hand. I have to look at the wavelet coefficients and then throw away all coefficients below a certain threshold to achieve some amount of compression with an acceptable loss in reconstruction accuracy.

For the DCT however, most of the energy of a signal is collected in its "low" frequency coefficients. What that means is that if you only keep the first lets say N coefficients (you don't need to store the locations of which coefficients you kept), you can get decent compression. In 2D, a "zonal" or "zig-zag" pattern is used to choose the coefficients to keep.

On a more fundamental level, the DCT is basically a linear combination of your signal as a bunch of discrete cosine waves. If you think of each cosine wave as an elementary building block, then the DCT tells you how much contribution to your image comes from a discrete cosine of a particular frequency.


Inverse of a summing function I am an experimental research engineer working on a Laser experiment. The details do not matter but due to practical limitations my data is influenced by the following transformation during data acquisition, and I was wondering if there is a mathematical method or an algorithm to "undo" its effect on the captured signal, that is to get the inverse.

I want to get the y in terms of s.
Thanks in advance

Looks like a convolution. S is the convolution of y with another signal of w+1 consecutive ones. So, S is the result of a moving sum filter. But the moving sum filter completely removes certain frequencies in the signal. For example, if you use w=3 and have the following signal y:

 y: ... 0  1  0 -1  0  1  0 -1  0  1  0  -1 ...

Then S[i] for any i (the moving sum over 4 neighboring samples) will always be zero, no matter how you shift or scale y. And that makes the problem of recovering y from s something between impossible and hard.

In case your data is finite, it is possible to write your equations as a matrix-vector product and solve it with the usual algorithms. You can ask tools like Matlab or GNU Octave about how noisy the inverse would be. This is something you should check for since the problem is ill-conditioned. Your equations, assuming y and s are finite can be written like this:

$$ A \cdot y = s $$

Here, A would contain mostly zeros and w+1 consecutive ones per row, at least for first rows. At the end, you have less ones because y is supposed to be finite. Then, asking for the condition number is a matter of typing

cond(A)

in Matlab or Octave. Example for 6 samples and w=2

octave:1> A = toeplitz([1 0 0 0 0 0],[1 1 1 0 0 0])
A =

   1   1   1   0   0   0
   0   1   1   1   0   0
   0   0   1   1   1   0
   0   0   0   1   1   1
   0   0   0   0   1   1
   0   0   0   0   0   1

octave:2> cond(A)
ans = 7.9053

Let's actually test it:

octave:3> y = [1 9 3 2 4 3]';
octave:4> s = A*y
s =

   13
   14
    9
    9
    7
    3

And now the inversion:

octave:5> y2 = A\s
y2 =

   1
   9
   3
   2
   4
   3

This inversion is actually easy to implement because the matrix exhibits a triangular structure. You can recover y sample by sample starting at the end and moving to the front.

This example was free of noise. Keep in mind that the condition number tells you how the relative error in s will be amplified (worst case) if you recover y from s. So, for example, if s contains a relative error of 3%, the recovered y may have a relative error as large as 3% * 7.9053 = 23.716%. Unfortunately, the condition number increases with the number of samples in this case. So, if you have to recover not six but 1000 samples and/or your w parameter is large, you're out of luck because the inversion will amplify the noise in s a lot, too.


Counting vechicles in an image I have been trying to implement an algorithm for successfully counting cars in an image. I have tried implementing A Method for Vehicle Count in the Presence of
Multiple-Vehicle Occlusions in Traffic Images

It estimates the background from a set of various images. I have looked into various other techniques for this purpose and all of these in one way or the other use either background estimation from a set of images or require a video. I have as an input, traffic images where the background (probably road in most of the papers) is hardly visible. Moreover, the images are from different areas so they don't have common backgrounds as well. How should I proceed in that case ?

I am thinking that if somehow I can match the structure of vehicles (car) then probably they can be matched. But I don't know if this is feasible and if it is how to proceed since the image contains multiple occluded vehicles as well.

Any hints or even research papers are also welcome.

A sample image is as follows:

As mentioned the object counting problem is very challenging. A good account of some common approaches is given in http://www.robots.ox.ac.uk/~vilem/NIPS2010.pdf.

Creating a SIFT feature database on training images would appear to be the natural path to go down. That in combination with some image segmentation may be the way to go.

Another path might be to look at HOG http://chrisjmccormick.wordpress.com/2013/05/09/hog-person-detector-tutorial/, which is similarly a feature detection algorithm, could be adapted for cars.


Cross correlation, when two signals are "uncorrelated"? Definition of cross correlation is

$$
\int_{\mathbb{R}}s_1^*(\tau)s_2(t+\tau)d\tau
$$

I wonder using such definition, when two signals are not correlated? My guess is the following...
defining
$$
s_{2,\tau}(t) = s_2(t+\tau) = \left(T_{-\tau}s_2\right)(t)
$$

We have

$$
\langle s_1,T_{-\tau}s_2 \rangle = 0
$$

for all $\tau$, namely $s_1$ is orthogonal to all the signals $s_{2,\tau}$. Is this guess correct?

This definition is taken from Wikipedia.

As assumption you can assume $s_1,s_2 \in \mathcal{L}^2(\mathbb{R})$, the definition should be then well posed by the Cauchy-Schwarz inequality.

Cross correlation is a measure of similarity between two signals, where one signal is allowed to be time-shifted. In this sense, the correlation is not a single number, but a function of the time shift. We say, "these two signal have a certain correlation $R(\Delta)$ for a time shift $\Delta$".

Intuitively, two signals that tend to have the same sign (both positive or both negative) for a given time shift $t$ are similar, and will have large correlation (positive or negative). Signals that are uncorrelated are just as likely to have opposing signs, and then the integral will be small.

As an example, consider two signals that are rectangular pulses of duration $T=1$, but one starting at $t=0$ and the other at $t=10$. Their product is zero (and thus their correlation for time shift 0 is also 0). However, intuitively we could say that they are similar, except for the delay. Thus, the correlation allows us to quantify that intuition: these two pulses are completely uncorrelated for some delays; however, when if the first is time-delayed by 10 seconds, then they become highly correlated.

Note that $R(\Delta)=0$ is a sufficient and necessary condition for the two signals to be "uncorrelated" (meaning they have zero correlation) with time shift $\Delta$.


Isolate High Variance vs Low Variance Sections of Signals I have the following one-dimensional signals:

The goal is to identify the sections with smaller variances compared to the rest (For Python, that is to get their index value in the array). In the first signal, there are two such sections (one at the left end and one at the middle). In the second signal, there is one such section at the middle.

What I've tried so far

I've tried computing the local stdev in each signals with window size 25 and produced the following plots:

From here on, I can use threshold to select sections with small stdev. However, one can see that both have different y-axis, and need different threshold values to isolate the low variance sections. Furthermore, applying a rolling window is notoriously computational expensive (I have to apply this on ~5,000,000 similar one-dimensional signals). So the question is, what is a quick method to identify sections that can be applied to both signals?

EDIT

I've put two sample signals at here and here.

What you did is the reasonable solution.
From here you can do 2 things to mitigate your issues:

Computation Efficiency
You can use online calculation of the Mean and the STD.
Remember when you move one sample to the right 24 samples are shared.
Hence if you take advantage of that you'll get a performance boost.
Threshold
Use adaptive threshold per signal if your classification of the STD is relative. Namely per signal you want to find zones with high variance relative to other zones in the same signal (And not an absolute measure).
You could use statistical properties of the STD of the signal like Mean / Median or even use adaptive Threhshold method like Otsu's Method.


How to reliably detect the state (opened or closed) of sliding doors from a statically-positioned video stream I have a mounted (static) camera that has a view of a set of sliding doors.  Since doors are rectangular in nature, and they don't have many degrees of freedom, this seems like a problem that should be relatively easy to solve reliably.  Here are some sample views (the door outlines in red I drew for illustration purposes, they are not detected by an algorithm):

First, here is a view of the sliding doors being closed:

Next, here is a view of the sliding doors being opened:

Finally, here is a view of the sliding doors opened but with occlusions:

My original thought was something around a hough transform, but I feel like there is likely a better solution to be had, given the knowledge that the doors can really only slide along one axis.

Reiterating my question, how can one reliably detect whether sliding doors are opened or closed from a video stream like the one in my images?

To better deal with occlusions, my idea would be to separate this problem into detecting if:

the 1st door is in position fully opened (1)
the 1st door is in position fully closed (2)
the 2nd door is in position fully opened (3)
the 2nd door is in position fully closed (4)

To tackle either of these problem, I would apply the following algorithm let's say with (1) :

1.1 Take a picture from the camera as a background image.

1.2 Crop it to keep only the part where the 1st door should be when it's fully opened, this is easy as the camera is static.

1.3 Do the same with the current image feed from your camera.

1.4 Calculate the absolute difference of those two croped images.

1.5 Calculate the average value of that absdiff  (optionnaly divided  by the surface of the croped part to normalize it).

1.6 Decide from that value if the door is fully opened or not ( low value mean it's probably opened, low value mean it's probably not), using a threshold value for example.

For steps 1.1-1.4, especially for updating the background  image, you could use a different background subtraction algorithm.

Once you can correctly decide, without much occlusions, if the statements (1-4) are true or not, you could, for dealing with occlusions, this works best assuming that both doors are always silmutaneously in the same open/close state:

5.1 keep the (normalized) absolute difference average of each (1-4) part

5.2 use some algorithm to decide if the doors are both open or closed

I'll illustrate step 5.2 with your last image as example, lets say 1st door is the one in the bottom of the picture (the only one we see well):

We have a high absdiff average value for both "Door 2 fully closed" and "Door 2 fully opened" as the man is in front of both places. So both  (3) and (4) are said as false(or true depends how you implement it), which is unlikely.
We have a high absdiff average value for (2) but a low one for (1).

-> So a possible algorithm for 5.2 could be using the difference between the absdiff average value from (1) and (2), and (3) and (4) as a threshold.

I think that, plus some noise removal techniques, will do the job.

Do not hesitate to ask me in comment if my solution doesn't appear clear to you.

Edit: didn't saw the date of the post, I don't know why it showed up so late in my feed, I let this answer in case it'll be useful to someone else.


How to set STFT-parameters to visualize woodpecker / Sawtooth Signal?  I always wanted to work with FFTs and Spectrograms to characterise sounds and their frequencies using python.

Lately I recorded a woodpecker at his work.

Woodpecker Wav File

To know where to look at I compared the sound with a generated sawtooth from Onlinegenerator.com. I think 17 Hz fits well.

17 Hz Sawtooth Wav File

Checking the code with the generated file worked good.

But sadly I cannot find the frequency in the DFT/FFT or spectrogram generated with the code below. No matter which settings I tried, I cannot visualise the frequencies below 20 Hz properly. No major signals cuts on or off. Only a signal at 23 Hz, which seems to be continuous.

from scipy import signal
import scipy.io.wavfile as wav
import matplotlib.pyplot as plt
import numpy as np

# Read WAV File
the_file = 'IMG_0710_short.wav'
samplerate, samples = wav.read(the_file)

#samples_left = samples[:,0]
samples = samples[:,1] # reduce to right channel

timevec = range(len(samples)) # Time vector for plot
timevec = [x / samplerate for x in timevec]
t_max =timevec[-1]

# Lineplot of Signal
dpi = 200
plt.rcParams['figure.dpi']= dpi

plt.plot(timevec,samples)
plt.title(the_file)
plt.ylabel('Amplitude [-]')
plt.xlabel('Time [sec]')

# STFT Settings
nperseg  = 0.5 * samplerate # Window Size 1
noverlap = nperseg*0.95 # Overlap
nfft     = 1 * samplerate # STFT 2
window   = 'hann' # Window Type

timespan = [0, 2] # Calculation Window
fromm = int(len(samples)/t_max*timespan[0])
too   = int(len(samples)/t_max*timespan[1])

f, t, Zxx = signal.stft(samples[fromm:too], samplerate, nperseg=nperseg, window=window, noverlap=noverlap, nfft=nfft)
t = t + timespan[0]
cmap=plt.cm.nipy_spectral
vmin = 10
vmax = 18
fig = plt.figure(figsize=(7, 5))
pcm = plt.pcolormesh(t, f, np.log(np.abs(Zxx)), cmap=cmap, vmin=vmin, vmax=vmax)

plt.ylim(0,50)
fig.colorbar(pcm)
plt.title(the_file)
plt.ylabel('Frequency [Hz]')
plt.xlabel('Time [sec]')

# FFT
n = len(samples) # length of the signal
k = np.arange(n)
T = n/samplerate
frq = k/T # two sides frequency range
frq = frq[range(int((n/2)))] # one side frequency range

Y = np.fft.fft(samples, norm='ortho')#/n # fft computing and normalization
Y = Y[range(int(n/2))]

fig, ax = plt.subplots(2, 1)
ax[0].plot(timevec,samples) # plotting the signal
ax[0].set_xlabel('Time')
ax[0].set_ylabel('Amplitude')
ax[1].plot(frq,abs(Y),'r') # plotting the spectrum
ax[1].set_xlabel('Frequency (Hz)')
ax[1].set_ylabel('|Y(freq)|')
ax[1].set_xlim(0,50)

How do I need to set the stft parameters in order to visualise the woodpecker spectrum correctly?
Or is there a fundamental problem with stft here?

Your first problem is that you have your axes mixed up.

Your "17 Hz" picture shows this quite clearly.  The red bars are the impulses from the signal.  The bars should indicate an impulse (wide range of frequencies) but your axes indicate they are across time.

This is what it should look like:

Notice that the bars are along the frequency axis rather than time as in yours.

This is your woodpecker recording:

It contains bars, too.  But, they are not as wide band and are much closer (in intensity) to the background noise.

The bars repeat at just above 20Hz - you can see that the frequency is around 23Hz in your last plot (the red line plot below the blue oscilloscope style plot.)

Your overlap of 0.95 amounts to a time resolution of about  just about 40Hz.  You need a higher overlap of maybe 0.99 to separate the bars.

You could downsample the recording down to like 11025 Hz sampling rate.  That would get you better frequency resolution and make calculating the spectrogram faster.  You can do that because there's not really much of interest in your signal above about 4500Hz.


$\lim \limits_{a \to 0} \frac{1}{a}[u(\frac{t}{a}+\frac{1}{2})-u(\frac{t}{a}-\frac{1}{2})]= \delta(t) $ It would be very helpful if someone could provide the intuitive reasoning and mathematical proof for:

$$F_a(t):=\frac{1}{a}\left[u\left(\frac{t}{a}+\frac{1}{2}\right)-u\left(\frac{t}{a}-\frac{1}{2}\right)\right]$$

$$  \lim_{a \to 0} F_a(t) = \delta(t) $$

where $u(t)$ is the unit step function and $\delta(t)$ is the unit impulse function.

The graph of $F_a(t)$ is a rectangle with base length $a$ and height $a^{-1}$. The area of the rectangle is therefore equal to $a \cdot a^{-1}=1$. The Dirac delta distribution is often defined as a certain (see below) limit of a sequence of suitable functions of integral 1 in the limit of vanishing support. This definition is already very intuitive, so that your intuition criterium should be met.

A more rigorous look uses a different definition of the delta distribution through the equation

$$\int_{-\epsilon}^{+\epsilon} f(x) \delta(x) dx = f(0)$$

for any open neighbourhood of radius $\epsilon$ around $0$ on which $f$ is smooth, as well as the condition $\delta(x\neq 0)=0$.

If you plug in your $F_a(x)$ in place of $\delta(x)$ we get

$$\int_{-\epsilon}^{+\epsilon} f(x) F_a(x) dx$$

Now the limit of $F_a$ for $a\to0$ is not to be understood of a limit in the usual function sequence sense. That limit does not exist in terms of a function. Rather, the limit is meant to be taken outside an integral under which $F_a$ appears. So taking the limit happens like

$$\lim_{a\to0}\int_{-\epsilon}^{+\epsilon} f(x) F_a(x) dx$$

In order to understand the limit, we can start at any finite but small enough $a$. If we choose $a<2\epsilon$ and use that $F_a$ and with it the function under the integral vanishes outside of $[-a/2,a/2]$ and offer no contribution to the integral, the boundaries can be adjusted accordingly.

$$\lim_{a\to0}\int_{-a/2}^{+a/2} f(x) F_a(x) dx$$

Within these boundaries $F_a$ is constant and equals $1/a$. So the integral simplifies to

$$\lim_{a\to0} \frac{1}{a}\int_{-a/2}^{+a/2} f(x) dx$$

The requirements for $f$ are that it be smooth on the open interval of our integral. So we can expand $f(x)=f(0)+x \cdot f'(0) + \mathcal{O}(x^2)$ which together with the linearity of the integral results in

$$\lim_{a\to0} \frac{1}{a}\left( f(0) \int_{-a/2}^{+a/2} dx + f'(0) \int_{-a/2}^{+a/2} x dx + \int_{-a/2}^{+a/2} \mathcal{O}(x^2) dx\right)$$

The first integral contributes $a$, the second integral vanishes because of symmetry of the boundaries and the odd function under the integral, and the third integral contributes $\mathcal{O}(a^3)$.

This gives

$$\lim_{a\to0} \frac{1}{a}\left( f(0) a + \mathcal{O}(a^3)\right) = \lim_{a\to0} \left( f(0)  + \mathcal{O}(a^2)\right) = f(0) $$

This result tells us that your $F_a$ works as a $\delta$ distribution if the limit is taken outside the integral. Or in other words, the "limit" of the function sequence is a delta distribution.

Edit: I just noticed I did not go into the second condition of $\delta(x \neq 0) = 0$. To show this for $F_a(x)$, just observe that for any $\hat{x}\neq0$ you can find an $\hat{a}$ so that $F_a(\hat{x})=0$ for all $a<\hat{a}$.

Edit 2: Why is the condition $\delta(x \neq 0) = 0$ necessary at all? This one is a bit subtle. The integral condition should uniquely define the distribution since it must hold for all test functions $f$ from the space of locally smooth functions. The problem is that the test functions are not enough to rule out distributions that are identical to $\delta$ only almost everywhere. For example, if we had a distribution $\bar{\delta}$ that coincides with $\delta$ with exception of $\bar{\delta}(1)=1$, then the integral condition would still be met because the set that supports the deviation is of measure zero. The requirement of vanishing $\delta$ for non-zero arguments excludes these cases. It is however often omitted because it does not present a problem unless you let the delta distribution act on other distributions.


Prove using superposition that an LTI system is linear Given an LTI system of the form:

$$ \frac{d}{dt}{x}(t) = {A} {x}(t) \ + \ {B} {u}(t) $$

where ${x}(t)$ is the state variable, ${u}(t)$ an input and ${A}$ and ${B}$ constants, how does one show, by superposition, that the system is linear? I've searched and searched and tried to do it myself and I am unable to do so. Seems more difficult that it seems. Most books seem to gloss over the fact.

I'm going to assume the system has the initial condition $x(0) = 0$, as otherwise the system isn't linear.

Suppose the input $u_1(t)$ yields the output $x_1(t)$ and the input $u_2(t)$ yields the output $x_2(t)$.

We need to prove that the input $u(t) = \alpha u_1(t)+\beta u_2(t)$ yields the output $x(t) = \alpha x_1(t) + \beta x_2(t)$ for any scalars $\alpha$ and $\beta$.

The given input and output pairs must satisfy the differential equation and the initial condition. Hence, $x_1'(t) = Ax_1(t)+Bu_1(t)$ and $x_2'(t) = Ax_2(t)+Bu_2(t)$ for all times $t$ and $x_1(0) = x_2(0) = 0$.

Now, add $\alpha$ times the first equation to $\beta$ times the second equation to get: $\alpha x_1'(t) + \beta x_2'(t) = \alpha(Ax_1(t)+Bu_1(t)) + \beta(Ax_2(t)+Bu_2(t))$

Then, rearrange stuff to get:

$\dfrac{d}{dt}[\alpha x_1(t) + \beta x_2(t)] = A(\alpha x_1(t)+\beta x_2(t)) + B(\alpha u_1(t)+\beta u_2(t))$

$x'(t) = Ax(t)+Bu(t)$

Also, we have $x(0) = \alpha x_1(0)+\beta x_2(0) = \alpha \cdot 0 + \beta \cdot 0 = 0$.

Thus, the signals $u(t) = \alpha u_1(t)+\beta u_2(t)$ and $x(t) = \alpha x_1(t) + \beta x_2(t)$ satisfy the differential equation and the initial condition. Hence, the input $u(t) = \alpha u_1(t)+\beta u_2(t)$ yields the output $x(t) = \alpha x_1(t) + \beta x_2(t)$ for any inputs $u_1(t)$ and $u_2(t)$ and any scalars $\alpha$ and $\beta$.

Therefore, the system is linear, as desired.


How to calculate Long Term Average Spectrum of an audio file? I'm working with audio files and I have to find their long term average spectrum. Although I had an option available in the software Praat but I don't understand the algorithm of finding the spectrum. I need to understand it so I know the actual working and also implement it through MATLAB!

Sounds like a conventional spectral analysis, just do an FFT over your entire  audio file and you obtain an averaged spectrum of the file. A Fourier transform can be viewed as a set of band-pass filters with integration.


Fourier descriptors: trying to classify objects Describing my background: I have around 33 items labeled. For example, 3 pictures of the contour of a basil plant, 4 pictures of the contour of earphones, 7 of a mug, etcetera.

I'm trying to calculate N Fourier descriptors of a random picture in order to compare them with the other items, making the average of the distances between classes and saying: "This is a mug!" or "this is a basil plant!".

In my experience, I take a random picture of a MUG:

it looks to my MUGS that I've classified (because it gives me a distance close to 0):

but if I take a picture of the MUG without the handle:

the closest averages of the norm of distances to this MUG are of BASIL plants or EARPHONES (instead MUG, which is boundary closer):

indeed I have been realising that for every random picture I use to get a low value for earphones and basil. It seems that I have to normalize somewhere... but... what should I do ??

Do anybody could give me a clue?

Thank you very much in advance.
Regards.

UPDATED
In response to @tbirdal , I have to say that that's what I do:

vector<double> vector;
if (dft_z.size() < 2)
    return vector; // handling errors
dft_z.at(0) = complex<double>(0,0); // First DC to 0 (invariance on translation)
double si1 = abs(dft_z.at(1)); // Module of the a(1)
for(int i = 0; i < dft_z.size(); i++) {
    dft_z.at(i) = complex<double>(dft_z.at(i).real() / si1, dft_z.at(i).imag() / si1); // dividing by the module
    vector.push_back(abs(dft_z.at(i))); // and putting into my vector the absolute value (to eliminate any phase)
}

HOW I'M GOING TO CONTINUE

I've seen that making the absolute value of every FD to remove its phase makes confusion on my shapes.
In addiction, I realised that only simple forms and contours are well matched. For example basil was too complicated, but the mug is not.

Thanks to @tbirdal who gives me an idea.

The normalization of the Fourier descriptors is performed as follows:

1) Set the DC component of the descriptors to $a(0)=0$.

2) Divide all the Fourier descriptors with the magnitude of the second one. i.e. :

$a(1)=r_1e^{j\phi_1}$

$a(k)=\frac{a(k)}{\lVert a(1) \rVert}$

3) Now, only the position of the starting point remains to be normalized. This is done by subtracting the phase of the second Fourier descriptor $\phi_1$ frin the phase of all Fourier descriptors and weighting by $k$; that is:

$a(k)=a(k)e^{-j \phi_1 k}$

Note that I did not describe rotational invariance here. In simple terms it is achieved by taking the magnitude of each Fourier coefficient.

Also, you have to chain code your contours and represent them in complex space to perform the Fourier transformation (the ordering matters if you don't perform step (2) ).

Fourier descriptors are not necessarily robust. They are good if you have some control of environment (such as lighting and appearance etc. ), but I don't expect them to work on cases such as a missing handle.


Finding optimal decision threshold for binary comm system I'm trying to solve the following exercise: Given a system that transmits bit $b$ with probability $p_b$ and $-b$ with $p_{-b}$ and the transmission is degraded by AWGN $N(0,\sigma^2)$. What is the optimal decision threshold $\mu$ to minimize the probability of error? Use Q functions to express the probability of error.
My approach:
$P_e = P_b P(-b|b) + P_{-b} P(b|-b)$, where $P(-b|b)$ means probability of detecting $-b$ given that $b$ was transmitted.
$P(-b|b)$ is the left tail of the distribution ($\Phi$) since I need to receive a value smaller than $\mu$ after sending $b$, but since the exercise asks in $Q$ -> $Q(-(\frac{\mu-b}{\sigma}))$
$P(b|-b)$ is already the right tail as I need to get something greater than $\mu$ after sending $-b$, so $Q(\frac{\mu+b}{\sigma})$.
To find the optimal $\mu$, I took the derivate of $P_e$ with respect to $\mu$ and set it equal to 0.
To facilitate notation, I will use $X = \frac{-(\mu-b)}{\sigma}$ and $Y = \frac{\mu+b}{\sigma}$.
Now I have $P_1 e^{-\frac{X^2}{2}} = -P_{-1} e^{-\frac{Y^2}{2}}$.
After some manipulation and taking $ln$ on both sides: $ln(\frac{p_1}{p_{-1}}) = \frac{\frac{Y^2}{2}}{\frac{-X^2}{2}} = \frac{2\mu b}{\sigma^2}$.
So, $\mu = ln(\frac{p_1}{p_{-1}}) \frac{\sigma^2}{2b}$ but the solutions I have don't have the $\sigma^2$ term, so I cannot find where I made a mistake.
As for the derivative of the $Q$ function, it should be $-\Phi'$, right? Which ends up being the pdf, so $e^{-\frac{z^2}{2}}$ (omitting the normalization factor as they will cancel out on this particular question)

Let P(0) ( P(1) ) is a probability of transmitting bit zero (one); P(e|0) ( P(e|1) ) is a probability of error when detecting bit zero (one).
The probability of erroneous detection is
$$
P(e) = P(e|0)P(0) + P(e|1)P(1) \tag {1}
$$

Let $V_0$ ( $V_1$ ) be a nominal signal voltage of bit zero ( one ) signal at the transmitter.
$$
P(e|0) = \int_T^{\infty}{{\frac {1} {\sigma \sqrt{2\pi}}}\exp\left(-(\nu-V_0)^2/2 {\sigma}^2\right) d\nu } \\
P(e|1) = \int_{-\infty}^T{{\frac {1} {\sigma \sqrt{2\pi}}}\exp\left(-(\nu-V_1)^2/2 {\sigma}^2\right) d\nu }
$$
where $T$ is a detection threshold (OP's $\mu$).

Differentiating $P(e)$ of eq.1 w.r.t. $T$, we arrive at
$$
-P(0){{\frac {1} {\sigma \sqrt{2\pi}}}\exp\left(-(T-V_0)^2/2 {\sigma}^2\right)} + P(1){{\frac {1} {\sigma \sqrt{2\pi}}}\exp\left(-(T-V_1)^2/2 {\sigma}^2\right)} \tag {2}
$$
To find an optimal threshold, we equate the expression in eq.2 to zero:
$$
P(0)\exp\left(-{\frac {(T-V_0)^2} {2 {\sigma}^2}}\right) = P(1)\exp\left(-{\frac {(T-V_1)^2} {2 {\sigma}^2}}\right) \\
T = {\frac {V_0+V_1} {2}} + {\sigma}^2{\ln{(\frac {P(1)} {P(0)})}}
$$
that is, exactly the OP's  result.
The OP's reference source may compute the optimal threshold for the system in which the probabilities of zero/one bit transmissions are equal. The variance ${\sigma} ^2$ is multipled by $\ln{(\frac {P(1)} {P(0)})}$, and it disappears from the expression for optimal threshold, when the zero/one probabilities in transmission are equal ($P(0)=P(1)=1/2$).


Sketching Phase Spectra Using Group Delay and Magnitude Spectra Informations I have only magnitude spectra and group delay information and I need to sketch phase spectra from this information. For example, group delay is given like this: $\tau_{g}(\omega) = c$ where c is a constant value.

I know the group delay is the derivative of minus phase function of the system. In this case, phase function can be written likse this: $\beta(\omega) = - c \omega + constant$, but how to determine the $constant$ here? Is it possible to find phase from these information?

If the magnitude spectrum is symmetric

$$M(\omega)=M(-\omega)\tag{1}$$

(as I assume), then your system is real-valued. The phase response of a real-valued system is asymmetric:

$$\phi(\omega)=-\phi(-\omega)\quad(\mod 2\pi)\tag{2}$$

This means that there can be two cases:

The phase goes through zero at $\omega=0$, i.e. the phase is given by $\phi(\omega)=-c\cdot\omega$, where $c$ is the constant group delay.
The phase is $\pm \pi$ at $\omega=0$, which means $\phi(\omega)=-c\cdot\omega\pm\pi$.
The phase jumps at $\omega=0$. This is only possible if the magnitude has a zero at $\omega=0$, i.e. $M(0)=0$. The phase jumps by $\pi$, which is simply a sign change of the (bipolar) amplitude function. In this case the phase is given by $$\phi(\omega)=\begin{cases}-\pi/2-c\cdot\omega,&\quad\omega>0\\\pi/2-c\cdot\omega,&\quad\omega<0\end{cases}$$


What is the theoretical upperbound of information that could be transmitted by a device per second? Suppose I would like to stream data from a device in a single second. What would be the theoretical upperbound of the amount of information (in bits, bytes, pixels...) that could be sent during this period?

What determines this upper bound?

There are many factors involved in understanding the theoretical limits to communication. What follows is just a brief introduction that only scratches the surface.

First, let's consider a simple scenario: there is no noise and no distortion of the signal being transmitted. We do allow for attenuation. Under these circumstances, you can transmit up to $R_p=2B$ pulses per second. $B$ is the bandwidth you have available. For example: if you have 10 kHz of bandwidth, you may transmit up to 20,000 pulses in one second. $R_p$ is called the Nyquist rate, after Harry Nyquist.

However, your question is about bits, not pulses. It turns out that you can transmit bits using pulses. Let's say you want to transmit the bit sequence $1001$. You could do this by transmitting pulses with amplitudes $+1, -1, -1, +1$. (The pulse shape does not matter, but they do need to "fit" in the bandwidth you have available.) If you do this, your bit rate $R_b$ is equal to your pulse rate, $R_b=R_p$.

You can transmit more than one bit per pulse, though. Consider allowing the pulses to have four amplitudes $+3,+1,-1,-3$. Then you could transmit the sequence $1001$ with only two pulses, of amplitudes, say, $+1,-3$. You're free to map bits to pulse amplitudes in any way you want. In this scenario, you have $R_b=2R_p$.

So, in theory, if there is no noise, then the pulse rate is limited to twice your bandwith, but the bit rate can be as large as you want, by allowing the pulses to take more and more different amplitudes.

Now consider a more realistic scenario, where noise is present. Now, the received signal is not just what you transmitted, but it includes noise too. Assume that you're transmitting with power equal to $S$ watts, and that the noise over your bandwidth $B$ has power $N$ watts. Now your bit rate is no longer as large as you want; it cannot be larger than

$$C=B\log_2(1+S/N).$$

$C$ is known as Shannon's channel capacity, after Claude Shannon. However, keep in mind that transmitting at, or near, capacity may involve a very complex coding system. In practice, a simple system will operate well below the channel capacity. Even very complex systems may not achieve the capacity.

What happens if you only have a simple system, and you attempt to communicate at a significant fraction of the capacity? In that case, you're likely to have errors in your communication. An error happens when the receiver mistakes a bit 1 for a bit 0 or vice versa. Depending on your application, a certain amount of errors can be tolerated. For example, in voice communications one error out of 1000 bits (a bit error rate (BER) of 0.001) is acceptable. Other systems are intolerant to errors: for example, you don't want your hard drive to corrupt your files, not even once a year.

In more complex situations the capacity may be much smaller and the systems that achieve high bit rates become much more complex. For example, in wireless communications the transmitted signal suffers not only from noise, but also from echos and Doppler shift.

In summary, as you see, there is no simple answer. Let me know if you need clarification, or more details.


MLE formulation -- confusion regarding the terms in the equation (Part1) If $v \sim CN(0,2\sigma^2_v)$ is a circularly complex Gaussian random variable which acts as the measurement noise in this model
$$y_n = x_n + v_n \tag{1} $$ where $x \sim CN(0,2\sigma^2)$, then is the following log-likelihood expression correct?

$$P_x(x_1,x_2,...,x_N) = \prod_{n=1}^N\frac{1}{2\pi \sigma^2_v} \exp \bigg(\frac{-{({y_n-x_n})}^H ({y_n-x_n})}{2\sigma^2_v} \bigg) \tag{2}$$

Confusion:

I am having a doubt whether there will be a 2 in the denominator in (3) from the term $2\sigma^2_v$ or not. Am I doing it correctly? Please help.

Hi: I don't have time to carefully look at everything right now and I'm unfamiliar with circular complex random variables. But, if it's similar to regular normal rv's and the two variables are independent , then their sum has mean zero and their variance is $ 2 \sigma^2 + 2 \sigma^2_{v}$. So, the likelihood  can be written having mean zero and that variance.  But it's not clear to me what you're estimating. There aren't any coefficients in the model ( implied one's I guess somewhere or maybe you are doing bayesian ? ) so what is actually observed ? $x_{n}$ or
 $v_{n}$. That has to be known if you want to write the likelihood correctly.
If $x$ is viewed as fixed and known and not a random variable, then what you wrote is missing a factor of 2 in the factor before the exp. The $2 \pi$ has to be multiplied by the variance which is $2 \sigma^2_{v}$. Also, it's not clear to me how X can be viewed as fixed since you said in the question that it was an RV also. If it is an RV and also not observed, then that has to be accounted for in the likelihood.  That's all I can tell you and even what I'm saying may not be correct if complex means something totally different. Good luck.


Wavelet time extraction I am trying to understand one simple thing in wavelets and I searched a lot, but did not get any idea.

I need to perform Multilevel 1-D wavelet decomposition on my signal. I have arrays

signal[] and time[]

( time[x] corresponds to the time at which signal[x] occurs ) . After I perform

wavedec(signal,level=1,'wname'))

I need to plot cA1 vs the time signal. But how do I know which time points should I select. If i change the 'wname' , I get different length cA1.

For example, for a signal of length 9285, when I do

[cA1,cD1]=pywt.wavedec(signal,db1, level=1,mode=pywt.MODES.zpd)
length of cA1=4643

[cA1,cD1]=pywt.wavedec(signal,db2, level=1,mode=pywt.MODES.zpd)
length of cA1 = 4644

How can we determine the variation of cA1 and how can we select the corresponding time. I use pywavelet

The length of the arrays you get back depends on the mode you use to do the decomposition and the wavelet length, which can vary like you see between db1 and db2:

len(cA) == len(cD) == floor((len(data) + wavelet.dec_len - 1) / 2)

If you switch to periodization mode mode='per', you get a different calculation that may be more consistent for you:

len(cA) == len(cD) == ceil(len(data) / 2)

Regardless, decomposition is generally dyadic, so your dyad width is len(signal)/len(cAn).

Once you go muiltilevel (use level='2', etc.), you'll have small variations in length as well. To account for this, use IDWT to reconstruct.


Linear classification and overlapping problem reduction I am study linear classification(Y=aX+b) But it is still weakness that is overlapping problem. Could you suggest to me how to reduce overlapping problem in linear classification or how to find these parameters that helpful for overlapping reduction. In my question, I don't considier SVM method because it is more difficult to implement. Thank you so much

It's very difficult to tell what you're asking exactly, but one thing to keep in mind is that not al data is linearly separable. You may have to obtain more features from data or try non-linear separation techniques. SVM is one of them, but since you're not considering it, you can try polynomial separation or logistic classifiers.


Pixel Wise Write into a Video File in MATLAB A video file (in .ptw  format) is read in MATLAB in the form of a 3 dimensional matrix (a structured format). My aim is to process the video in time domain on each coordinate pixel individually and eventually write the processed time domain output in a video file. I am looking for a technique for the same.

Please see - The processing are of type cross-correlation, filter in fourier domain, that requires the data to be transformed to 1-d matrix.

I assume your matrix is 3D where the first 2 dimensions are Width and Height and the third is time (Gray Scale Video).

If you can write your processing as a Filter you can use MATLAB's filter() function.
The function has the dim property hence you can set it to 3 and it will apply on the filter over the third dimension.
The limitation is applying the same filter to whole data.
If you need to have different filter for each pixel then you should iterate over the filter bank.

Another solution is to use the functions permute() and reshape() to rearrange the data to have the Time Axis in the dimension which fits your needs.


How to design very narrow band, very sharp low pass filters - Only DC needed? I am using a DSP processor to sample signals at 108kHz and I only want to get the DC part out of it. Such a high sampling rate is used because I want to use oversampling to reduce quantization noises.

Two filters are to be designed:

One has 10Hz stop-band frequency; new data are read every 0.1 second;
One has 0.1Hz stop-band frequency; new data are read every 10 second;

What I did

I used several stages of decimation and downsample the signals to 300Hz. Then I use FIRs to filter the signal:

For the 10Hz bandwidth: 128 stage filters are used
For the 0.1Hz bandwidth: 768 stage filters are used.

My question

Performance is the most important to me: The sharper low-pass filters are, the better.

Could I ask for suggestions in implementing such a filter?
Is it possible to design a 0.01Hz low pass filter? If so, the new data can be read every 100 seconds?
I heard that using IIRs could be a solution, is it?

Thank you very much!

Is averaging the best DC filter?

I tried to compare with matlab:

Setting 1:
%Use Matlab fir1 function

    size = 1080000;
    h = fir1(size, 0.00000000001);
    fvtool(h)

Setting 2:
%Build an averaging FIR

    size = 1080000;
    h = ones(1,size);
    h = h *1/size;
    fvtool(h)

The theory should be both two settings have similar performance. However, it seems that fir1 has narrower stopband than the simple average filter. It seems that the average DC filter will include more noises.

If all you want to do is get an estimate of the DC component of your data then you can simply obtain the average value of your signal over a unit of time. This can be easily shown by looking at the expression for calculating the Fourier coefficients for 'f=0' (DC component).

The easiest way to do this is via a moving average filter (http://en.wikipedia.org/wiki/Moving_average) which would work over a buffer.

In terms of an FIR filter, this looks like 'h = ones(1,Nbuff)./Nbuff' (so, 'h=[0.5,0.5]' for a 2 sample buffer). 'Nbuff' stands for the order of the filter. This can be applied directly on the buffer you obtain from your CODEC through a convolution function as provided by your processor's command set.

You could also use an 'online algorithm' (an IIR filter essentially) which maintains an updated estimate of the mean of the signal on a per-sample basis. For this purpose please see: http://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Online_algorithm You could also design a "typical" IIR filter that can be as strict as you like, provided that your processor and sampling period can accommodate the computations required.

(Of course you can adapt between working in terms of frames and working in terms of individual samples)

Just a note, I have a feeling that you are struggling with the concept of Group Delay. The cut off frequency of your filter does not exactly mean that "new data can be read every 100 seconds". New data will be fed to your filter constantly, at the sampling rate of the system, but the phase delay caused on each harmonic component of your signal between the input and the output of the system might vary. How this varies, depends on the type of the system (whether it is FIR / IIR) and its form. Here is a starting point: http://en.wikipedia.org/wiki/Group_delay_and_phase_delay

Hope this helps.


GPS Code Phase usage During Cross-correlation process:

xcorr=ifft(fft(a)fft(b)∗)

we have frequency offset and code phase offset.

How is it used later in the process of finding a user location?

To make it simple, finding user location is done by using four pseudo-range measurements (approximative distance from user/receiver to satellite) to resolve a four unknown equation. As information, these unknown are receiver's (X,Y,Z) coordinates and it's clock offset.

These pseudo-range measurements are done by tracking stage of the receiver. Acquisition results (code phase and Doppler shift) are used to initialize delay and phase tracking loops. I let you the following simple explanation of the receiver behaviour : https://gssc.esa.int/navipedia/index.php/Generic_Receiver_Description

Please let me know if my explanation are not sufficient,

Math


what is the effect of variation of duty cycle of sampling frequency on a reconstructed signal I have been asked  to study an experiment which says "what is the effect of variation of duty cycle of sampling frequency on  a reconstructed signal"

I am not an electronics/electrical engineer I am CS guy.This experiment probably is a part of Analog to Digital Converion techniques, by now I have been able to understand Nyquist theorem and reconstruction is based on this.But for above experiment I am not clear with what to do or study so far what ever I have googled has not given me fruitful results let me know if some one has an answer to above problem.Following is the link which I have studied and could follow http://mhakhan.tripod.com/07sampling.pdf  a bit for the above question,what is exactly the thing which I should search?

First establish the following convention:

$$DutyCycle = \frac{SampleClockOnTime}{SampleClockPeriod}$$

In the following discussion,

$$SampleClockOnTime = T_p$$
$$SampleClockPeriod = T_s$$

You can model the effects of $T_p$ on your recovered signal as a pulse amplitude modulation signal recovery problem.  The sample clock signal is a pulse with finite duty cycle.  The reconstruction signal starts by modulating the sample clock signal with the data values of your sampled signal.  Literally you can think of multiplying your samples values by a unit amplitude pulse (your sample clock) to create a physical time sequence.

In the frequency domain, you end up with a reconstructed signal that has altered frequency characteristics (roll off at the high frequencies) that are a direct consequence of the sample clock pulse width (aka duty cycle).

To complete the recovery of your signal, you will pass the modulated pulse signal through a low pass filter to obtain only the baseband components of the signal.  The resulting signal with have frequency magnitude characteristis of the form:

$$|F'(f)| = |F(f)||\frac{T_psin(pifT_p)}{T_s*f*T_p}$$

Where:
F'(f) is the recovered signal and
F(f) is the spectrum of the origianl sampled signal

The effects of the sample clock pulse width is represented by:
$$\frac{T_psin(pifT_p)}{T_s*f*T_p}$$

If you plot this, you will see the amplitude rolls off as frequency approaches the sample clock frequency and the roll off is becomes more pronounced as $T_p$ increases.

You can look for more information on this by searching on the topic of "Sample and Hold effect".  Most of the work I've seen on this topic is in the context of switched capacitor filters (which is how I was exposed to it).


Constant used to determine the minimum samples in Compressed Sensing A key result in compressive sensing states that, with high probability on the random
draw of an $\ m × N$ Gaussian or Bernoulli matrix $\ A$, all $\ s$-sparse vectors $\ x$ can be reconstructed from $\ y = Ax$ using a variety of algorithms provided

$\ m ≥ Cs\ln(N/s)$,

where $\ C > 0$ is a universal constant (independent of $\ s, m,$  and $\ N$). This bound is in fact optimal.

In practical terms, how does one interpret/find $\ C$?

source: http://www.cis.pku.edu.cn/faculty/vision/zlin/A%20Mathematical%20Introduction%20to%20Compressive%20Sensing.pdf

The constant C is usually left unspecified as it can be very hard to calculate.  It is independent of the signal dimension, and results in Compressed Sensing are generally only concerned with big-O results, so it falls out when we say something like $m = \mathcal{O}(s \ln{(N/s)})$.  The important part in this statement is the strong dependence on the sparsity level, and that the signal dimension $N$ only logarithmically affects the number of measurements.

To give an example of where such a constant would come from:  Let $V_d$ be the volume of the unit $\ell_2$ ball in $\mathbb{R}^d$.  Some simple cases are that $V_2 = \pi r^2$ and $V_3 = \frac{4}{3} \pi r^3$.  Note that both of these results are of the form $V_d = C r^d$, and this pattern generalizes for finite $d$.  In CS, it is $r^d$ that we are interested in, as it captures the dependence on the signal dimension.

So $C$ will vary depending on the algorithm as different inequalities may be used in the proofs.  I don't know of any publications that specifically calculate C.


What do you call it when you add zeros between the elements of a vector? What do you call it when you add zeros between the elements of a vector?

let's say you have x =[1 2 3 4]

and x_2 = [0 1 0 2 0 3 0 4]

what do you call this process?

In the context of filter banks, simply upsampling (without filtering, or with the trivial all-pass filter, as Wikipedia relates it to interpolation), denoted by an up-arrow, as illustrated below fro a 2-fold or order-2 upsampling:

This is described in Filter banks: Decimation, Interpolation, and Modulation, and the  interpolation filter associated with upsampling is called synthesis filter.

As noted by @Jason R, one also finds "expansion". Expansion is also found as a streching for continuous signals, while compression is the effect in the dual Fourier domain. Sometimes, one can find "zero-insertion" or "zero-stuffing".

Actually in image processing some use "upsampling" as a proxy for pixel increase with filtering.


Finding the best principle component The title might be unclear, but the problem is this. I have a signal sampled 1500 times with a rate of 60/s, and a sensor array 512 units large. There is a lot of noise, echo and other frequencies being picked up, but I am interested in only one. First I do a spike removal, then a bandpassfilter (butterworth) around the frequency range I suspect the signal is hidden. I then do a PCA to find wich of the 512 sensors picks up any systematic variation. Then I search for the best fit sinuswave in each of the top sensors (say 5 out of 512).

So main question, how to determine which of the PCA's is the one picking up the true signal without me knowing excactly what the frequency of the signal is? Second question, does the above steps seem reasonable? I am no expert in this, but experiments seems to indicate that with high STN objects of measurements (say a pendulum) it is clearly visible which PCA fits best (amplitude, hz, residuals of fit), but with low STN (say heartbeat), it is not so clear.

Sorry for lengthy text. Thanks alot for any answers!

edit: Running spectral analysis (using different Methods) gives different results depending on wich PC i am analysing. Might there be a spectral Method that does multivariate samples?

edit 2: So while waiting for answers I have chosen to filter out principle components based on the correlation coefficient from a sinus wave fit (r>0.5) and % of data explained (p>0.15). So in other words, the signal i am looking for has to be present in at least 15 out of the 512 sensors (if I understand PCA correctly), and the best fit sinewave has to explain at least 25% of the variation. This works well with test setups in a noisy environment. Question: Is my approach sensible? Should or shouldnt I forego FFT and spectral estimation? Besides visualization, does frequency estimation based on the FFT gain anything over a iterative sinewave fitting, given that I know the signal is sinusoidal?

Thanks alot

By doing PCA, the principal components you will get will not correspond to a single recording, but rather to a mix between them. PCA is a feature extraction method, whereas what you are looking for, seems to me, as a feature selection problem.

Also, if you have so many simultaneous recordings, all affected by the same sources of noise, why not perform active noise cancelation?


Max useful filter tap count for given fixed point bitwidth I have a 10 bit ADC sample stream in which I would like to apply a digital band pass filter. Is there a theoretical filter tap count that would define the threshold of "usefulness"?

I guess my thought process is: There are only 10 bits, so at some point (via increasing the tap count) the attenuation will be so great that even a maximum amplitude frequency component will get attenuated down to a fraction of a bit, and rounded down to 0. At that point, the effects of additional taps could not be seen.

Since you can widen your bit width with each operation, no, there's no such maximum useful width.

Take GPS receivers as an example: many of them citation needed use 2 bit ADCs and yet are useful, because the signal, already hidden in noise before it reaches the ADC, only gets "visible" by massive processing gain.

A typical step with massive processing gain would be a very sharp digital filter. Those can get very long.

There's a lot of experienced audio and video folks on here. They can tell you much more about the astronomical lengths high-quality filters can get.


How to determine the frequency of a waveform? I am getting data from a sensor which gives 64 samples per second. If I take 64 as the sampling frequency, what would be the frequency of the signal?
I have applied a Fourier transform on the signal. I got a double sided spectrum with X-axis ranges from 0-64, if I take a positive side of the spectrum I want to determine the frequency of those points.

Image1: plot of samples for 1 sec
Image2: FFT of 1 sec signal

Question 1: the sampling frequency would be $64$ Hertz
Question 2: the frequency indices range from $0$ to $63$

In the following code, I tried to mimic your data. At each run, a different observation is produced, based on a sine with a different amplitude distortion and offset.

Using the FFTR.m code that displays one half of the spectrum for real signals, and removing the average, you can get the following figure:

% Simulated data to mimic your signal
timeSamp = 1/64;
nSample = 64;
xAxis = linspace(0,1,nSample)';
dataClean = sin(2*pi*xAxis*6);
distAmpl = 1/10*( medfilt1(rand(nSample,1),3)+0.15);
dataDist = data.*distAmpl + 0.01*(rand+10);

% Half-axis FFT for a real signal
 [fftR,fftAxe] = FFTR(dataDist-mean(dataDist),timeSamp);

% Display
subplot(2,1,1)
plot(xAxis,dataDist,'.-');grid on;axis tight
xlabel('Time');ylabel('Amplitude');
subplot(2,1,2)
plot(fftAxe,fftR,'.-');
xlabel('Frequency');ylabel('Amplitude');
axis tight;;grid on


Is the Kalman Filter a Best Linear Unbiased Estimator (BLUE) for Heteroscedastic Noise? According to the Gauss-Markov Theorem, a ordinary least squares estimator is BLUE if the noise entering a system is uncorrelated with zero mean and is homoscedastic (has a constant finite variance). I am aware that a Kalman Filter applied to a system with additive noise of known mean and variance but non-gaussian distribution is  BLUE. Does this imply that the noise has to be homoscedastic ? Or does the KF have a trick up its sleeve ?

Kalman filter is the best linear estimator regardless of stationarity or Gaussianity. Also in the Gaussian case it does not require stationarity (unlike Wiener filter). In the linear Gaussian case Kalman filter is also a MMSE estimator or the conditional mean.


Map the Amplitude Envelope from [dbmin,dbmax] to [-inf,dbmax] I would like to map the envelope of a sound (which has max dB level = dbmax, and min dB level = dbmin)   from [dbmin,dbmax] to [-inf,dbmax].
Here is an example :

It's like some sort of dynamic compression (or rather dynamic range 'expansion' to be precise), but here the mapping is linear from [dbmin,dbmax] to [-inf,dbmax].

What's the simplest way to do this with an  x[n] array? (my input is a soundfile, not an envelope. Is there a simple method for working on envelope (envelope follower ?))

Thank you.

PS : not realtime processing. Thus I can pre-compute dbmax, dbmin.

May be you can try the following:

 dbnew = dbmax * (1 + log((db-dbmin)/(dbmax-dbmin)))

dbnew = -inf when db = dbmin and dbnew = dbmax when db = dbmax. Please let me know if it works.


Error analysis for uniform quantization of uniform input I was told that

  error energy $$\delta^2_e=\Delta^2/12,\quad \Delta=(B-A)/2^R$$
  signal energy $$\delta^2_x = (B-A)^2/12$$
  Here signal bounded:$A\le x[n]\le B.$ $[A,B]$ is split into $2^R$ subintervals.

But i don't understand how the signal energy was derived. More precisely, i don't understand why signal energy equals to variance.

EDIT: Ah, probably because we center signal around $\frac{B-A}{2}$.

This assumes that both quantization error and original signal are uniformly distributed and have zero mean.

That's a good assumption for quantization noise, but it feels very arbitrary for the signal.

The variance will be proportional to the power (not energy), ONLY if the signal is mean free. Power is defined roughly as the "mean of the squares"  and variance as the "mean of the squares after subtracting the signal mean".


What is the interpretation of Fourier Transform containing only imaginary part? The FT of a unit step function is taken as:

$$
X(\omega) = \int_0^\infty e^{-j\omega t}dt = \frac{-1}{jw}e^{-j\omega t} \Biggr |_{0}^{\infty} = \frac{j}{\omega}
$$

The transform only has the imaginary component and nearly (actually) all literature that I have read deem this transform as non-existent because of the missing real component.

Unfortunately, I have not found an intuitive explanation that describes why this transform is invalid/non-existent.

The integral doesn't exist in the conventional sense. How did you evaluate the limit $\lim_{t\to\infty}e^{-j\omega t}$?

If one allows distributions (generalized functions), such as the Dirac impulse (and its derivatives), as Fourier transforms then the class of functions that can be transformed becomes much wider, and the unit step is one of the functions that has a Fourier transform in that sense. However, its transform is not given by the expression in your question. The Fourier transform of $u(t)$ is

$$\mathcal{F}\{u(t)\}=\pi\delta(\omega)+\frac{1}{j\omega}\tag{1}$$

Take a look at this answer for a derivation of $(1)$.

A simple sanity check also shows that the expression in your question cannot be the transform of $u(t)$. A purely imaginary Fourier transform corresponds to an odd time domain function, which $u(t)$ is not. The expression in your question is actually a valid Fourier transform, but of a different function. Its inverse transform is $x(t)=-\frac12\textrm{sgn}(t)$, which is obviously an odd function.