Timbre vs sound via oscilloscope My understanding is that a distinct sound will generate waves of identical amplitude and frequency, but that most things (such as people's voices, music intsruments, etc) produce multiple distinct "sounds" (waves) at the same time.
My understanding is that the base wave generated by something that does this (something like a human voice) is called the fundamental, and that all the other "secondary" sound waves/signals generated are referred to as the harmonics. And that the set of waves as a whole (the fundamental and any accompanying harmonics) is referred to as the timbre.
Hence, my understanding is that the "sound" we hear with our ears is really the timbres of the things making the sounds (some type of aggregation/composition of its fundamental plus any other harmonic waves generated).
So to begin with, if my understanding above is incorrect or misled, please begin by correcting me! Assuming my understanding is more or less correct however...
In this video of a digital oscilloscope in use you can see (around 0:40 seconds in) that when the author speaks, only a single wave is produced. Whereas, given my understanding of "sound" as being the timbre, I expected to see multiple waves displayed on the screen (one for the fundamental and any others for the respective harmonic waves).
My question: is this oscilloscope displaying the author's "timbre" wave (some type of aggregating or composition of his voice's fundamental plus any accompanying harmonic waves)? Or is it somehow filtering out the harmonics and only showing his fundamental wave? Or is my understanding of sound waves inherently flawed?

Sorry, your understanding is mostly wrong.

Sound is a pressure & velocity wave in a fluid medium, i.e. air.
A human speaker will generate a 3 dimensional sound field around them. The sound pressure and particle velocity will be different at each location.
A microphone measures the "sound" at one location. Depending on the type of microphone it either picks up the pressure, the particle velocity or a mixture of both. Let's stay with sound pressure to keep it simple.
The microphone converts the local sound pressure into voltage. That is the signal that you see on the scope in the video.
Human speech consists of different types of sound. There are three major types
Voiced sounds. These are sounds that have an audible pitch. This includes vowels and nasals, such as "nnnn", "mmm", etc. Everything you can "sing" or "hum" to make a melody. They can be as long as you want.
Fricatives: these don't have any pitch and are more like noise. Example "ssss", "ffff", "shshsh" etc. They can be as long as you want.
Plosives: these are sounds that are short and just a single spike. "p", "k", "t" etc.

Your whole notion of fundamental + harmonics only applies to "voiced sounds", i.e. things that have a pitch. These can indeed be modelled by a fundamental and harmonics and, for example, the location of the harmonics determines whether a vowel is an "aaaa" or an "eee".
What you see on the scope is the sum of the fundamentals and the harmonics (and everything else: noise, reverb, air condition, ...). Technically it's not a wave, it's a signal.
Physically speaking: a wave is something that propagates over both over space and time.
The microphone has already resolved space, since it samples the sound wave only at one location and the dependency on space is lost. While it's still in air, it's a wave, but once it got picked up by a microphone it's just a signal. This is partially just "semantics", many people would still call it a "waveform" since it's a signal derived from a wave.
EDIT: answering the questions from the comments

(1) You speak of particle velocity; my understanding is that sound travels at 343 m/s in air. Are you saying this is not a constant, and that different waves/signals produce sound traveling at different speeds (lets say air is the fluid medium)?

No. Particle velocity and wave propagation speed are different things. Particle velocity is simply the speed with which  air molecules wiggle back and forth. The higher the intensity, the more they wiggle and the higher the velocity. Wave propagation speed is constant and doesn't change with intensity. It's a function of the compressibility and density of the media.

(2) My understanding is that pitch is the frequency of a wave/signal;

No. Pitch is perceived. If a human can tell whether something is "higher" or "lower" it has pitch. For signals like sine waves and periodic signals, pitch correlates highly with the period of the signal. Other signals have no discernable pitch.

wouldn't fricatives and plosives (which are sounds) also have amplitudes and frequencies, hence, pitches?

No. They don't have pitches since they are not sufficiently periodic.
(3) What would be the correct definition for "fundamental" and "harmonics" here?
Any periodic signal can be expressed as a sum of the sine waves that are integer multiples of the fundamental frequency. That's the whole concept behind the Fourier series.

(4) when you say "What you see on the scope is the sum of the fundamentals and the harmonics..." what would you call this signal? Is there a special name for it?

Not really. It's just a signal.

(5) whatever the name of this signal is, I'm intrigued by your use of the word "sum" here; that its the "sum" of all these other things. Again, I would have expected each distinct "sound" to form its own signal on screen. Any chance you could elaborate as to how all of these things gets "summed" into a single signal?

I think you have it backwards. The sound signal just is what it is: it's air molecules wiggling back and forth and the microphone detects whether the pressure goes up and down. The concepts of fundamentals and harmonics (etc) can be used to analyze the signal or break it down into different components, but it's not the way it's generated or propagates. It's more of a math tool to make predictions and analysis easier.


Use wavelet transform to extract a waveform of certain frequency how to apply wavelet transform method to extract alpha waveform( higher and lower band frequency known) from a given known signal ?

The mother wavelet is a bandpass filter. This means we have to apply the CWT(in this case), calculate the correspondence between scales and frequencies, zero out everything outside the known frequency band for the alpha wave and invert the transform.

data = Import["your_sample_eeg_data"];

ListLinePlot[data]

cwd = ContinuousWaveletTransform[data, MorletWavelet[], SampleRate -> 200, Padding -> 0]

freq = (cwd["SampleRate"]/(#1 cwd["Wavelet"]["FourierFactor"])) & /@
   (Thread[{Range[cwd["Octaves"]], 1}] /. cwd["Scales"]);

ticks = Transpose[{Range[Length[freq]], freq}];

WaveletScalogram[cwd, Frame -> True, FrameTicks -> {{ticks, Automatic}, Automatic},
 FrameLabel -> {"Time", "Frequency(Hz)"}, ColorFunction -> "SunsetColors", ImageSize -> 500]

So, we are interested in the 5th octave(the region between 10.51 and 5.25, of course, roughly speaking)

mwd = WaveletMapIndexed[#1 0 &, cwd, Except[{5, _}]]

WaveletScalogram[mwd, ColorFunction -> "SunsetColors"]

It's like performing a brain surgery with a chainsaw ;__;

ListLinePlot[InverseContinuousWaveletTransform[mwd]]

Those are the basic steps I would perform if I have to extract a signal within a certain frequency band.


Audio Separation of .wav signal What I am attempting to do is take a signal (.wav audio file) which contains calls from different Bat species. Here is a graph representation of the signal:

What I have done so far is:

1) Emphasise the signal

2) Compress the signal

This is the output I got plotted:

I need to remove the signals that do not warrent consideration. At the moment, the signal is in a 1D representation, and, after the Separation has been completed the signal will be a 2D representation containing JUST the calls, thus I can perform analysis on each of the calls.

This is my suggested methodology to do this:

1) Split the signal into blocks of: 256x100 (so they overlap)

2) Calculate the RMS (Root-mean-square) of each of the blocks

3) Compare (2) to a threshold value, and, if the block warrants consideration then this can be classed as a call.

But the problem is this:

1) Will splitting the signal into blocks of 256x100 go over the calls and thus making some of the calls lost when calculating the RMS?

2) How do I therefore pick a suitable threshold value? Would finding the peaks in the signal be an accurate assumption for this?

I hope someone can help me :)

If you are using Matlab, the plots looks like Matlab, you can use indexing directly to get a threshold value:

% 'time' is the time series
% 'signal' is the signal itself
idx = signal > 0.2;
plot(time(idx), signal(idx));

You are not clear on what you're trying to achieve, but there could be some information found in the frequency spectrum, so that you could do for example a fourier transform on the signal.


Finding autocorrelation matrix of an autoregressive process AR(1) Having that $\ v(n) = [x(n),x(n-1),x(n-2))]^T $, and being $\ x(n) $ an autoregressive process AR(1) with known variance $\ \sigma_v^2 $ and transfer function $\ H(z) ={ 1 \over {1-0.7z^-1}}$, how could I obtain the autocorrelation matrix $\ R_{vv} = E[v(n), v^T(n)] $ ? Any help is appreciated!

Hi: You're expression will be a 3 by 3 matrix with the 1's along the diagonal and the correlations in the respective places.

$\begin{matrix}  1 & \sigma_{12} & \sigma_{13}  \\
                 \sigma_{21} & 1 & \sigma_{23}  \\
                 \sigma_{31} & \sigma_{32} & 1
\end{matrix}$

The 12 and 21 elements are $\rho$   because the lag one autocorrelation of an AR(1) with parameter $\rho$ is $\rho$.

The 23 and 32 elements are also $\rho$  for the same reason as above:  The lag one autocorrelation of an AR(1) with parameter $\rho$ is $\rho$.

The 13 and 31 elements are $\rho^2 $  because the lag 2 autocorrelation of an AR(1) with parameter $\rho$ is $\rho^2$.

The transfer function implies that $\rho = 0.7$. This is because

$ y_t = \rho \times y_{t-1} + \epsilon_t$ has transfer function $H(z) =
\frac{1}{(1-\rho z^{-1})}$.


How to find the inflection point in a noisy curve? I have the the noisy curve defined by numpy 2D array:

As you can see, it has the first flat segment, then rise, peak and decay phases. I need to find the starting point of the rise phase, marked here by the red dot. What algorythm can I use for that? I need it to be as noise-tolerant as possible.

I think Jazzmaniac's suggestion concerning Savitkzy-Golay filters makes a lot of sense. However, I do not suggest to simply smooth your data using these filters. What many people are not aware of is the fact that there are differentiating Savitzky-Golay filters, which smooth and compute the derivative at the same time. Such filters are implemented in Matlab: check this

You can use the smoothed derivative to detect your desired starting point. It will still require some heuristics and tweaking, but I do not see any other more straightforward way to do it.


A basic question regarding frequency analysis of an EEG signal Assume that an EEG signal is sampled at $f_s = 300$ Hz then a 10000-point segment of it is selected, called $x[n]$. The corresponding 10000-point DFT is then computed and called $X[k]$.
Assume further that famous EEG frequency bands are known as the following intervals (in Hz):
$ \delta=[1,4],$ $ \theta = [4,7],$ $\alpha = [7,12],$ $ \beta = [12,30]$.
Please help me decide which of the following 4 statements are true, and which are false?

Any 5000-point of X[k] is enough to obtain the entire frequency information.

My attempt: According to the famed Nyquist criterion, the maximum frequency we can have is $f_m = f_s / 2 = 150$Hz. The corresponding index $k$ can be found this way:
$$ \frac{k f_s}{N} = f_m = \frac{f_s}{2} \rightarrow k= \frac{N}{2} = \frac{10000}{2} = 5000.$$
This implies that the maximum frequency happens at the 5000'th frequency index. But it doesn't mean any 5000-point segment gives the frequency information so this statement is incorrect.

$X[k]$ for 1000<k<1500 is enough to get the $\alpha$ band information.

My attempt: At $k=1000$ we get
$$ f = \frac{  k f_s }{  N } = \frac{  1000 \times 300 }{  10000 } = 30 \textit{Hz}$$
and at $k= 1500$ we have
$$ f=\frac{ 1500\times 300  }{10000   } = 45 \textit{Hz}$$
which are out of the $\alpha$ band range, so this statement is incorrect as well.

$X[k]$ for 8000<k<10000 is enough to get the $\alpha,\beta,\delta, \theta$ bands information.

My attempt: For the very same seasons just discussed, $k=8000$ corresponds to $f= 240$ Hz and $k= 10000$ corresponds to $f= 300$ Hz  which are out of all bands' ranges. As a matter of fact, 240 and 300 Hz are above the maximum frequency $f_m$ allowed by Nyquist's criterion, so I don't really know how to interpret this. Maybe I'm mistaken?

A lowpass filter with frequency $ \frac{\pi}{2}$ removes power line frequency at $60$ Hz.

My attempt: I think this is trivially correct because $\frac{\pi}{2}$ Hz is about 1.57 Hz and the lowpass filter zeros out anything above this frequency, including 50 Hz.

The first $5000$ points (well, $5001$ to be exact) contain all frequency information, because the signal is real-valued and, consequently, the spectrum is conjugate symmetric:
$$X[k]=X^*[N-k]\tag{1}$$
where $N$ is the DFT length. But from $(1)$ we cannot infer all values of $X[k]$ from any arbitrary collection of $N/2$ (consecutive) points of $X[k]$. So you're right.

You got that right.

This is a bit tricky, and you'll need Eq. $(1)$. Think about which values of $X[k]$ you can derive from the index range $k\in[8001,9999]$. It should turn out that in fact you got everything you need.

They should specify what they actually mean by $\pi/2$. I would expect this to be a normalized frequency in radians, i.e., $2\pi f/f_s$, where $f_s$ is the sampling frequency. If that is the case, then $\pi/2$ is just half the Nyquist frequency. Draw your own conclusions.


Who first understood the importance of poles? Who first understood (or at least published papers on) the importance of poles in understanding transfer functions in the frequency domain?

If I had to guess, I'd suggest Nyquist or Bode but I know very little about the history of this particular subject.

If you consider poles of an integral transform domain to be important to the solution of differential equations: (as usual,) Euler did it first, 1753.
One "importance" of poles is that they're part of a very useful representation for linear systems. They must've appeared when people started looking at functions as built from generating functions, so that'd date that .
So, both Bode and Nyquist are probably 100 years too late: this would emerge from looking at systems like Laplace did, around 1785.
If you think of poles as something that is relevant to your analysis of holomorphic functions, then clearly Cauchy would be your guy, starting with "Sur les intégrales définies", 1814, where he looked at simple poles within an area of integration.
When you really look into where the residue theorem (which really concerns itself with integrals over functions with poles, and is immensely important to a lot of things that we can do with the Fourier transform) was formed, you'd be looking at Lindelöf's "Le calcul des résidus et ses applications à la théorie des fonctions", 1905.
Finally, the Nyquist stability criterion was published when he was working at Siemens, 1930.

     Laplace     Lindelöf
  Euler | Cauchy    |
    |   |  |        |  Nyquist
    v   v  v        v  v
|---------.---------.---------.-->
1         1         1         2
7         8         9         0
0         0         0         0
0         0         0         0


Bandwidth of a finite impulse response filter In a machine learning application, a model learns to apply a filter h(n) via convolution to a 1-dimensional input signal s(n) (e.g. sampled at Fs = 16 kHz). The filter size (filter length?) is limited to a fixed number, e.g. k = 100, so the model cannot learn a filter that exceeds this duration but it is free to insert zeros and make the filter "shorter".

Does this mean that the filter size limits the duration of the "slowest" signal component that falls into the passing bandwidth and hence Fs/k = 160 Hz is the minimal bandwidth that h(n) can represent (assuming it's a band pass filter)?

Yes, the length of a filter limits the "frequency-selectiveness". If you remember modern physics, that's mathematically very related to the Heisenberg uncertainty principle: something that is temporally short can't be very sharp in the frequency domain, and vice versa.

You can find one example of such a constraint in this answer, though the formula only applies to specific filter constraint.

Generally, filter design is one of the oldest application fields of evolutionary / ML algorithms. You'll be delighted to find that even some of the very "classical", non-ML algorithms for bandpass filter design make extensive use of a loss / reward function while optimizing a filter.

If you're doing modern machine learning for filer design, there's a lot of modern literature, but I'd really start by reading up on the Remez exchange (or just Remez) algorithm, to give you a starting direction. Understanding that relatively simple algorithm and its convergence will give you a solid base!

Especially when later writing a report on it, that would be the most basic benchmark for your ML approach that I could think of. Compare yourself to that!


What is the variance of DFT of Fourier coefficient of difference of a vector of white noise? Consider $\big\{x[0], x[1], \ldots, x[N-1]\big\}$. Suppose,
\begin{cases}
x[n] \sim \mathcal N\left(0, \sigma^2\right)\\
\big\langle x[n], x[n-1]\big\rangle = \frac12 & \forall \ n\\
\big\langle x[n_1],x[n_2]\big\rangle = 0 & \text{whenever}\quad\lvert n_1-n_2\rvert > 1
\end{cases}
I need to compute $\operatorname{Var}\{X[k]\}$ where,
$$
X[k] = \sum\limits_{n=0}^{N-1} x[n] \ e^{-i 2 \pi \frac{nk}{N}}.
$$
The definition of variance I am using is
$$
\operatorname{Var}\big\{X[k]\big\} = \operatorname{E}\big\{\big\lvert X[k]\big\rvert^2\big\} - \operatorname{E}\big\{\big\lvert X[k]\big\rvert\big\}^2
$$
The second term is zero. But I am failing to compute the variance. I am not sure how to do it.

This appears to be a homework/quiz question; we don’t provide full solutions but the following hint should help:
All samples in a white noise process are independent and identically distributed (IID) so it is helpful to understand how the variance grows as you sum N IID random variables: the sum of $N$ IID random variables with variance $\sigma^2$ will have a variance of $N\sigma^2$. Note that the exponential in the DFT equation only changes the phase of each sample without changing the amplitude. Importantly, all the samples were independent from each other and identically distributed; how does this change when we change the phase of those samples?
If we roll 10 fair die, the face value of each die will be independent. If we then turn over each die by same or different rotations, would the new values also be independent?


How to explain IDFT of square root of a DFT? I have this MATLAB code:

w = [1 2 3 2 1];
W=fft(w);
V=sqrt(W);
v=ifft(V);

And I need to explain how to get to $v[n]$ without using the direct formula for DFT - I need to use the properties of DFT.

I've tried solving the problem using circular shift and circular convolution but I can't seem to understand the result.

According to MATLAB, $v[n]=\textrm{[1 1 1 0 0]}$.

Since $V(\omega)V(\omega)=W(\omega)$, in the time domain we expect $v[n]*v[n]=w[n]$. That is, the convolution of $v[n]$ with itself should become $w[n]$.

$w[n]$ is in the form of a triangular signal. A triangular signal of odd length $N$ can be constructed by convolving a half-rectangular signal of length $(N+1)/2$ with itself. So $v[n]$ is a rectangular pulse with $3$ elements equal to one and the remaining $2$ elements equal to zero.


Question regarding transfer functions and prerequsities for finding the real impulse response The transfer function of a system is given by:

$$\large H(s) = \huge \frac{V_{out}(s)}{V_{in}(s)}$$

In digital domain the principle is of course the same, just replace laplace transform with z-transform and voltages with amplitudes. The transfer function, when converted to time domain, solves the impulse response of the system when conditions are met.

My question is about the prerequisites for finding the real impulse response using this method (other than LTI properties). Obviously, if a signal with a shorter length than the impulse response is used for the analysis, the real impulse response can't be solved for. However, is it enough to have a signal the length of the impulse response, regardless of what the input is? (Using dirac delta as input will of course solve the impulse, but what about other cases such as $sin(t)$?)

I can at least answer the question when it comes to fourier transform and related frequency response $H(w)$. In that scenario the signal indeed needs to be zero padded properly or we will have circular convolution. I am assuming the case is then the same for laplace transform. Thus, is it correct to say that for the $H(s)$ to yield the correct transfer function, the time signal representing the input needs to be zero padded by the length of the impulse (-1)?

If I understand your question correctly, it is about measuring the impulse response or frequency response of an actual system.

One of your questions appears to be "can I use any input signal to measure the system's frequency response (or impulse response)?" The answer is clearly no. As an extreme example, take a sinusoidal input signal $x(t)=\sin(\omega_0t)$. An LTI system with frequency response $H(\omega)=|H(\omega)|e^{j\phi(\omega)}$ will have a response $y(t)=|H(\omega_0)|\sin(\omega_0t+\phi(\omega_0))$, where $\phi(\omega)$ is the system's phase response. Consequently, with a sinusoidal input signal you can just measure $H(\omega)$ at a single frequency (namely the frequency of the input signal). That's why an impulse is a good input signal for measuring the properties of a system: it contains all frequencies, not just a few.

One common practical method to measure a system's impulse response (or frequency response) is the sine sweep method, where the frequency of the input signal is swept over the whole frequency range of interest. This method and a few others are compared in this paper.


How to transform data to uniform distribution (uniform percentiles)? Given the list of data points (normalized in [0,1] range), I plot the histogram of values and compute percentiles (shown as x ticks).

How to find a transformation of data values so the histogram is approximately uniform. Which would, in turn, make percentile values also uniformly distributed.

Hi: You can calculate the empirical cumulative distribution of the data. By this, I mean, given some observation in the sample, $x_i$, calculate $P(X < x_{i})$ by calculating the proportion of observations that are  less than $x_{i}$ (i.e. the percentiles ). Then, do this for all the $x_{i}$ so that you have the cumulative distribution of the $x_{i}$.

Then, $P(X < x)$ is uniform for a given value of $x$.

In fact, it seems like you already did this but the percentile values should be on the vertical axis and the values of the data should be on the x-axis.

Note that page 14 of this PDF explains the concept more clearly than I have.

Example Implementation

Below is a quick-and-dirty attempt to illustrate this answer.  The image below shows the original histogram of the Gaussian, the empirical cumulative distribution function of that data, and then the histogram of the converted data.

R Code Below

par(mfrow=c(3,1))
# First, generate some Gaussian numbers.
gaussian <- rnorm(1000,0.0,0.05)
gh <- hist(gaussian, breaks=1000)

empirical_cumulative_distribution <- cumsum(gh$counts)/1000

plot(gh$mids, empirical_cumulative_distribution)

uniformize <- function(x) {
  ans_x <- x
  for (idx in seq(1,length(x))){
    max_idx <- max(which(gh$mids < x[idx]))
    ans_x[idx] <- empirical_cumulative_distribution[max_idx]
  }
  return(ans_x)
}

uniform2 <- uniformize(gaussian )
hist(uniform2, breaks=100)
par(mfrow=c(1,1))


FFT of a $N$-length real sequence via FFT of a N/2-length complex sequence I would like to optimise FFT algorithm by only using a FFT of size $N/2$ to transform a real sequence of length $N$. There are already many questions and many excellent answers on the matter (here for example).

I've been trying to follow the linked answer to implement the splitting algorithm in Matlab and I checked against the standard FFT Matlab routine. I found a shift of 1 sample in the spectrum (sine input) and a weird peak at a (seems to me) random frequency. I would really appreciate a review of my code or simply a hint to debug the problem.

As a side note, I'm trying this route because I need a memory efficient FFT for an embedded device.

close all;
clear all;
clc;

%% Parameters
sample_rate = 50; % Hz
frequency = 10; % Hz
N = 256;
p = 1:N;
t = (p-1) / sample_rate;

f = sample_rate*(0:(N/2 -1))/N; % frequency vector

%% Generate real sequence
data = sin(2*pi*p*frequency/sample_rate);

figure();
plot(t, data);
title('Input');

%% Split even-odd
% Even indeces are odd because
% Matlab starts indexing from 1
idx_even = 1:2:N;
idx_odd = idx_even+1;
data_even = data(idx_even);
data_odd = data(idx_odd);
data_cpx = complex(data_even, data_odd);

%% Perform FFT of length N/2
X = fft(data_cpx, N/2);

%% Perform Butterfly
k = 0:(N/2 -1);
W = exp(-2*pi*1i*k/N);

X_even =     0.5*(X(k+1) + conj(X(N/2 - k)));
X_odd  = -1i*0.5*(X(k+1) - conj(X(N/2 - k)));
G = X_even + W .* X_odd;

%% Perform FFT of length N to check results
X_N = fft(data, N);

%% Plot
figure();
plot(f, real(X_N(1:128)));
hold on
plot(f, real(G));
title('Real');

figure();
plot(f, imag(X_N(1:128)));
hold on
plot(f, imag(G));
title('Imag');

figure();
plot(f, abs(X_N(1:128)));
hold on
plot(f, abs(G));
title('Mag');

figure();
plot(f, angle(X_N(1:128)));
hold on
plot(f, angle(G));
title('Phase');

The following Matlab Excerpt corrects your mistake. You had a problem with the DFT length N of the composed sequence... Note that I used a modulus operator for simplicity, you should replace it with indexing for efficiency.

%% Perform Butterfly

%k = 0:(N/2 -1);
%W = exp(-2*pi*1i*k/N);
%X_even =     0.5*(X(k+1) + conj(X(N/2 - k )));
%X_odd  = -1i*0.5*(X(k+1) - conj(X(N/2 - k )));
%G = X_even + W .* X_odd;

k = 0:N-1;
W = exp(-2*pi*1i*k/N);
X_even =     0.5*(X(mod(k,N/2)+1) + conj(X( mod(- k,N/2)+1 )));
X_odd  = -1i*0.5*(X(mod(k,N/2)+1) - conj(X( mod(- k,N/2)+1 )));
G = X_even + W .* X_odd;

fN = sample_rate*(0:N-1)/N;   % new frequency vector data length N...

Mathematically speaking, the line where you made the implementation error is the Even-Odd combination stage:

$$ X[k] = X_e^{N/2}[k] + e^{-j\frac{2\pi}{N}k } X_o^{N/2}[k] ~~~,~~~\text{ for } k=0,1,...,N-1 $$

Where $X_e[k]$ and $X_o[k]$ are $N/2$ point DFTs of the even and odd indexed samples of the original length $N$ real sequence $x[n]$. And $X[k]$ is the N-point DFT of $x[n]$. Note that frequency index $k$ ranges through $0,1,..,N-1$. But those $X_e[k]$ and $X_o[k]$ are half length (but periodic!) sequences. Hence for $k> N/2-1$ they repeat. Remeber that DFT seqeuences are inherently periodic: $X_e[k+N/2] = X_e[k]$. In order to use the base length of those half-length periodic sequences, you can utilize the modulus operator:

$$ X[k] = X_e^{N/2}[(k)_{N/2}] + e^{-j\frac{2\pi}{N}k } X_o^{N/2}[(k)_{N/2}] ~~~,~~~\text{ for } k=0,1,...,N-1 $$

where $(k)_{N/2} = \mod(k , N/2)$ is the modulus of $k$ wrt $N/2$.
Argument of every DFT sequence should circularly wrap around...


Amplitude and Phase of FFT gives messy results - Matlab I'm trying to calculate the spectrum magnitude and phase of a signal, but when I plot them I get pretty messy results:

The signal is displayed on the right, then I have the magnitude and phase of the fft in the centre and to the left accordingly.

I would have expected to see just few peaks in both the spectrums but instead there are a lot of peaks and it gets really messy.

This is my code:

T = mean(diff(t))
fs = 1/T

figure;
subplot(1,3,1)

plot(1000*t,sig);
title ('Time domain');
xlabel('Time(ms)')
subplot(1,3,2)
plot((0:9999)/10,abs(fft(sig))); % magnitude spectrum
title('Spectrum');legend({'Magnitude'});
xlabel('Frequency (Hz)')
subplot(1,3,3)
plot((0:9999)/10,angle(fft(sig))); % phase spectrum
title('Spectrum');legend({'Phase'});
xlabel('Frequency (Hz)')

Am I missing something or is it supposed to be like this?

I tried already rounding off to zero values that are close to zero, but that did not change anything.

If you round off tiny magnitudes to zero (to remove numerical noise), you should also zero the matching phase result bins.  The phase of noise is usually just random.


Interpretation of fft plot I am trying to understand how to interpret a spectrogram image plot: (1) what does the yellow band and its corresponding color code with a high positive value indicate? There are so many yellow bands in comparison to blue. The legend shows that the yellow bands have higher positive values, so is that higher frequency?
(2) I used the following code to generate the spectrogram. The X axis contains 50 data points because the method splits the signal into segments of length 50. I arbitrary used 10 into the spectrogram function. Is there a specific rule as to how to determine the number of points?
x(1) = rand;
for n = 1:500
    x(n+1) = 4*x(n)*(1-x(n));
end

specto=spectrogram(x,10,0,256); %using 256 fft points with no overlap

logspecto=(20*log10(abs(specto)));

I'm assuming you're using Matlab. The Matlab function and its parameters are:
spectrogram(x,window,noverlap,nfft)

The color indicates the power at a given frequency. Typically, this is given in dB. In your plot, the power is concentrated in the low frequencies a the bottom of the plot. The higher frequencies at the top of the plot have less power.

The value of 10 you used is the width (number of samples) of your windowing function. It takes a vector of n samples and applies weights to the samples such that the ends of the vector taper to zero. This is done to limit leakage (a loss of frequency resolution/introduction of artificial frequencies) during the DFT process. You choose the width of your window according to the frequency resolution you desire (more samples give a finer frequency resolution).

I hope this helps.


How could I remove the central red line How can I remove the central line; there is someone told me use high-pass filter, is that right? how to use it in matlab

Just apply:

[B,A]=butter(3,2*f0/fs,'high');
y=filter(B,A,x);

where fs is your sampling frequency, f0=5 (your red band-width, as shown approximatedly in your spectrogram), x is your nx1 column vector data and y is your new filtered data.

A | B is the denominator | numerator of your highpass filter.

Hence the original picture:

spectrogram(y);

will show the 0-5Hz band in a pure blue band, a sign of deleting everything in there...

Best regards....


Why the correct probability for a square M-QAM is equal to the square of correct probability of ($\sqrt M$)-PAM? The question is why:

$P_c (\text{square} $ $M-QAM) = (1 - P_{\sqrt M})^2$

where: $P_{\sqrt M}$ is the error probability for $\sqrt M$-PAM signal.
Thanks.

The decision in M-QAM is actually two independent sub-decisions:

(i) what is the horizontal coordinate of the constellation point nearest to the received signal point?

and

(ii) what is the vertical coordinate of the constellation point nearest to the received signal point?

Consequently, the decision is correct if and only if both
sub-decisions are correct. But the probability of a sub-decision being correct is just $(1-P_{\sqrt M})$, and so the probability that both
are correct is $(1-P_{\sqrt M})^2$ by independence of the events.  If
the I and Q branches of the receiver had mismatched gains, the independence would still hold but the two probabilities being multiplied
would be different (corresponding to the different SNRs in the two
branches).


I don't really understand what are you trying to achieve with vector B that has 6 elements. Pre-emphasis filter is defined as:

$$y[n]=x[n]-\alpha x[n-1]$$

Where $\alpha$ is $0.97$ in your case.

You can see that there are only two elements in vector B, right? So the code should be something like that:

B=[1, -0.97];
xd = filter(B, 1, data, [], 2);

You can verify the frequency response:

freqz(B, 1)


Capacity of cascade binary symmetric channels Let's imagine that we have interconnected in cascade $L$ binary symmetric channels each with the same transition probability $p(y|x) \in \{p, q=1-p\}$, where the output of each BSC is connected to the input of the next.

The overall channel is also a BSC channel, but what happens as $L \to \infty$?

A binary symmetric channel (BSC) can be characterized by its complemented probability $p$. Its well-known capacity is

$$C = 1 - H(p) = 1 - (-p\log(p) - (1-p)\log(1-p))$$

where $H(p)$ is binary entropy function:

A $L-$concatenated BSC, which is also a BSC characterized by $p_L$, can be visualized as in the figure below

The complemented probability $p_L$ is derived

\begin{align}
p_L &= p_{L-1} (1-p) + (1-p_{L-1}) p \\
&= p + (1-2p)p_{L-1}\\
\implies 1 - 2 p_L &= (1-2p)(1-2p_{L-1}) \\
\implies 1 - 2 p_L &= (1-2p)^L \\
\end{align}

Thus
$$p_L = \frac{1}{2}(1 - (1-2p)^L)$$

If $p=0$ then $p_L = 0$, $H(p_L) = 0$ and $C = 1$.

If $p=1$ then $p_L = 0$ or $p_L = 1$ depending on $L$ is pair or impair; $H(p_L) = 0$ and $C = 1$.

If $0 < p < 1$ then $\lim_{L \to \infty} p_L = 0.5$; $H(p_L) = 1$ and $C = 0$.

Conclusion: if the unit BSC is not certain $(p \neq 0, 1)$, the capacity of infinitely-concatenated BSC tends to $0$.


4th order butterworth filter using frequency response formula I know this may be a kind of basic question, but I have my head wrapped over it for some time and hadn't found a solution.
I want to generate the frequency response of a 4th order butterworth filter programatically for a given sampling rate and cutoff frequency. I've used the matlab filter response and my filter code response to compare. In my code i've used:

$$
H(j\omega )=\frac{1}{(1-3.4142136\hspace{1.2mm}\omega^{2}+\omega^{4})+j2.613126(\omega - \omega^{3})}
$$

And got a close comparison between the two. What am I missing to make my code response like the matlab one? adjust the even function symmetry on the magnitude and odd function symmetry on angle at nyquist frequency and rolloff of both using the H(jw) programatically.

Here is my matlab code and results

clear all;
clc;
n = 4;
fs = 100;
t = 1:1:100;
lo = 35;
in = ones(1,length(t));
in = ifft(in);
[b,a] = butter(n,lo/(fs/2),'low');
outA = filter(b,a,in);
outB = ones(1,length(t));
for x = 1:length(t)
    wl = x/lo;
    outB(x) = (1/((1-3.4142136*(wl^2)+(wl^4))+(2.613126*(wl - wl^3))*1i));
end
fourierabsA = abs(fft(outA));
fourierabsB = abs(outB);
fourierangA = angle(fft(outA));
fourierangB = angle(outB);
subplot(2,2,1);
stem(t,fourierabsA);
title('Matlab Filter Magnitude')
subplot(2,2,2);
stem(t,fourierabsB);
title('My Filter Magnitude')
subplot(2,2,3);
stem(t,fourierangA);
title('Matlab Filter Angle')
subplot(2,2,4);
stem(t,fourierangB);
title('My Filter Angle')

Matlab's fft function returns the two-sided spectrum. To get your graphs to agree, start by only looking at the first half of the spectrum.

% note *2 to double the length
in = ones(1,length(t)*2);
in = ifft(in);
[b,a] = butter(n,lo/(fs/2),'low');
outA = filter(b,a,in);
fftA = fft(outA);

% snip out the first half
fftA = fftA(1:length(fftA)/2);

fourierabsA = abs(fftA);
fourierangA = angle(fftA);


Comparing two sound files in embedded-systems (C language) what would be the most efficient ("cheap" in terms of system resources etc.) way to compare 2 sound files?

Here is the scenario. Sound A comes as an input. There is sound B stored in memory. We need to see if sound A is identical to sound B. We know that both sounds can be copies of the same sound file, but there is a chance that sound B was processed to reduce the size.

What would be the most efficient ways to compare these sounds and see if they are similar (or rather the same).

ADDED

"Cheap", architecture is Cortex ARM M0. 16kb ram, 256kb flash.
"Identical" means that these 2 sounds come from the same sound file
(same exact sound effect). THe difference may be that sound B (stored
on the microcontroller side) may be processed to reduced size.
Basically imagine we have 10 sound files, we copy N7 to MCU (maybe
process it to reduce size) and then play them and expect MCU to
recognize the sound file when it's played.

The kind of sound - technically it can be any sound (within the frequencies recognizable by humans). Sound effects (explosions, weapons firing, tires screeching), people saying words, melodies played by musical instruments. But the common thing is that they are all identical (instruments playing melodies would be same instruments same room, same time, recorded with the same microphone etc.).

Worth mentioning (to avoid confusion), that the sound is being sent to MCU from a sound output, not as a file.

Matched filtering, or cross-correlation, is a pretty cheap and easy way to determine the similarity between signals. Otherwise you could try implementing the algorithm Shazam uses.


What is this kind of distortion called? Today while watching a movie I noticed the following distortion:

Note. The red ring is me painting over to show the area of interest.

There are certain areas in which I can see non-smooth shades of different colors. What is this distortion called and what is causing it?

This is the video meta information in case it helps:

Video #1
      ID : 1
      Format : AVC
      Format/Info : Advanced Video Codec
      Format profile : High@L4.0
      Format settings, CABAC : Yes
      Format settings, ReFrames : 4 frames
      Codec ID : V_MPEG4/ISO/AVC
      Duration : 43mn 1s
      Bit rate : 5 109 Kbps
      Width : 1 920 pixels
      Height : 1 080 pixels
      Display aspect ratio : 16:9
      Frame rate mode : Constant
      Frame rate : 23.976 fps
      Color space : YUV
      Chroma subsampling : 4:2:0
      Bit depth : 8 bits
      Scan type : Progressive
      Bits/(Pixel*Frame) : 0.103
      Stream size : 1.54 GiB (85%)
      Default : No
      Forced : No
      Color primaries : BT.709
      Transfer characteristics : BT.709
      Matrix coefficients : BT.709

Your image with increased contrast:

Figure 1. Denis Kulagin's figure with gain applied.

That is kind of color banding due to quantization. Quantization errors can be made worse by repeated decoding and encoding, because quantization levels and things that are quantized may not match between the encodings. In video encoding, it is not directly the pixels that are quantized, but things like discrete cosine transform coefficients for 8x8 or so blocks of chroma and luma values of pixels. JPEG images compressed with a low quality setting exhibit color banding for the same reasons:

Figure 2. Top: 60 KB PNG file, Bottom: 2 KB JPEG file

Also color space conversions may introduce quantization errors and banding. For example sRGB and Rec. 709 transfer characteristics differ at worst by about a factor of two at low intensities.


pixels in the same image are correlated but how can I visualize it? I am new in image processing.I know that pixels in the same image are correlated with each other but I don't know how can I draw the relationship b/w these pixels in the form of graph?.

If you would like to show that in natural images, neighboring pixels tend to have similar values, you can use the cooccurrence matrix. The values in the matrix will be distributed mostly along the diagonal.

In the snapshot below, the second image is random and non-correlated. For that image the cooccurrence matrix is pretty much a constant.  I'm not showing the matrices themselves, but a plot of their values such that white means '0'.


How to plot poles, zeros and impulse response of a transfer function I want to plot poles and zeros as well as impulse response of a transfer function that is given below:
$$H(z) = -0.0625z^4 + 0.25z^3 + 0.625z^2 + 0.25z -\frac{0.0625}{z^4}$$
For poles and zeros, I used MATLAB function residuez. However, I got error which is similar in impulse response plot.
code:
 b = [-0.0625 0.25 0.625 0.25 -0.0625];
 a = [0 0 0 0 1];
 [R,p,C] = residuez(b,a); % for poles
 y = impz(b,a); %for impulse response
 figure(1);
 stem(y);

Could someone help.
Thanks

Two hints:

The standard way of writing a transfer function in the $\mathcal Z$ domain is a rational function in $z^{-1}$, not in $z$. Try this
Use tf2zp() instead of residuez()


Difference between Fourier-Transform and FFT of rectangular pulse I'm trying to find a link between the Fourier-Transformation of aperiodic Signals and the FFT of them.
So to start with a basic example, let's take a rectangular pulse with width 0.1s and amplitude of 1 shifted by 0.05.
Using correspondency, i can calculate the expected spectrum:
$X(f) = 0.1 \cdot sinc(0.1f) \cdot e^{j 2 \pi f \cdot 0.05} $

But now, when i generate the Signal with the follwing Matlab-Code:

f_abt = 50e3;
x=0:1/f_abt:1;
y=zeros(1,length(x));
for ii=1:length(x)
    if x(ii)<=.1
        y(ii)=1;
    end
end

And calculate the spectrum of it, the result is depending on the length of the signal. So, when i compute the one-sided spectrum from the signal generated above (1s duration), i get:

Then, when i put the signal length to 2s (everything else unchanged):

x=0:1/f_abt:2;

I get the following spectrum:

I guess the difference comes from the FFT-Algorithm i use. When doing FFT, i normalize the Values by Nfft, so it makes complete sense that my amplitudes change when i change the Signal length.

My question is: How do i get the right spectrum and how do i know it is right, e.g. when i can't calculate it 'by hand' using correspondencies? I'm having issues finding the connection between my "real", time limited signal and its FFT and the "theoretical" rectangular pulse.

Code I use for calculation of one-sided spectrum:

function [f_xa, mag, phase] = calc_fft_f(ta, xa)
N_a = numel(xa);
fft_xa = fft(xa);
P2_norm = fft_xa/(N_a);
if (mod(N_a,2))
    P1_norm_single = P2_norm(1:ceil(end/2));
    P1_norm_single(2:end) = 2*P1_norm_single(2:end);
else
    P1_norm_single = P2_norm(1:(end/2)+1);
    P1_norm_single(2:end-1) = 2*P1_norm_single(2:end-1);
end

mag = abs(P1_norm_single);
phase = rad2deg(angle(P1_norm_single));
Fsa = 1/(ta(2)-ta(1));
f_xa = Fsa*(0:(length(mag)-1))/N_a;

end

Thanks in advance!

Assuming that the relevant portion of a continuous-time signal $x(t)$ is inside (or has been shifted to) the interval $[0,T]$, the DFT of a sampled version of the signal approximates the continuous-time Fourier transform (CTFT) in the following way:

$$\begin{align}X(f)&=\int_{-\infty}^{\infty}x(t)e^{-j2\pi ft}dt\\&\stackrel{\textrm{truncation}}{\approx}\int_{0}^{T}x(t)e^{-j2\pi ft}dt\\&\stackrel{\textrm{sampling}}{\approx}\sum_{n=0}^{N-1}x(n\Delta t)e^{-j2\pi f n\Delta t}\Delta t\tag{1}\end{align}$$

with $T=N\Delta t$.
From $(1)$ with $\Delta t=T/N$ and with $f=k/T$, a sampled version of $X(f)$ can be approximated by

$$X\left(\frac{k}{T}\right)\approx \Delta t \sum_{n=0}^{N-1}x(n\Delta t)e^{-j2\pi k n/N}=\Delta t \cdot X_d[k]\tag{2}$$

where $X_d[k]$ is the length $N$ DFT of $x_d[n]=x(n\Delta t)$.

Note that for time-limited signals, the truncation error can be made zero, and for perfectly band-limited signals, the sampling error can be made zero. Since a signal cannot be time-limited and band-limited at the same time, there's always at least one of the two errors present. In practice, you usually have to deal with both types of errors.

The following Matlab/Octave code shows an example:

Fs = 1e3;   % sampling frequency
Ts = 1/Fs;
T1 = 0.1;
T2 = 2;
tgrid = 0:Ts:T2;
N = length(tgrid);
x = zeros(1,N);
x( find( tgrid <= T1 ) ) = 1;
fgrid = (0:N-1)*Fs/N;

% analytic continuous-time Fourier transform
X = T * sin( T*fgrid*pi ) ./ (T*fgrid*pi) .* exp( -1i*pi*fgrid*T );
X(1) = T;

% DFT approximation
X2 = fft(x,N) * Ts;

plot( fgrid,abs(X),fgrid,abs(X2),'r' )
axis([0,Fs/2,0,T]), grid on
legend('analytic','DFT')


Questions about discrete signal energy calculation I'm looking at the definition of signal energy (e.g. Wikipedia, cnx.org). For discrete signals, it's defined like the following, where $x(n)$ holds the signal:

$
Energy = \sum_{n=-\infty}^{\infty} |x(n)|^2
$

So my questions:

For a windowed, finite signal like double signal[256], the sum is from 1 to 256 (or 0 to 255 in a program) rather than $-\infty$ to $\infty$, right? (I don't even know how I would sum over infinity.)
Why does the energy formula have the absolute-value operator $|...|$? The result of taking the absolute value is squared anyways to produce a positive value, so taking the absolute value seems to be pointless. Is it because the $x(n)$ can be complex, so the absolute value of a complex number would be the scalar from Pythagoras' theorem?

Yes. No need to sum up an infinite number of zeros.
Yes, it is because x(n) could be complex. If we didn't take the absolute values of the complex numbers (Euclidean norm), a signal containing
$[\dots,0,0,a+ia, a-ia,0,0,\dots]$
would have an energy of zero (instead of $4a^2$) although it contains non-zero samples. For real numbers however $|a|^2 = a^2$ and the absolute value doesn't matter.


How to identify useful elements of audio recording and ignore lulls? I assume this is a relatively simple DSP task, but I'm having trouble finding any information on how to approach my problem.  Hopefully someone with more experience and brain power can help me out:

I have many (thousands) audio records of bird calls.  I am looking for a technique that will identify the timestamps within each audio file during which the bird is chirping.  Here's an example of a high quality audio file:

Ultimately, I wish to export a table of the frequencies created by the bird in each file, for further analysis.  I've already created a Python program that will allow me to manually select data points of interest, but I'd like to automate the process.

The simplest process I can think of would be to set a minimum power threshold for data points of interest, but the audio files have varying degrees of background noise.  Some of the audio files may be unusable due to background noise, which is fine, but they need to be flagged as such.

So the two questions are:
1) How to identify low-quality audio files
2) How to identify the useful sections of high-quality audio files

Thanks!

EDIT: I don't know why i thought this question was matlab based, but my solution is written with matlab in mind. a quick search showed my the python scikit has otsu's method for automatic thresholding (this is what matlab uses). Other than that I think most of the code is relatively python safe and can easily be translated. My apologies for my error

from skimage import filter
val = filter.threshold_otsu(camera)

END EDIT

Could you threshold the amplitude of the recording? say this is our recording

where we are just plotting

plot(abs(sound))

can we say the chirping is probably above 2? and assume everything below magnitude=2 is just noise? If so, and assuming you have the imageprocessing toolbox, you could

normalize the audiofile amplitude (between 0-1) needed for next step
get automatic thresholding of file imageprocessing toolbox function
find indicies that exceed threshold and say these are bird chirps

the only problem is with sound even with a strong signal, it will usually cross 0 so you would have discontinuities even in the actual bird chirping. Maybe you can add some sort of hysteresis to the masking operation. Like
if the current sample is below the threshold BUT a large percentage of the previous few samples were above the threshold, we will assume this sample is also an interest point

roughly your code would be something like this (im not near a comp with matlab, so this may have slight syntax errors). This is also not very optimized, but maybe it's a springboard for you or others

%% step 1
%we are using magnitude only
my_min = min(abs(soundfile));
my_max = max(abs(soundfile));

%gets sound file magnitude between 0-1 only
normalized_sound_mag = 1/(my_max-my_min) * (abs(soundfile) - my_min);

%% step 2
%get a threshold, we needed the magnitude between 0-1 for this function to work
sound_thresh = graythresh(normalized_sound_mag);

%% step 3
%this vector will store all sample indexes that are interesting
idx_sound = 0;

num_for_hysteresis = 10;        %the number of samples to use for hysteresis
hysteresis_percentage = 0.8;    %percentage of samples that must be above threshold
                                %so say in our set of 10, 2 samples are below threshold
                                %since .8 are above it, we say they are all of interest

%because the way we implement hysteresis we have to skip the first few samples
for (i=num_for_hysteresis:1:length(normalized_sound_mag))
    %if we are above thresh, without a doubt save the index
    if (normalized_sound_mag(i) > sound_thresh)
        idx_sound = [idx_sound, i];
    else
        %hysteresis check prev samples, creates a boolean vector. 1 means value above thresh, 0 means it was below
        samples_above_thresh = normalized_sound_mag(i-num_for_hysteresis:i) > sound_thresh;

        %nnz, gets the number of nonzero elements in a matrix
        num_prev_samples_above_thresh = nnz(samples_above_thresh);

        %if the number of samples in the prev window met our criteria, this current sample
        %is probably of interest as well
        if (num_prev_samples_above_thresh > hysteresis_percentage * num_for_hysteresis)
            idx_sound = [idx_sound, i];
        end
    end
end

%idx_sound should now have all the indicies f interest, these can now be used
%directly on the original soundclip

we could also use the frequency power and threshold, rather than the sound amplitude, the basic outline would still be the same.


What are the ranges of values for each audio feature/descriptor? I need to normalize some audio feature measurements I've done with the software Wolfram Mathematica.

Coming from a non signal processing background, I have some difficulties understanding what several audio features are actually measuring and thus the range of values in which each measurement is expressed.

I can't seem to find any documentation regarding the "physical" range of these descriptors, both in the Wolfram software and the MPEG-7 standard.

This is the list of features analyzed, using the names given to them by Mathematica:

{"Max", "MaxAbs", "Min", "MinAbs", "MinMax", "MinMaxAbs", "Mean",
"Median", "StandardDeviation", "Total", "Power", "RMSAmplitude",
"Loudness", "CrestFactor", "Entropy", "LPC",
"PeakToAveragePowerRatio", "TemporalCentroid", "ZeroCrossingRate",
"ZeroCrossings", "FundamentalFrequency", "Formants",
"HighFrequencyContent", "MFCC", "SpectralCentroid", "SpectralCrest",
"SpectralFlatness", "SpectralKurtosis", "SpectralRollOff",
"SpectralSkewness", "SpectralSlope", "SpectralSpread",
"ComplexDomainDistance", "ModifiedKullbackLeibler", "Novelty",
"PhaseDeviation", "SpectralFlux"}

I can suppose some of the ranges, e.g. the one of Spectral Centroid, are similar to the range of the audible spectrum, but I have no clue about the more exotic ones.

Is there any trusted source I can rely on?

Unfortunately, the original standard is not available for free. Furthermore, freely available information on the MPEG-7 descriptors seems to be scarce. Therefore, my recommendation is to get [1] which contains in-depth information about all descriptors, including their computation. Once you understand what each descriptor is supposed to describe, and how it is computed, you'll get an idea about its range of values.

The following resources may also be of interest:

An online computation of a number of MPEG-7 audio descriptors for a wave file

References

[1] H.-G. Kim, N. Moreau, and T. Sikora, MPEG-7 audio and beyond: audio content indexing and retrieval. Chichester, West Sussex, England; Hoboken, NJ, USA: J. Wiley, 2005.


Why the bandwidth of the system with the minimum-order observer is higher than the full-order observer? In Ogata's Modern Control Engineering (5th Ed., page 786) it says:

  The bandwidth of the system with the minimum-order observer is higher
  than that of the system with the full-order observer, provided that
  the multiple observer poles are placed at the same place for both
  observers.

Can anyone explain me why this is so?

I think the key of that statement is the sentence

  [...] provided that the multiple observer poles are placed at the same place for both observers.

Suppose that you have a Bode plot, with poles at $\omega_1,\omega_2,...,\omega_n$. If you leave those poles as they are (i.e. you don't change their locations) and you add more poles to the system, then the bandwidth will decrease.

Why? Well, suppose we add a single pole at $\omega_{n+1}$. In that case, the slope of the magnitude as a function of the frequency will decrease for all the frequencies such that $\omega>\omega_{n+1}$. If the magnitude was increasing, it will increase slower or will become constant  (referring to the asymptotic Bode plot). If it was decreasing, it will decrease faster.

If you add even more poles and not just one, then, following the same reasoning, the bandwidth will decrease more abruptly.


Kernel vs impulse response: what's the difference? I am reading The Scientist and Engineer's Guide to
Digital Signal Processing by S.W. Smith. In Chapter 9: Applications of the DFT in the 3rd paragraph he writes "kernel (impulse response)"
What does it mean? What is meant by kernel? Is it same as impulse response?

A kernel is a more general concept, but an impulse response is a special case of a kernel. One usage of the term kernel is to describe an integral transform:
$$y(t)=\int_{-\infty}^{\infty}x(\tau)K(\tau,t)d\tau\tag{1}$$
The function $K(\tau,t)$ is called the kernel of the integral transform.
If you compare $(1)$ to the convolution
$$y(t)=\int_{-\infty}^{\infty}x(\tau)h(t-\tau)d\tau\tag{2}$$
then you see that the kernel of the convolution $(2)$ is given by
$$K(\tau,t)=h(t-\tau)\tag{3}$$
Note that the term kernel is also used to describe linear transformations in the discrete domain:
$$y[n]=\sum_{m=-\infty}^{\infty}x[m]K[m,n]\tag{4}$$
Discrete-time convolution is again a special case of $(4)$ with $K[m,n]=h[n-m]$.
In signal processing we regularly use linear transforms such as $(1)$ and $(4)$. The most well-known examples apart from convolution are the Fourier transform, the Laplace transform, the $\mathcal{Z}$-transform, and the Hilbert transform (which is in fact just a convolution).


Twisted convolution - simplifying multiple integral via substitution I need help with understanding a proof from a paper. Following is the description:

Let
$$r_{1}(t) = \int \int h_{1}\left(\tau, \nu\right) e^{j 2\pi \nu (t - \tau)} s\left(t - \tau\right) d\tau d\nu$$
$$r(t) = \int \int h_{2}\left(\tau, \nu\right) e^{j 2\pi \nu (t - \tau)} r_{1}\left(t - \tau\right) d\tau d\nu$$
Substituting $r_{1}(t)$ in $r_{2}(t)$  we obtain after some algebraic
manipulations,
$$r(t) = \int \int h\left(\tau, \nu\right) e^{j 2\pi \nu (t - \tau)} s\left(t - \tau\right) d\tau d\nu$$
where, $f(\tau, \nu)$ is given by
$$h(\tau, \nu) = \int \int h_{2}\left(\tau', \nu' \right) h_{1}\left(\tau - \tau', \nu-\nu' \right) e^{j 2\pi \nu' \left(\tau - \tau' \right)} d\tau' d\nu'$$

I am unable to use appropriate algebraic manipulations to arrive at the final expression for $r(t)$. Can anyone help me? The first step is arrived at by substituting $r_{1}(t-\tau)$ in $r_{2}(t)$. But I am unable to proceed further. (Incorporated corrections from Tanya Choudhary's answer)
$$r(t) = \int \int h_2(\tau,\nu)e^{j 2\pi\nu(t-\tau)} \int \int h_1(\tau_1,\nu_1)e^{j 2\pi\nu(t-\tau-\tau_1)}s(t-\tau-\tau_1)d\tau_1 d\nu_1 d\tau d\nu$$
Substituting $\tau = \tau', \tau_1 = \tau - \tau', \nu = \nu', \nu_1 = \nu - \nu'$,
$$r(t) = \int \int h_2(\tau',\nu')e^{j 2\pi\nu'(t-\tau')} \int \int h_1(\tau - \tau',\nu - \nu')e^{j 2\pi\nu(t-\tau)}s(t-\tau)d\tau' d\nu' d\tau d\nu$$

Since these are two different integral equations, during substitution you cannot use the same variables for both the integrals.
On substituting $r_1(t-\tau)$ to $r(t)$ we get,
$$r(t) = \int \int h_2(\tau,\nu)e^{j 2\pi\nu(t-\tau)} \int \int h_1(\tau_1,\nu_1)e^{j 2\pi\nu(t-\tau-\tau_1)}s(t-\tau-\tau_1)d\tau_1 d\nu_1 d\tau d\nu$$
Now substitute $\tau = \tau', \tau_1 = \tau - \tau', \nu = \nu', \nu_1 = \nu - \nu'$


What kind of edge detection use paint.net? What kind of edge detection is used Paint.NET program? It allows me to choose the angle of the edge not like the one used by opencv (Canny).

Here an example:

The edges doesn't have single pixel width, so it's probably not canny (no non-maxima supression). Could be a simple high-pass filter, whose pass-band is masked out in the Fourier domain in specific angles.


Why non gaussian variables are independent I have read a lot about ICA. I think I could not find the answer that why non-Gaussian Variables are independent.  What I understand, Central Limit Theorem states that

  distribution of sum of independent variable tends toward more Gaussian than its original random variables.

$
s_i
$ are the number of original independent sources in ICA; whereas ICA model is $x=As$. So we can define $y=w^Ts$. The main goal is to find the unmixing matrix $w$ that maximize the non-Gausaanity. So my question is what is non-Gaussanity here and why its necessary to maximize it to extract the original sources.

Any Enlightenment Please!

The model ICA uses says that there exist some unknown, statistically independent sources, $s_i$ that are non-normally distributed (their distributions are something other than Gaussian):

$$
s_i \sim S(\mu_{s_i}, \sigma^2_{s_i})
$$

where $S$ is some (possibly) known but non-Gaussian distribution with mean $\mu_{s_i}$ and variance $\sigma^2_{s_i}$.

Then it is assumed that what you can actually measure is the addition of these:
$$
\mathbf{x} = \mathbf{A}\mathbf{s}
$$
where $\mathbf{s}$ is the vector of the $s_i, i=1,\ldots,N$, $\mathbf{A}$ is an $N\times N$ matrix (usually assumed to be invertible) and $\mathbf{x}$ is the vector of the actual measurements.

Because the $x_i$ will just be weighted sums of the $s_i$:
$$
x_i = a_{i1} s_1 + a_{i2} s_2 + \cdots + a_{iN} s_N
$$
the central limit theorem says that the distribution of $x_i$ will be closer to Gaussian than the distribution of the $s_i$ (provided certain non-restrictive conditions on the true distribution of the $s_i$ are  met).

Then all ICA really tries to do is to find a $\mathbf{W} = \mathbf{A}^{-1}$ the inverse of $\mathbf{A}$.

There are any number of measures for "non-Gaussianity".  One of the simplest (?) is to use kurtosis as the measure of how far the sample of a random variable is from Gaussianity.

So, to attempt to answer:

  So my question is what is non-Gaussianity here and why its necessary to maximize it to extract the original sources.

There are several different ways one could measure non-Gaussianity. For example, kurtosis can be used as it is known that the kurtosis of a Gaussian is 3. Any distribution with a kurtosis different from 3 is therefore "non-Gaussian" to some extent.

The reason we want non-Gaussianity is because we assumed that the original sources are non-Gaussian.


Solving a Convolution Problem of a 1D Signal I'm finding in trouble trying to resolve this exercise. I have to calculate the convolution of this signal:

$$y(t)=e^{-kt}u(t)\frac{\sin\left(\dfrac{{\pi}t}{10}\right)}{({\pi}t)} $$

where $u(t)$ is Heavyside function

well I applied the formula that says that the convolution of these two signal is equal to

$$Y(f)=X(f)\cdot W(f)$$

where $X(f)$ is the Fourier transform of the first signal and $W(f)$ is the Fourier transform of the second signal

well Fourier transform of $e^{-kt}u(t)$ is $$X(f)=\dfrac{1}{k+j2{\pi}f}$$

I have to make second signal as equals as possible to $\text{sinc}\left(\dfrac{t}{10}\right)$

so I do this operation:

$$\dfrac{\sin\left(\dfrac{{\pi}t}{10}\right)}{\left(\dfrac{{\pi}t}{10}\right)}{\left(\dfrac{1}{10}\right)}$$ this is equal $${\left(\dfrac{1}{10}\right)}\text{sinc}\left(\dfrac{t}{10}\right)$$

right or not?

Even though I realize that this is a very late response, I will nevertheless try to answer this question because I find it instructive and also because the number of upvotes suggests that this question is of general interest to the community.

As already suggested in the question, let's define two signals $x(t)$ and $w(t)$ as
$$x(t)=e^{-kt}u(t),\quad k>0\\w(t)=\frac{\sin(\pi t/10)}{\pi t}$$

One possible interpretation of the convolution $(x*w)(t)$ is that an exponentially damped signal $x(t)$ is filtered by an ideal lowpass filter with impulse response $w(t)$. In the question it was also correctly pointed out that convolution in the time domain corresponds to multiplication in the frequency domain. The Fourier integral of $x(t)$ can be easily computed:

$$X(j\omega)=\int_{0}^{\infty}e^{-kt}e^{-j\omega t}dt = \frac{1}{k+j\omega}$$

The Fourier transform of $w(t)$ should be familiar because it is an ideal lowpass filter. In the question there was some confusion concerning the definition of the Sinc function. I suggest to simply remember the impulse response of a unity gain lowpass filter with cut-off frequency $\omega_0=2\pi f_0$ without making use of any of the definitions of the Sinc function:
$$h_{LP}(t)=\frac{\sin \omega_0 t}{\pi t}\tag{1}$$

Comparing (1) with the definition of $w(t)$, we see that $w(t)$ is simply a unity gain lowpass filter with cut-off frequency $\omega_0=\pi/10$:
$$W(j\omega)=u(\omega+\omega_0)-u(\omega-\omega_0)$$ where I've used the step function $u(\omega)$ in the frequency domain.

In order to find the time function $y(t)=(x*w)(t)$ one can compute the inverse Fourier transformation of $Y(j\omega)=X(j\omega)W(j\omega)$:

$$y(t)=\frac{1}{2\pi}\int_{-\infty}^{\infty}X(j\omega)W(j\omega)e^{j\omega t}d\omega=
\frac{1}{2\pi}\int_{-\omega_0}^{\omega_0}\frac{1}{k+j\omega}e^{j\omega t}d\omega$$

Unfortunately, there is no closed form solution of this integral using elementary functions. It can be evaluated numerically using the exponential integral $\text{Ei}(x)$, or, alternatively, the sine and cosine integrals $\text{Si}(x)$ and $\text{Ci}(x)$. So I don't think that the purpose of the exercise was to actually compute the convolution, but its purpose was probably to come up with a qualitative description of what is going on (exponential signal filtered by an ideal lowpass filter).

Nevertheless, I thought it would be instructive to take a look at the signal $y(t)$, so I evaluated it numerically for the parameters $k=0.05$ and $\omega_0=\pi/10$. The following figure shows the result:

The green curve is the input signal $x(t)$ and the blue curve is the filtered signal $y(t)$. Note the (non-causal) ripples of $y(t)$ for $t<0$ caused by the ideal (non-causal) lowpass filter. If we increase the cut-off frequency of the lowpass filter, the distortion of the input signal should become smaller. This is shown in the next figure where I've increased the cut-off frequency by a factor of 10, i.e. $\omega_0=\pi$ (instead of $\pi/10$):


Determining Type and Bandwidth of a filter Given a filter, if it is given as an equation such as:

$$f(x,y) = \left(\frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2}\right) \exp\left(-\frac{x^2 + y^2}{\sigma^2}\right)$$

Or in a kernel such as:
$$
\left[\begin{array}{rrr}
0 & 1 & 0 \\
1 & -4 & 1 \\
0 & 1 & 0
\end{array}\right]
$$
How much can I find out about the filter. Specifically can I find out if it is a highpass/ lowpass/ bandpass filter (if so how?), and can I find out the bandwidth of the filter?

If it helps, the context is image processing.

Thanks

As others have mentioned, performing a 2D FFT on the kernel will give you the frequency response of the filter. However, it's worth mentioning that 2D filters can be analyzed using the Z-transform, which may or may not provide deeper insight, depending on the filter (and what you want to know).

For example, given the kernel you specified, the corresponding difference equation would be

$$y(n_1,n_2) = x(n_1+1,n_2) + x(n_1,n_2 + 1) + x(n_1-1,n_2) + x(n_1,n_2-1) - 4x(n_1,n_2).$$

Its Z-transform is

$$ Y(z_1,z_2) = z_1X(z_1,z_2) + z_2X(z_1,z_2) + z^{-1}_1X(z_1,z_2) + z_2^{-1}X(z_1,z_2) - 4X(z_1,z_2),$$

which, after rearranging, yields the following transfer function for the filter:

$$ H(z_1,z_2) = \frac{Y(z_1,z_2)}{X(z_1,z_2)} = z_1 + z_2 + z_1^{-1} + z_2^{-1} - 4.$$

To determine the magnitude response, just plug in a pair of complex exponentials and simplify, as follows:

$$ \begin{aligned}
H(e^{iw_1},e^{iw_2}) & = e^{iw_1} + e^{iw_2} + e^{-iw_1} + e^{-iw_2} - 4  \\
                     & = (e^{iw_1} + e^{-iw_1}) + (e^{iw_2} + e^{-iw_2}) - 4 \\
                     & = 2\cos{w_1} + 2\cos{w_2} - 4 \\
|H(e^{iw_1},e^{iw_2})| & = 2\sqrt{ (\cos{w_1} + \cos{w_2} - 2)^2  } .
\end{aligned}
$$

Evaluating the magnitude response at extreme frequencies will give you a sense of highpass vs. lowpass for the filter. For example,

$$ |H(e^{i0},e^{i0})| = 2\sqrt{(\cos{0} + \cos{0} - 2)^2} = 2\sqrt{(1+1-2)^2} = 0,$$
and
$$ |H(e^{i\pi},e^{i\pi})| = 2\sqrt{(\cos{\pi} + \cos{\pi} - 2)^2} = 2\sqrt{(-1-1-2)^2} = 8.$$

Of course, the transfer function can be evaluated directly to plot the magnitude response of the filter as well. Here's an example using numpy:

import numpy as np
import pylab as py
from mpl_toolkits.mplot3d import axes3d

def H(z1,z2):
    return z1 + z2 + 1./z1 + 1./z2 - 4.0

n = 100
w1 = w2 = np.linspace(-np.pi,np.pi, n)
mag = np.zeros((n,n))

for i1 in xrange(0,n):
    for i2 in xrange(0,n):
        z1 = np.exp(1j*w1[i1])
        z2 = np.exp(1j*w2[i2])
        mag[i1,i2] = np.abs(H(z1,z2))

fig = py.figure()
ax = fig.add_subplot(111, projection='3d')
X, Y = np.meshgrid(w1,w2)
ax.plot_surface(X, Y, mag, cmap='bone', alpha=.5)
py.show()

Keep in mind that if you use the kernel 2D FFT technique, the resulting magnitudes won't necessarily be centered at zero like the above plot.


Prove that time-reverse yields inverse filter for deconvolution? In reading literature on the construction of impulse responses from sine sweeps (e.g. papers by Farina), I see it stated over & over that the way the impulse response is constructed from a (exponential) sine sweep is to take the original sine sweep signal $S(t)$, time reverse it ( $S(T-t)$ ), scale it somewhat in magnitude,  and convolve this with the recorded signal to obtain the impulse response of the (time-invariant, linear) system such as a room's reverberation..  i.e., the claim is that $S^{-1}(t)$ is proportional to $S(T-t)$.

Can anyone explain how to prove this, or point me to a reference that actually shows a 'proof' of this?  The closest I've seen is "it can be shown that...", and a few hand-waving comments about how time-reversal inverts the phase -- but it's not obvious that this 'does the job' for producing a proper inverse from a function with nonlinear time dependence.

I've been trying to show this explicitly (in the continuous-time case) for about 7 hours now, and have a few pages of integrals written out, even with a little help from WolframAlpha :-(.   The exponential term in the frequency makes it hard to compute the Fourier Transform explicitly (you get Gamma functions with two arguments), let alone try to take the reciprocal of the FT and then take the Inverse FT of that.

Thanks!

This is not a homework problem, I just want to understand why it works.  (I know it works because I've coded it up, but it still seems kind of like magic.)

PS- Thought I'd try a linear sine sweep as it might be easier, but you also get a frequency-squared term which again makes it...not at all obvious...that taking $t$ -> $-t$ is going to invert all the phase terms.

P.P.S.- Here are my notes so far; the calculation starts on page 4.  Dropbox link to PDF

This is not a general result. It specifically holds, iff the Fourier transform of the signal is of unit magnitude. In that case you can write the signal in frequency domain as

$$S(\omega)=\exp(i \phi(\omega))$$

for a real function $\phi$. The time reversed signal (I'll ignore the constant time offset, it doesn't change anything relevant) is then the conjugate

$$\bar{S}(\omega)=\exp(-i \phi(\omega))$$

Convolving the two by multiplying the frequency domain representation gives

$$\exp(i\phi(\omega))\times \exp(-i\phi(\omega)) = 1$$

Which is the response of the identity system and demonstrates that the signal and its time reversed are the responses of inverses.

This result is only exact in this special case. Sweeps can belong to this class of signal, for example if they have been generated as an allpass filter response. You can also do that by constructing the sweep entirely in frequency domain.

A sweep that is synthesised in time domain using a modulated sinusoid is typically not of that form, but approximates it well. Especially near the endpoints of the frequency range you will get significant deviation from the unit response.


What Does Normalizing Gradient Means? Since my algorithm needs to detect long, coherent edges, the estimated local magnitude P for each pixel is normalized with respect to the magnitudes in it’s local neighborhood W (i.e., a 5x5 window).

$P_{normalized} = \dfrac{P - \mu _w}{\sigma_w}$

where $\mu _w$ and $\sigma_w$ denote the average and the standard deviation of edge magnitudes in the pixel neighborhood $w$.

This is given in the paper "GradientShop: A Gradient-Domain Optimization Framework for Image and Video Filtering" by Bhat, and I need to implement the algorithm described there. I am not getting, what does normalizing the gradient mean exactly ?

Since you're working local it is suggested for you to compare things normalized to their relative surroundings.

The gradient is a vector (2D vector in single channel image).
You can normalize it according to the norm of the gradients surrounding this pixel.

This is what suggested above.


About the proof of an equality related to the DFT [sampling the DTFT to obtain the DFT] This wiki page about the DTFT says that the DFT can be obtained from the DTFT by sampling the latter in one cycle at $N$ points:

  When the DTFT is continuous, a common practice is to compute an arbitrary number of samples (${{N}}$) of one cycle of the periodic function (${X_{1/T}}$):\begin{align}
\underbrace{X_{1/T}\left(\frac{k}{NT}\right)}_{X_k} &= \sum_{n=-\infty}^\infty x[n]\cdot e^{-i 2\pi \frac{kn}{N}} \quad \quad k = 0, \dots, N-1 \\
&= \underbrace{\sum_{N} x_{_N}[n]\cdot e^{-i 2\pi \frac{kn}{N}},}_{DFT}\quad \scriptstyle{\text{(sum over any }n\text{-sequence of length }N)}
\end{align}
  where $x_{_N}$ is a periodic summation:
  $$x_{_N}[n]\ \triangleq\ \sum_{m=-\infty}^{\infty} x[n-mN]$$

So we essentially have the following equality
$$\sum_{n=-\infty}^\infty x[n]\cdot e^{-i 2\pi \frac{kn}{N}}= \sum_{n=0}^{N-1} \sum_{m=-\infty}^{\infty} x[n-mN] e^{-i 2\pi \frac{kn}{N}}. $$

I didn't manage to prove it, is this a known identity?

The right hand side is the periodic repetition of the left, tells the last equation.

Starting from Fourier transform $$X(\omega)=\sum_{n=-\infty}^{\infty}x(n)e^{-j\omega n}$$ and $X(\omega)$ is periodic with period $2\pi$. So only the samples in the fundamental frequency range are sufficient.

Now we take $N$ equi-distance samples in the iterval $0\leq\omega\leq 2\pi$, so the spacing will be $\dfrac{2\pi}{N}$

Let $k$ be the sample index and $k=0,1,\dots,N-1$ and we evaluate $\omega$ at $\left(\dfrac{2\pi k}{N}\right)$

$$X\left(\dfrac{2\pi}{N}k\right)=\sum_{n=-\infty}^{\infty}x(n)e^{-j2\pi kn/N}$$

which can be written as

$$X\left(\dfrac{2\pi}{N}k\right)=\sum_{l=-\infty}^{\infty}\sum_{n=lN}^{lN+N-1}x(n)e^{-j2\pi kn/N}$$

Now change the index $n$ to $n-lN$ which results in

$$X\left(\dfrac{2\pi}{N}k\right)=\sum_{n=0}^{N-1}\left[\sum_{l=-\infty}^{\infty}x(n-lN)\right]e^{-j2\pi kn/N}$$

The sequence inside the square brackets is the periodic repetition of $x(n)$ of period $N$

$$x_p(n)=\sum_{l=-\infty}^{\infty}x(n-Nl)$$

Thus $$x_p(n)=\begin{cases}
               x(n),\;\;\;\;0\leq n \leq N-1\\
                     \\
               0,\;\;\;\;\;\;\;\;\;\;\;     \mathrm{elsewhere}\\
            \end{cases}$$


Z transform stability What is the causality & stability status for three cases shown (aso in attached photo) ?

$$H(z) = \frac{z(z-1)}{(z+1)(z+\frac{1}{3})} $$

for three possible regions of convergence as:

a-) |z| > 1

b-) |z| < 1/3

c-) 1/3 < |z| < 1

Am i ok as below?

a) Causal & marginally stable

b) anti-causal & unstable

c) non-causal & marginally stable

Yes, all answers given by you are fine.

Causal: If ROC is outside the outermost pole

Left sided: If ROC is inside the innermost pole

Stable: If ROC contains the unit circle (marginally stable if it touches unit circle)


How to Draw a Line at Non Integer Pixel Coordinates on an Image I have lines (or segments really) I want to draw into a new image I.
Let's define the intensity of the pixel at row = r and col = c as I(r,c).  However, the lines are defined by 2 points whose coordinates are not integers.

For example, one line might be from point A (with row = 2.2 and col = 3.5) and to point B with (row = 10.8 and col = 30.4).

What is the algorithm or method for filling in the intensity values of image I at the integer rows and cols?

A long time ago I wrote an article on drawing lines with the Wu-Line algorithm. This can be used to draw lines from non-integer coordinates. The effect is very nice to see when animated. You can still find it in the Wayback Machine:

Wu Anti-aliased Lines from The good-looking textured light-sourced bouncy fun smart and stretchy page

There's also some pseudo-code to get you started.


Fractional spaced equalizer + timing (clock) recovery I am trying to figure out how timing recovery is used after equalization (or, being more correct, when the equalizer is in the timing recovery loop), having a T/2 fractional spaced equalizer and using the Gardner algorithm as timing error detector.
The usual situation is to have the equalizer output at 1sample/symbol, by having downsampled, or combining the outputs of the even and odd filters. The latter is my case (case (a) in figure). But Gardner algorithm works at 2samples/symbol, so I wonder what I must do with the samples coming out from the even and odd filters. Must I multiplex the outputs of the two filters so I have a 2samples/symbol, as depicted in (b)?

Regards,

The answer is to NOT down-select to one sample per symbol until after using the Gardner Timing recovery since the TED requires 2 samples per symbol. If the equalizer is running at 2 samples per symbol, that is perfect for use with the Gardner; why would you down-select to one sample per symbol after the equalizer? You can downselect after timing recovery since at that point you have determined the precise sample location.

Also note that the fine timing adjustment can also operate on just 2 samples per symbol using polyphase filters! Below I show an example for a CDMA waveform before and after timing recovery using a Gardner TED at 2 samples per symbol with then with a loop filter adjusted timing using a polyphase resampler. There was no analog change in the sampling clock position nor an actual interpolation to a higher sampling rate (the polyphase filter design is based on interpolation but no actual higher sampling rate is used!).


Real and imaginary part of signal after fft I will take two sine waves i) 0.1 AMPLITUDE with no phase shift ii)0.4 amplitude with some phase shift 20 degree (Considering Both same frequency).

In time domain - If i divide signal 1 by signal 2 by taking RMS value,it was nearly 0.25
In frequency domain - I took fft of both the signals and divided signal 1 by signal 2, i got the result in real and imaginary part.

If i take the RMS value of frequency domain output(real and imaginary), i am not able to get the same result as time domain.

So what should i do to get the same result in both the case.

Thanks

Here is one way to do it properly:

import numpy as np

#=======================================================================
def main():

        t = np.arange( 0, 200 )
        N = len( t )

        signal1 = 0.1 * np.sin( .1 * t )
        signal2 = 0.4 * np.sin( .1 * t + .35 )

        rms_time1 = np.sqrt( np.sum( signal1 * signal1 ) / N )
        rms_time2 = np.sqrt( np.sum( signal2 * signal2 ) / N )
        q_time = rms_time1 / rms_time2

        print q_time, rms_time1, rms_time2

        dft1 = np.fft.fft( signal1 )
        dft2 = np.fft.fft( signal2 )

        dft1_real = np.real( dft1 )
        dft1_imag = np.imag( dft1 )
        dft2_real = np.real( dft2 )
        dft2_imag = np.imag( dft2 )

        rms_freq_real1 = np.sqrt( np.sum( dft1_real * dft1_real ) / N )
        rms_freq_imag1 = np.sqrt( np.sum( dft1_imag * dft1_imag ) / N )
        rms_freq_real2 = np.sqrt( np.sum( dft2_real * dft2_real ) / N )
        rms_freq_imag2 = np.sqrt( np.sum( dft2_imag * dft2_imag ) / N )

        rms_freq1 = np.sqrt( rms_freq_real1 * rms_freq_real1 \
                           + rms_freq_imag1 * rms_freq_imag1 )

        rms_freq2 = np.sqrt( rms_freq_real2 * rms_freq_real2 \
                           + rms_freq_imag2 * rms_freq_imag2 )

        q_freq = rms_freq1 / rms_freq2

        print q_freq, rms_freq1, rms_freq2

        print rms_freq1 / rms_time1
        print rms_freq2 / rms_time2
        print np.sqrt( 200 )

#=======================================================================
main()

The output is:

0.246138603222 0.0699023232816 0.283995774603
0.246138603222 0.988568136262 4.016306761
14.1421356237
14.1421356237
14.1421356237

Followup:

This is a more straightforward way to calculate the RMS of the DFT bins.  I coded the sample above to more closely align with how I thought the OP was doing it.

        sumsquares1 = np.real( np.sum( dft1 * dft1.conjugate() ) )
        sumsquares2 = np.real( np.sum( dft2 * dft2.conjugate() ) )

        rms_freq1 = np.sqrt( sumsquares1 / N )
        rms_freq2 = np.sqrt( sumsquares2 / N )

        q_freq = rms_freq1 / rms_freq2

        print q_freq, rms_freq1, rms_freq2


Is R suitable for digital signal processing While asking a question about representing large time series in R I was discouraged from using R for digital signal processing. I understand that R is geared towards statistics. However, a signal is just a set of measurements. R also deals with time series, but they are usually at different time scales. Nevertheless, e.g. when predicting price evolution it is a matter of understanding the model behind it and possibly predicting it. That's quite close to what I usually want in physics.

In "traditional" digital signal processing (DSP), it is more common to analyze the signal in terms of its frequency components and their interplay or coherence (which is something like frequency dependent covariance of the amplitudes of frequency components) between signals.

I'm interested in R mostly because of some packages on CRAN with advanced scientific routines (which aren't as developed in Python) like signal for general DSP, wavelets and biwavelet for wavelet transforms and timsac for bispectrum. These packages should make this more "traditional" approach possible. However, I get the feeling that time series are generally analyzed in a different way in R as is apparent from the time series task view on CRAN. Are these more statistical models perhaps unsuitable for analyzing non-stationary time series?

So I ask whether there is something that makes R worse for DSP than e.g. Python+SciPy (which I currently use the most) or MATLAB. Is it perhaps some underlying implementation that makes operations slower? Or is it just that it would be necessary to change one's approach to analyzing signals to a more statistical one when using R? And would such an approach be perhaps limited in some way in comparison to the "traditional" DSP approach?

Note:

I tried asking this question on SO, but I was directed to rather ask here.

I'm interested in objective answers with sources, not just opinions. E.g.

which internal representations makes it slower
what are the limitations of the statistical models used for time series analysis, e.g. how they cope with non-stationary signals, n-wave interaction, etc.

Since the bulk of R’s DSP capability comes from the signal package which was ported over from the open source project Octave (itself influenced by MATLAB), there's no intrinsic limitation of R.

What you have picked up on, are ecosystem preferences. We learned MATLAB in university, picked up numpy/scipy/sklearn at work, so R isn't the first weapon of choice.

That said, there are some interesting "bridge" projects that you may be interested in exploring: https://stackoverflow.com/questions/5630441/how-do-rpy2-pyrserve-and-pyper-compare.

There is also the excellent interactive python (iPython) project, now called Jupyter, that recently introducing support for R notebooks: https://www.continuum.io/blog/developer/jupyter-and-conda-r. So tools & their inter-operability are converging!


Is system $h(t) = u(t)-u(t-1)$ LTI or not? $$y(t) = \int_{0}^{\infty} x(\tau)h(t-\tau) \ \mathrm{d}\tau$$
where $h(t) = u(t)-u(t-1)$.

I have some problem with understanding whether system LTI or LTV just looking the impulse response of the system. Is there any short-cut to identify the system time-variance looking the impulse function?

If a system is described by an impulse response $h(t)$, then by definition it is an LTI system. Because impulse response of the form $h(t) = \mathcal{T} \{ \delta(t) \} $ only exists for LTI systems.


How do you apply Kalman Filter to track a signal? The example that I've seen on state estimation involves deriving the ABCD matrix of a physical system (i.e. falling object) and tracking that object.

I would like to use Kalman Filter for signal denoising applications (specifically EEG signal). How could I apply KF in this situation since ABCD matrix cannot be derived for a signal.

The power of the Kalman filter lies in the effect that it predicts the next "state" of the signal/object by using an internal model of the process. That why it is very effective for physical processes, because they can (often) be modelled with quite large precision.

The fact that you point out that it is impossible to derive a model for your signal renders the effect of using (and even the possibility of using)  Kalman filter usesless/impossible.


What kind of information can we obtain from curvature for exemplar-based inpainting? I just want to know what kind of information does Curvature suggest to us? I've seen an article that used from the curvature of an Image. I can't understand what's the purpose of the author or it's better to say that what kind of information Curvature has!!!

Can anyone help me to handle this problem?

The link of the article:
Article link

The article defines isophote curvature as:
$$
K(p) = \bigtriangledown \left[ \frac{ \bigtriangledown I_p}{ \left|  \bigtriangledown I_p \right|} \right]
$$

I haven't read the paper before, but the rationale is state in III.A:

  The filling order is very crucial which directly affects the quality of image inpainting. However, during filling procedure we find the data terms may be zeros for many pixels, which make the priorities being close to zeros, although the confidence terms may be very big. This would result in incorrect filling order.[...] by introducing the curvature term,
  we effectively avoid the priority being close to zero while the
  data term is zero.

The foundation for the method used in this paper is in Object removal by exemplar-based
inpainting, so be sure to read that first.


Two-Sided Frequency Spectrum I am trying to make FFT simulation in Matlab by generating noise added
two sinus waves in 60Hz and 100Hz.

After adding the noise into these signals then I have applied the FFT
as I put my Matlab code below.

But I am in difficulty interpreting the FFT plot which shows two spectral peaks in the plot diagram below.

Could you please explain how to interpret this spectrum?
Why we see two peaks in the spectrum? How can we reduce the spectral to show only the frequencies of 60Hz and 100Hz?

%%%% Noise_Added_Two_Sinus%%%
>> f1=60;
>> f2=100;
>> fs=512;
>> t=0:1/fs:2-1/fs;
>> x1=2.4*sin(2*pi*f1*t);
>> x2=0.96*sin(2*pi*f2*t);
>> y=x1+x2+randn(size(t));
>> F=fft(y);
>> plot(abs(F))

First of all, you should totally fix the $x$-axis :)

Now, as you only have two sine waves, you should expect to have peaks at their exact frequencies (that is, at $\pm 60\textrm{ Hz}$ and $\pm 100\textrm{ Hz}$). The FFT function in MATLAB gives you the Discrete Fourier Transform of your signal within $0$ and $F_s$, where $F_s$ is your sampling frequency -and given that you represent your signal with samples, even if you are not aware of it, you have one. Anyway, to make visualization easier, we focus on the interval $[-\frac{F_s}{2},\frac{F_s}{2}]$; that way, your spectrum will be centered around the $0$ frequency -use the function fftshift for that. So, now, once you fix the $x$-axis according to your sampling frequency, you will find that the peaks are at the frequencies you have defined for your sine waves.

When it comes to the little ripples in the lower part of your spectrum, that's just disturbance caused by the noise. Since you are adding noise to the signal by the randn function, your noise follows a Gaussian distribution, translating into white noise for your signal. That is why your entire spectrum is showing sings of being slightly disturbed.

By the way, after fftshift, just get rid of the first half of your vector and you'll end up with positive frequencies.


Equation of wave in Linear frequency modulation? We know that we have the following equation for wave.
$$g(t)=A\cos(\omega t+\theta_0)=A\cos(2\pi ft+\theta_0)$$
The equation of frequency with respect to time will be:

$$f(t)=\frac{Bt}{\tau}+f_c-\frac{B}{2}$$
Then:

\begin{align}
g(t)&=A\cos(2\pi (\frac{Bt}{\tau}+f_c-\frac{B}{2})t+\theta_0)\\
&=A\cos(2\pi(f_c-B/2)t+\mathbf{2\pi}(B/\tau)t^2+\theta_0)
\end{align}

Then why that bold 2 is omitted in this answer?
$$f(t)=A\cos(\theta(t))=A\cos(2\pi(f_c-B/2)t+\pi(B/\tau)t^2+\theta_0)$$

If you have a signal

$$g(t)=\cos(2\pi \hat{f}(t)t)\tag{1}$$

then the function $\hat{f}(t)$ is not the instantaneous frequency of $g(t)$ (unless $\hat{f}(t)$ is constant).

If you want an instantaneous frequency $f(t)$, then the equation

$$\frac{\phi'(t)}{2\pi}=f(t)\tag{2}$$

must be satisfied, where $\phi(t)$ is the phase of the signal $g(t)$. So in order to obtain the phase $\phi(t)$, you have to integrate the desired instantaneous frequency $f(t)$. For

$$f(t)=\frac{Bt}{\tau}+f_c-\frac{B}{2}\tag{3}$$

you get

$$\frac{\phi(t)}{2\pi}=\frac{B}{2\tau}t^2+(f_c-\frac{B}{2})t+\frac{\theta_0}{2\pi}\tag{4}$$

So the signal with the desired instantaneous frequency is

$$g(t)=\cos(\phi(t))=\cos\left[\pi(Bt^2/\tau+(2f_c-B)t)+\theta_0\right]\tag{5}$$

Also take a look at this related question and its answer.


What does the frequency axis of a Power Spectral Density mean? I have never really understood what the frequency axis meant when we plot the Power Spectral Density(PSD).

Does it correspond to frequency as we get after we take the Fourier Transform of a time domain signal. This doesn't make sense to me because the frequency here(when we do FT of a time domain signal) corresponds to the time variable "t".

To find the PSD we take the Fourier transform of the Autocorrelation function which is a function of either time difference $$\tau = t_2 - t_1$$ i.e. $$R_{xx}(\tau) $$ (for the Wide Sense Stationary case) or $$R_{xx}(t_1,t_2)$$ for other cases.

Here the autocorrelation is a function of either time difference or two time variables. If this is the case then how can we say that the PSD gives the power distributed in the frequency domain?

How do I understand this. If someone could give both a mathematical explanation and an intuitive explanation I would be really grateful.

Thank you.

The fact that the frequency variable of a power spectral density (PSD) equals the one of a Fourier transform of a "normal" time-domain signal can be seen more easily by considering the following definition of the PSD of a wide-sense stationary (WSS) random process $x(t)$:
$$S_x(\omega)=\lim_{T\rightarrow\infty}E\left\{ \frac{1}{T}\left| \int_{-T/2}^{T/2}x(t)e^{-j\omega t}dt \right|^2 \right\}\tag{1}$$
Note that $(1)$ is just the limit of the expectation of the scaled squared magnitude of the Fourier transform of a truncated version of $x(t)$. So the frequency variable $\omega$ is just the frequency variable of the (squared magnitude of the) Fourier transform.
The Wiener-Khinchin theorem says that the autocorrelation function and the PSD of a WSS random process form a Fourier transform pair.
This answer shows that the inverse Fourier transform of $(1)$ indeed equals the autocorrelation function of $x(t)$.


Computing narrowband spectrograms using MATLAB I'm having an assignment of computing narrowband spectrograms using MATLAB.
And I totally have no idea about the code.
Can someone help me writing the code and explaining what's going on to me.
Lots of thanks!!

Followings are the requirement:

B.1. Load the recorded .wav file using the function “wavread”. Let y be the variable
storing the signal;

B.2. Perform short-time Fourier transform on y as follows:
(i) Divide the signal into many short-time frames. Each frame is 50 msec long. Frame
shift is 1 msec;

(ii) Apply 1024-point FFT on each frame;

(iii) Store the power spectrum of each frame in a column vector (note that only the first
512 FFT coefficients are needed).

(iv) Append the power spectrum vectors of all frames to form super matrix S. The size
of the matrix will be 512 × no. of frames;
(v) Plot the spectrogram using the following MATLAB commands

No one is going to give you a full working code - the point is that you should learn and understand how to do it. Following guidelines for you:

As given, use the wavread function to load the recording (please mind that in the newest 2014a version this function is deprecated, you should use audioread instead). MATLAB has wonderful help with very good examples, so simply use it.
Having your signal loaded, together with its sampling frequency, you need to split it into chunks (many short-time frames). This can be done in two ways, either you are using a loop construct (for, while - whichever you prefer), or (my favourite) use the MATLAB buffer function (please consider using 'nodelay' option in it, as first frame will have zeros in it up to overlap length). Although you require the Signal Processing Toolbox license (I'm sure your University has it). This function allows you to specify the length of frame, overlapping and other things. In result you will obtain the matrix which columns are your signal frames. Having the sampling frequency from wavread, it is very easy to calculate corresponding length of your frame in samples for given number of miliseconds
Now you can iterate over the columns of your matrix and calculate the FFT for each of them. You might want to use the fft function to get it. In extra parameter you can specify the length of your signal (to apply zero padding eventually). So let's say that now you have FFT for first frame. Now take the magnitude of it (abs), and choose first 512 samples. Store this result into fist column of your super matrix $S$. Optional step: before calculating the FFT, I suggest you to multiply your frame by the window function, i.e. hamming window.
Now you have the matrix $S$, so plotting can be done by using various functions (pcolor, imagesc, and so on - pick whichever you prefer).

I suggest you to compare the final result with MATLAB built-in spectrogram function.

Good luck!


In Kalman filters why is it necessary to transform the systems dynamics matrix to the state transition matrix? Previously when I have implemented Kalman filters I have used the transformation

$$
\mathbf{A(t)} = \mathcal{L}^{-1} \left( s \mathbf{I} - \mathbf{F} \right) ^{-1}
$$

to calculate the state transition matrix $\bf{A}$ from the system dynamics matrix $\bf{F}$.

Where the equation describing the state of the system is

$$
\vec{x_t} = \mathbf{F_t} \vec{x_{t-1}} + \mathbf{B_t} \vec{u_t} + \vec{w_t}
$$

and my prediction step calculation is

$$
\vec{x_{k|k-1}} = \mathbf{A} \vec{x_{k-1|k-1}} + B\vec{u_k}
$$

For example for simple sinusoidal Kalman filter started from the following equation of motion:

$ \ddot{x} = -\omega^2x $

So my matrix equation of motion was:

$$
\begin{bmatrix}
  \dot{x} \\
  \ddot{x} \\
\end{bmatrix}
=
\begin{bmatrix}
  0 & 1 \\
  -\omega^2 & 0 \\
\end{bmatrix}
.
\begin{bmatrix}
  x \\
  \dot{x} \\
\end{bmatrix}
$$

and my system dynamics matrix was

$$
F = \begin{bmatrix}
  0 & 1 \\
  -\omega^2 & 0 \\
\end{bmatrix}
$$

I then performed the transformation $
\mathbf{A(t)} = \mathcal{L}^{-1} \left( s \mathbf{I} - \mathbf{F} \right) ^{-1}
$ to get the state transition matrix like so:

$$
A(t) = \begin{bmatrix}
  \cos(\omega t) & \dfrac{\sin(\omega t)}{\omega} \\
  -\omega \sin(\omega t) & \cos(\omega t) \\
\end{bmatrix}
$$

Where I then replaced $t$ with my sample time $T_s$ where $T_s = \dfrac{1}{F_s}$ where $F_s$ was my sample frequency.

I then used this to extract my sinusoidal signal.

My question is did I need to perform this transformation? And if I did, why did I have to do so?

I was reading the paper linked here where they attempt the explain the Kalman filter in simple way and they did no such transformation.

As I said in the comments, there should be no need at all to do such a transformation if all of your equations start in discrete-time.

The equation
$$
 (s \mathbf{I} - \mathbf{F} )^{-1}
$$
is just the Laplace transform of the solution for the state equation:
$$
\dot{\vec{x_t}} = \mathbf{F} \vec{x_t} + \vec{w_t}
$$

However, I really can't see how $\mathbf{F}$ works in:
$$
\vec{x_t} = \mathbf{F_t} \vec{x_{t-1}} + \mathbf{B_t} \vec{u_t} + \vec{w_t}
$$
because this is a discrete-time equation, and the equation you're looking at generating the sinewave is continuous-time.

Can you clarify?


Help Solving a Difference Equation I have an exam in my signals and systems class in a couple of days, and I'm unsure how to go about solving this practice problem. The question is as follows:

Consider a discrete time system whose input and output are related by the following difference equation.

$$  y[n] =
\begin{cases}
(1/4)y[n-1] + x[n],  & \text{if $n>=0$} \\
(1/4)y[n-1] + x[n] + x[n+1]*\cos(n*u[n] + \pi/2), & \text{otherwise}
\end{cases}
$$

Use recursion to find an expression for the zero-state output, y[n], when the input is given by
$$x[n] = δ[n + m]$$
where m is an integer

I'm having a bit of trouble with the fact that the equation extends to before 0; usually to solve one of these by recursion we'd have something like y[-1] = 0, which would allow us to go ahead and compute more ys ... but in this case that seems to be impossible.

The only thing I managed to notice was that cos expression is irrelevant because of the unit step function; it's all 0 before 0, and that whole term goes away.Also recursively the y[n-1] term seems to divide by infinity ... but I don't know where to take it from there.

You seem to have figured out that the system is actually given by

$$y[n]=\frac14 y[n-1]+x[n],\quad \forall n\tag{1}$$

Since you're looking for the zero-state output, you can simply compute the system's impulse response $h[n]$ by applying $x[n]=\delta[n]$ as an input signal (with $y[n-1]=0$), and compute the response to any other input signal using this impulse response. For $x[n]=\delta[n+m]$, the response is simply $h[n+m]$.

The recursion for the impulse response is

$$\begin{align}h[0]&=1\\
h[1]&=\frac14\\
h[2]&=\left(\frac14\right)^2\\
&\vdots\end{align}\\$$

I'm sure you can take it from here and derive a general expression for $h[n]$.


Which one can be accompanied by linear filters? I have a matrix :
\begin{bmatrix}
1 & 2 & 3\\
1 & 4 & 5\\
2 & 6 & 7
\end{bmatrix}
After doing operation 1 ,
I get
\begin{bmatrix}
0 & 0 & 0\\
1 & 2 & 3\\
1 & 4 & 5
\end{bmatrix}
However operation 2 on same matrix will give :
\begin{bmatrix}
0 & 1 & 2\\
0 & 3 & 4\\
1 & 5 & 6
\end{bmatrix}
Clearly, operation 2 reduces each matrix value by 1, while operation 1 shifts one row down.
Now I see that operation 1 is simply subtraction and thus can be applied using linear filter. But I am not using to linear filters on matrix. Can anyone guide me through this?

TL;DR: Time shift in doable with LTI systems. Subtracting $1$ is non linear.
In discrete systems, integer shifts can be implemented with filters that are zero everywhere, except at the shift location.  This can be implemented via circulant matrices, akin to convolution. This can be made compatible with Operation 1 (taking border effects into accounts).
The second operation adds one to every value. This is not linear (but affine), and therefore does not match with linear filters.


Conversion of FFT to PSD I would like to convert FFT data which is set of 396 values for 396 frequencies (from 4Hz to 399Hz by 1Hz) in unit of g's (gravitational akcceleration) to PSD in unit of (g^2)/Hz, is it possible? Thank you for your help.

If X is your fft data then you need to take Y=abs(X).^2; to convert it. Your frequency axis is f=4:399;. You can plot the result then using plot(f,Y).


Changing voice characteristics so it sounds like someone else voice I have a recording of speech, and I want to manipulate the voice, so a listener
will perceive the manipulated voice and the original one as belonging to two
different people. The manipulated voice should sound very natural, not computer-generated.

Just to clarify, the manipulated voice does not have to be modeled after the voice of someone; I'm not asking about mimicking a voice, as discussed here.

I have already used the SoundTouch library for pitch shifting, and it produces very nice results, but I want to change not only the pitch.

I saw this impressive video of a software called Voice Sculptor (couldn't find the company site), which demonstrate such a manipulation.

So I'm asking on a conceptual/algorithmic level how this can be done
And also if you know how to achieve this with specific tools.

Last note: I'm new to DSP, so I will appreciate any explanation of why and how.

You may use the same pitch shift technique, but a little bit modified. If, classically, in pitch shift algorithm, each frequency component changes by multiplication of its frequency on some pitch-shift coefficient k, you may try:

Changing a little bit the amplitude of each spectral component.
Multiply each frequency component not by k, but k+v(f), where v(f) shows variation dependent on frequency f. Such variation must be of much smaller value than k. I mean k >> v(f).