What is spectrum permutation? This paper refers to the concept of "spectrum perumutation" on page 4. I'm having a difficult time understanding what it means to permuate a signal. This is mainly due to the fact I don't understand the "invertible mod n" means.

What exactly happens to a signal when we permutate it? What does it look like in the time domain and what does it look like in the frequency domain?

Reminder what a permutation is: The exchange of indices in a sequence.

That's exactly what they're doing here.

Invertible in this context (only read the paragraph so far) means

"the multiplication with $\sigma$ under mod $n$ is invertible, i.e. given the result of that multiplication and $\sigma$, you can caclulate the original value in all cases."

Let's do an example:

let's say we're  $x\in \mathbb Z
\setminus 13$, and you multiply by $\sigma = 7$, you get $x\sigma = 4 \mod 13$. You know that $x=8$, undisputably! (because 7·8  = 56, and 56 mod 13 = 4, and it's only that easy because this is a prime field).

That works for any possible result $xz\in \mathbb Z\setminus 13$, because $\sigma$ is coprime to 13; so, $\sigma$ is invertible for $(\cdot,\mathbb Z \setminus 13)$.

I picked that 13 because it's prime, and prime numbers are by definition coprime to any smaller integer.

Let's do an counter-example of that:

For FFTs, the length is usually not a prime number, but can be any natural integer. $n=2^m, m \in \mathbb N$ is a very common FFT length (because you can decompose it easily using simple butterflies); then, not every $n$ is invertible in the sense above: for example, if your multiplication $x\cdot4 = 0 \mod 32$, you can't say whether $x=0$, $x=8$ or $x=16$; it's not invertible.


Why do we assume zero mean noise in sensor data? I am reading a paper on measuring respiratory patterns from video data. In defining the model, the authors formulate the problem mathematically as:

  $x_i(t)=h_i(t) \ast g(t) + n_i(t) $

Where $n_i(t)$ is a zero mean noise. My question is why do we model the pixel data as having zero mean noise? What prevents it from having a mean?

The brutally honest answer here is: The noise is considered zero-mean because that's what the author decided to do. Without looking deeper into the signal model employed, it's impossible to answer.

However, for many systems this makes a lot of sense physically, since the processes leading to a noise realization are very often zero-mean in nature. For example, the thermal noise over any resistor should be zero-mean, simply because, well, that's how these random fluctuations work.


Re-implementing the "spectrogram" function from matlab I am trying to make a spectrogram viewer without using the spectrogram command. For this purpose, I used a sin function. I broke up the input signal into 256 segments, multiplied each segment with the Hamming window, took the FFT of each slice and then stored them to obtain spectra of the signal.

My problem is that I couldn't get the exact spectrum. It is probably due to the fact that I could not divide input signal into overlapping segments?

This is my code: Can someone help me where I am doing it wrong?

Fs = 1000;                    % Sampling frequency
T = 1/Fs;                     % Sample time
L = 10240;                     % Length of signal
NFFT=10240;                   % number of fft points
t = (0:L-1)*T;                % Time vector
x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);
X = hann(L).*x';
figure
plot(X)
k=256;

for idx = 1:length(x)/k
     %take a slice from the input vector
    slice =  x(1+(idx-1)*k:idx*k);
    %multiply it with the window and transform it into frequency domain
    spectrum =  fft(slice'.*hamming(k),k);
    %get the spectrum magnitude at each of the 256 frequency points and store it
    mag_spectrum(:,idx*1024) = abs(spectrum).^2;

end

figure
w=hamming(k);
plot(w)
figure
plot(t,X,'k');
xlabel('Time (sec)','FontSize',14);
ylabel('\it{x(t)}','FontSize',14);

% spectrogram(x,128);
% SET PARAMETERS

figure
plot(slice); % To look for it is windowed or not which is not in my case?

%figure
%plot(mag_spectrum);
figure
spectrogram(X,256); % Real spectogram

Differences between your implementation and what the spectogram function is displaying.

This looks wrong and you don't preallocate the mag_spectrum matrix (which is bad practice).

mag_spectrum(:,idx*1024) = abs(spectrum).^2;

mag_spectrum(idx, :) = abs(spectrum).^2; % Use this instead

mag_spectrum(:, idx) = abs(spectrum).^2; % Or this... depending on how you want the data to 'waterfall'
You don't have the proper scaling (see reference paper and code below). Also, I believe the spectogram is displayed in dB.
You're saving the two-sided FFT. This is unnecessary and the spectogram defaults to the one-sided FFT. Use the proper scaling accordingly.

Here's some code used to compute the average PSD of an input sequence using a provided window and some other parameters. You can modify the last few lines to return the entire spectogram by getting rid of the call to 'mean' and properly handling the last step of reindexing.

Performing the overlap is pretty easy if you convert the data vector into a matrix using reshape, and then proceed to concatenate that matrix with a shifted version of itself.

A good explanation of some of the subtleties of the PSD is presented in this paper, especially regarding scaling.

[px f] = welch(randn(1,32768), 1024, ones(1,1024), 512, 1024, 0);
plot(f, 10*log10(abs(px)))

function [Px, f] = welch(x, fs, win, overlap, Nfft, oneside)
% x=vector of data samples
% fs=sampling frequency (in Hz)
% win=vector containing the window
% overlap=number of points to overlap
% Nfft=length of FFT to use
% oneside=a logical variable indicating whether
%    scaling is for 1-sided or 2-sided PSD (true if 1-sided)
% Px=estimate of the power spectral density
% f=vector of frequencies (in Hz) for plotting
if nargin < 6
  error('Insufficient input arguments')
end

x = x(:);
win = win(:);

Nw = length(win);
Ns = length(x);

% Compute the frequency vector and scale factor
if oneside == true
  sfactor = 2 ./ (fs * sum(win.^2));
  f = fs * (0:Nfft/2).'./Nfft;
else
  sfactor = 1 ./ (fs * sum(win.^2));
  f = fs * (-Nfft/2:Nfft/2-1).'./Nfft;
end

% Clip excess data
x = x(1:floor(Ns/(Nw-overlap)) * (Nw-overlap));

% Perform overlap by appending a shifted version
% of the data matrix
y = reshape(x, Nw-overlap, []);
y = [y; circshift(y(1:overlap, :), [0 -1])];
y = bsxfun(@times, y, win);

% get rid of last column
if size(y, 2) > 1
  y = y(:,1:end-1);
end

% Compute the PSD
Px = sfactor * mean( abs(fft(y, Nfft)).^2, 2);

% Return full or onesided result
if oneside == true
  Px = Px(1:end/2+1);
else
  Px = fftshift(Px);
end


advantage of using impulse response function for LTI systems? i have a discrete-time LTI system $L$ that takes input signal $x[n]$ and gives the output signal $y[n]$. since $L$ is LTI, $y[n]$ can be derived as a sum of shifted and scaled impulse responses of $L$. but why is that a useful fact and can you show a practical example where this saves us work?

this is my attempted answer and i would like to know if it's correct: suppose $L$ is an accumulator,

$y[n] = \sum_{k=-\infty}^{n}x[n]$

then by LTI properties, this can be rewritten as a sum of shifted and scaled impulse responses:

$y[n] = \sum_{k=-\infty}{\infty}h(n-k)$

where $h(n-k)$ is the impulse response function.

what does this gain us, intuitively? it seems that this means that if we could only measure the response of $y(n)$ to a single (unit) impulse, we could derive through shifting, scaling and summation $L$'s output to any arbitrary input signal $x[n]$. so we wouldn't have to measure $L$'s output to a more complicated signal, a single pulse would be enough. is this right?

if so how would it work in this example? suppose the input signal $x_0[n]$ was:

$x_0[0] = 1, x_0[k] = 0, \forall k \neq n$

then clearly $L[x_0[n]] = 1$. how can we use this to derive the response to a new input signal $x'[n]$:

$x'[0] = 1, x'[2] = 1, x'[3] = 2$

the result should be 4. how do you formally derive this answer using only the system response to $x_0$?

A sum of shifted and scaled impulse responses would look more like

$ y[n] = \sum_{k=-\infty}^n x[k] \cdot h[n-k] $

The impulse response of an accumulator is

$ h[n] = u[n] $

where $u[n]$ is the unit step function,

$u[n] = \begin{cases}
1 & n \ge 0 \\
0 & \mathrm{otherwise}
\end{cases}$

so for your example case of

$ x[n] = \{1, 1, 2\} $

The result will be

$ \begin{align}
y[0] &= \sum_{k=-\infty}^0 x[k] \cdot h[-k] \\
     &= 1 \cdot 1 + 1 \cdot 0 + 2 \cdot 0 \\
     &= 1 \\
y[2] &= \sum_{k=-\infty}^2 x[k] \cdot h[2-k] \\
     &= 1 \cdot 1 + 1 \cdot 1 + 2 \cdot 1 \\
     &= 4
\end{align}$

and hopefuly you can see how $y[1]$ will be 2, and $y[n], n> 2$ will be 4, as expected.

Your equations are incorrect for the sum of scaled and shifted impulse responses, and for your first input signal $x_0$.


How to reconstruct a time-discrete signal German: Sie betrachten das zeitkontinuierliche Signal $s(t)$ mit einem digitalen Oszilloskop mit Abtastperiode $T_\mathrm{s}$ und messen das zeitdiskrete Signal $\tilde{s}[n] = s(n\:T_s)$

Nehmen Sie an, dass das Abtasttheorem eingehalten wird.

Können Sie $s(t)$ aus $\tilde{s}[n]$ rekonstruieren? Begründen Sie!
Skizzieren sie $s(t)$ für einen Bereich von $- 3\:T_\mathrm{s}\leq t \leq 3 \:T_\mathrm{s}$

English: You consider the time-continuous signal $ s (t) $ with a digital oscilloscope with sampling period $ T_\mathrm{s} $ and measure the time-discrete signal $\tilde{s}[n] = s(n\:T_s)$
Assume that the sampling theorem is met.

Can you reconstruct $ s (t) $ from  $\tilde{s}[n]$? Explain!
Sketch $ s (t) $ for a range from $- 3\:T_\mathrm{s}\leq t \leq 3 \:T_\mathrm{s}$

Sadly, that's all that's given. This is an exam problem from university and I don't have any clue how to get a useful solution.
Maybe some one is able to reconstruct it for me and tell me how it works?

Consider things you know about the sampling theorem:

If it is met, the original signal can always be reconstructed perfectly!
The signal has to be bandlimited for the theorem to be met!

Further, think about the spectral and time domains and their relation to each other:

If a signal is bandlimited in spectral domain, it has to be infinite in time domain.

Then, thinking about the given signal:

The given sampled signal has just zeros before and after the one none-zero value.
It is not periodic.
One can conclude, that in the original signal, there have to be none zero values between all (or at least some) sampling points, otherwise it would not be infinite.

So we have to ask ourselves the question:

Which signal is bandlimited, has zero crossings at regular intervals, a non zero value in the middle and is not periodic?

If you had any lectures on signal theory, the answer should be pretty obvious.


Allpass group delay formula and effect on Karplus-Strong Some introductory questions:

What is the unit of the group delay in a digital filter (is it samples?)
What is $\Omega$ in the equations below (is it $\Omega=2 \pi\frac{f}{f_{sample}} ?) $

If I got any of these wrong, then there is no point reading the rest of this question. Just correct me.

What I want to achieve

I am using a Karplus-Strong algorithm to simulate a steel string. I noticed, that I need to detune the harmonics to higher frequencies, so it really sounds like a steel string (see here). To achive this, I added an allpass to the delay in the feedback loop.

$$H(z) = \frac{a + z^{-1}}{1 + az^{-1}}$$

However the allpass adds an additional delay and detunes my string (lower). So I tried compensating for that additional delay by subtracting the group delay of the allpass from the main delay.

The result was: it still detunes, but in the opposite direction.

The allpass delay:

Matt L.  gave me this formula for the group delay

$$\tau_g(\Omega) = \frac{1 - a^2}{1 + 2a \cos(\Omega) + a^2} $$

$$\Delta t = \frac{\tau_g}{f_{sample}}$$

But robert pointed out that it is the phase delay which determines the pitch and gave me this formula:

$$ \tau_\phi(\Omega) \ = \ \frac{1}{\Omega} \arctan\left(\frac{(1-a^2)\sin(\Omega)}{(1+a^2)\cos(\Omega) + 2a} \right) $$

But alas, when I subtract either of the two from the total delay in the feedback loop (i.e. the regular delay plus the allpass) it appears I am overcompensating. There is no problem with posive $a$ which leaves the fundamental undelayed, but for negative $a$ (e.g. $a=-0.93$) the string is about half a semitone too high.

Edit (1)

I had an error in the formulas which is fixed by now (thanks to Matt L.). Originally my string was tuned too low, after the fix it is too high.

Edit(2)

Meanwhile (thanks to Matt L. and robert) I have little doubt, that my calculations are basically correct. I ran my allpass through scilab and the group delay is reported as apx. 30 samples. With the formulas from Matt and Robert I get values around 31 samples. Hoever, to tune the string correctly (by ear) I must compensate for 24 samples only.

Interestingly Miguel Negrão who did some research on plucked steel strings wrote in his thesis: Generally a negative value of a would have been used, which would have delayed the lower frequencies more than the high frequencies (phase delay), and would require a re-tuning of the main delay in the string loop to keep the string at the same fundamental frequency. In the SupaStrings system case, after some listening tests, a positive value of
a was used instead (non-physical) which keeps the low partials more a less in the same frequency but detunes the high partials down to a lower frequency.

It seems like he faced similar difficulties and couldn't solve them, even though he is quite a brilliant guy.

For the Karplus-Strong algorithm, your primary concern is for phase delay,
$$ \tau_\phi(\omega) \ = \ - \frac{\phi(\omega)}{\omega} $$

rather than group delay.
$$ \tau_g(\omega) \ = \ - \frac{d \phi(\omega)}{d \omega} $$

where
$$ \phi(\omega) \ \triangleq \ \arg \left\{ H(e^{j \omega}) \right\} $$
and $H(z)$ is the transfer function of your filter, and
$$  H(e^{i \omega}) \ = \ \left| H(e^{j \omega}) \right| e^{j \phi(\omega)} $$

For phase delay, it is important to unwrap the phase, which is not necessary for 1st-order APFs.

This is because, for a sinusoid, it's the phase delay that is the actual delay, in units of time, that the filter delays a simple sinusoid.  The group delay is the amount the envelope of the sinusoid is delayed.

For linear-phase filters, $\tau_\phi(\omega)=\tau_g(\omega)$ for all $\omega$ and for any filter $\tau_\phi(0)=\tau_g(0)$.  Phase delay and group delay are not generally the same value in other circumstances.

Now for the 1st-order APF with its pole at $p$:
$$ H(z) \ = \ \frac{1-pz}{z-p}$$
the phase delay in sample units is
$$ \tau_\phi(\omega) \ = \ \frac{1}{\omega} \arctan\left(\frac{(1-p^2)\sin(\omega)}{(1+p^2)\cos(\omega) - 2p} \right) $$
and the group delay in sample units is (same expression MattL gets)
$$ \tau_g(\omega) \ = \ \frac{1-p^2}{1+p^2 - 2p\cos(\omega)} $$

the delay (phase or group) at DC is $\tau_\phi(0) = \tau_g(0) = \frac{1+p}{1-p}$ and we normally have $-1 \le p \le 0$ to get a delay of between 0 and 1 sample.

Now, in the Karplus-Strong algorithm, the frequency of every harmonic (perhaps better called a "partial"), including the 1st harmonic (which we'll call the "fundamental") is the reciprocal of the total loop delay (phase delay) at the frequency of that partial.  That is a circular mathematical relationship.  For the $k$th partial (and assuming normalized sampling frequency $f_s=1$):
$$ f_k \ = \ \frac{k}{n + \tau_\phi(2 \pi f_k)} $$
where $n$ are the integer number of samples in the delay line that is in series with the APF.  Given your fundamental frequency $f_1$, the value of $n$ is $$ n \ = \ \left\lfloor \frac{1}{f_1} \right\rfloor \ = \ \mbox{floor}\left(\frac{1}{f_1}\right) $$
 (Be sure to add 1 sample that is necessary in a digital filter with feedback.)

Now, often we approximate using the DC delay, but to design this exactly, you need to somehow solve that frequency-dependent phase delay at the frequency of each partial.  That's a bitch.  So you would need to solve
$$ \frac{k}{f_k} \ = \ n + \tau_\phi(2 \pi f_k) \ = \ n + \frac{1}{2 \pi f_k} \arctan\left(\frac{(1-p^2)\sin(2 \pi f_k)}{(1+p^2)\cos(2 \pi f_k) - 2p} \right) $$
for $f_k$ given every integer $k$.


What would be the best way to estimate velocity from a video? I am not from an EE background, so I am looking for answers to some open-ended questions.

I want to calculate how fast water is flowing during a flood. With or maybe even without debris in it. I have footage video from which I can extract sec by sec images of changing landscape.

It can be considered analogous to a car velocity estimation problem, except the "object" in the case of flood becomes much more complicated.

So far I know about Optical Flow: http://www.mathworks.com/discovery/optical-flow.html

What are the latest advances in the field that might help me to accomplish my goal?

Any tips at all?

You can probably address that by tracking debris across frames, assuming that they move at the speed of water and aren't too far apart in successive frames.

Maybe start with the old stuff: Shi, J.; Tomasi, C., "Good features to track," in Computer Vision and Pattern Recognition, 1994. Proceedings CVPR '94., 1994 IEEE Computer Society Conference on , vol., no., pp.593-600, 21-23 Jun 1994.


Bandwidth required to transmit a binary sequence Does a digital binary number have any bandwidth? I explain suppose T want to express the number 10  in binary digits:  1 0 1 0 would that digital number have any BW in frequency domain? I  would say the question has no meaning as frequency is the dual in time and a number 5 whether expressed in decimal or analog is not a quantity changing with time. am I right?

It's not meaningful, in my opinion, to talk of the "bandwidth" of a constant number.  Unless you want to characterize that constant number as DC (because the value does not change in time).  And the bandwidth of DC is zero.

But a stream of bits toggling on and off, yes, a bandwidth (which is the rate of change of the bits or "bitrate", if we're considering two-sided bandwidth or half the bitrate if we're considering the single-sided bandwidth) is meaningful.

In fact, a very important fundamental theorem from Claude Shannon relates bandwidth, signal-to-noise ratio, and rate of information together.


Filter design for frequencies in range 0-60Mhz First of all excuse me for asking a seemingly vague question. But i am going nowhere in the problem i am trying to solve.

I am working on a powerline communication system where i need to upconvert in 28-56Mhz range and later downconvert it to 0-28Mhz so that it works in the coupler range.So i need to design a band pass filter in 28-56Mhz and a low pass filter in the 0-28Mhz range. I need to design an analog filter for the same. I have no experience of filter design. I have been looking online for resources but wherever i go there is an explanation on chebyshev, butterworth and  elliptic filters.Then you have 1st order,2nd order and higher order filters equations. My problem is from all this i am  unable to decide which type of filter to use and what order to use. So i need to understand practical design constraints and examples of how to design given certain constraints and then boiling down to a solution. So could anybody help me in how to approach the filter design problem and suggest either examples or resources to go forward.

As you stated there are some parameters associated with filters such as type of filter (butterworth, elliptic, chebyshev I,chebyshev II, bessel,...), order of filter, etc. (Note that there are lots of considerations in filter design bases on the application)

Let's first talk about the order of filters:
Filters with higher orders let you have a filter with better quality. For example in lowpass filter if you want to have a sharper transition band you need to design a filter with high order filers. Order of a passive analog filter is determined by the number of elements you use in filter (resistor, inductor, capacitors). This was all I knew about the order of the filter.

Now let's talk about the most applicable kinds of filter types. The most important difference between filter types is the band they have ripples. Note that the priority of filters is due to the application and when we say one filter is better than the other we must state that in what aspects (order, ripples, poles, zeros, stability, group delay, phase response, minimum phase, ...) this is better than that one. Picture below compares 4 important types of filters designed with the same order. As you can see, Elliptic filter has ripples both in pass-band and stop-band while it has the sharpest transition, and butterwoth filter has no ripples but has least sharpness among filters of the same order.


understanding bases and frames for Gabor transform For the 2D discrete Gabor transform, why is it that we cannot use a set of orthonormal basis for its representation, instead we have to use frames for representing it?

The frames end up giving you an orthogonal representation, just not exactly a basis representation that you're talking about. Simply by definition that's not what a Gabor filter does. It's a time-frequency representation of your signal, while "orthogonal basis vectors" wouldn't have a time-specific property -- they stretch for as long as your signal stretches.

If you want to use orthogonal basis transforms, use Fourier, cosine, empirical mode decomposition and similar transforms.

If you want time-frequency representation, your can use short-time Fourier transform, wavelet transforms (of which Gabor can be considered a special case) and other similar transforms.


What is the difference between the spectra of white Gaussian noise and of a Dirac impulse? Perhaps it is a basic question, but

can anyone explain the difference between the spectra (or power spectral density) of white Gaussian noise and a Dirac delta?

They are both constant over the whole frequency range, however there are is a clearly striking difference between these signals.

Why do they have the same spectra?

They have very different phase spectra. For a dirac the phase of the spectrum is zero for white noise, it's uniformly distributed over  $[0, 2 \pi]$

They are both constant over the whole frequency range

That's only true for the spectrum of the Dirac. For white noise, it's more complicated. The FFT of a chunk of white noise does NOT have constant spectrum, it's actually white noise as well (the real part and the imaginary part). It only becomes constant if you average every over a large enough number of chunks (or a long enough period of time) or if you apply some sort of spectral smoothing.


Key detection Dance Music inaccuracies - matlab I have been working on a project with a vide variety of functions for DJ'ing. I am nearing the end, however, the accuracy rate on the key detection part of it is not pleasant when fed tracks of Dance music. With classical music, the algorithm returns a very better accuracy rate.

First of all I must explain how I never studied Signal Processing at University and I had to learn the concepts on my own, so excuse my lack of extended knowledge in the subject (hence why I am asking here).

I spent most of my time researching and at the end I came up with an algorithm that is a mixture of techniques suggested in literature.

Summary of Implementation: The audio is split into 5.5 second sections and translated to mono, downsampled to 11025hz and translated into the frequency domain using STFT technique. The STFT result in hz is translated into real music notes by a Chromagram. The chromagram is a 12 dimensional chroma vector in which each bin represents a real pitch class i.e. the first bin represents A, the second represents A#, third represents B etc. Each chroma vector is then correlated against 24 key templates. A coefficient is obtained for every key and a weighting system rewards the highest correlating keys when the difference between the other correlating keys is significant, and penalises when the difference is minute.

I perform the STFT and Chromagram steps with the help of a package by Dan Ellis, that I found online.

Now in the original code the parameters are: w (4086) is the basic STFT DFT length (window is half, hop is 1/4)

In the paper written by the same person there is this section: http://imgur.com/0EVyjFc, where he explains how the STFT is performed.

First question: I am not understanding why he divides the window by half before passing it to the STFT function?

Now, I am using an FFT length of 8192 samples unlike the 4086 of Ellis, after performing parameter testing which gave me best results at that length.

Before the hanning window is calculated these are the values of the parameters..

sr =

11025

N =

8192

W =

4096

 H =

2048

The second question: what could be the reason for such inconsistent and inaccurate results with dance music? Is it perhaps some characteristic of dance music that makes it difficult for the chromagram to match the notes in hz from the STFT to real music notes? Perhaps because of the Kicks, or something to do with the lower frequency ranges? What do you suggest can be done to the audio prior to performing STFT? may be some kind of filter? Can you suggest trying any different parameter in the STFT? or in general any other suggestion / extra step to make the final detection better? Suggestions that would alter the main steps of the key detection completely would be a bit out of reach, however, as I dont have much time left to submit the project.

Thanks in advance!

A 12-dimensional chromagram may not estimate pitch very well, given synthesized semi-tonal sounds that include lots of energy at multiple odd harmonics, as those harmonics will often end up in the wrong pitch class slots.  If there's lots of tonal bass, low-pass filtering just above the primary range of the bass instrument/synth's fundamental might help getting rid of some of the harmonics and overtones.

If the drums are not tuned to the key (probably not), then they will interfere with any key centered pitch class information.  Perhaps time-domain filtering out of all percussive events, and only analyzing the tones between the drum beat attacks might help in this case.  Even a simple low-envelope amplitude gate function might be suitable (assuming the sound isn't so compressed that there is no lower amplitude in any extent portion of the envelope!).

And, as per the previous answer, are you sure the sounds of your "dance music" actually contains a musical key? (as opposed to producing some cognitive or auditory illusion of tonality).


Frequency Modulation to arbitrary Audio signals I understand how FM synthesis works, but I was wondering how to FM arbitrary samples or part of samples. I.E any frame of audio samples.

Further more the modulator can itself be an arbitrary frame of samples itself.

How is this done in an algorithmic manner?

I have, in the past, used delay modulation on arbitrary audio signals.  but the modulation waveform itself was not arbitrary, but was a sinusoid.  but we could speed this modulation oscillator up to an audio rate.

Now that I think of it, I tried hooking that modulation oscillator up to a pitch detector.  This is a quarter century ago, so I don't remember everything.

I s'pose you could derive that delay modulation signal from the input and likely LPF it before applying it to the instantaneous delay of the signal going out.


Correlation Using FFT / IFFT (Convolution in Frequency Domain) in Java I try to find about the delay between two audio files using Cross Correlation in Java. I've already done this algorithm so far that i get a idea about how many samples is the delay.

FFT x1 -> Zero Padding to length: x1.length() + x2.length()
FFT x2 -> Zero Padding to length: x1.length() + x2.length()
X3 = FFT{x1} multiply FFT{2}*
x3 = IFFT{X3}

To validate my results, i use octave. So i can see that my result is in generally correct, but the position of my graph isn't correct.

For testing i use the same files, so that this represent auto correlation function. In octave graph the maximum correlation is in the center of output vector, which seems to be right. -> so the delay is 0, because same files

In my version, the graph is at the beginning and the end.

So if i use two different signals, in octave the maximum correlation moves left/right from center. In my algorithm its going from beginning to middle or end to middle.

So it seems my graph is calculated correct, because the peak moves the amount of samples that is the difference between the files. But the position isnt correct.

It looks like i have to modify my output vector: left half <-> right half

Any ideas whats my problem ?

Edit:

public void execFastCrossCorrelation() throws FFTException
{
    // the used block size is N(x1) + N(x2)
    int blockSize = (x1.length+x2.length-1);

    // create & calc FFT with zero padding & get result

    // FFT x1
    FFT f1 = new FFT(x1, blockSize, 44100);
    f1.execFFT();
    Complex[] c1 = f1.getOutputData();

    // FFT x2
    FFT f2 = new FFT(x2, blockSize, 44100);
    f2.execFFT();
    Complex[] c2 = f2.getOutputData();

    // create output for cross correlation by using fast convolution : c3 = c1 multiply c2*
    Complex[] c3 = new Complex[c1.length];
    for(int i=0; i<c1.length;i++)
    {
        c3[i] = c1[i].multiply(c2[i].conjugate());
    }

    // create & calc IFFT  & get result
    IFFT f3 = new IFFT(c3, c3.length, 44100);
    f3.execIFFT();
    Complex[] corrResult = f3.getOutputData();

// search the maximum corrleation coefficient and his index
    double maxVal = 0;
    int maxIndex = 0;

    for(int i=0; i<corrResult.length;i++)
    {
        if(corrResult[i].abs() > maxVal)
        {
            maxVal = corrResult[i].abs();
            maxIndex = i;
        }
    }

So FFT class is done by using apache-commons-math3. Simply use the input vecotr and create a FFT for the given blockSize (rounded up/down to power of two and zero padded if necessary).

Edit:

So now i tried to flip the output (means left halft -> right half and right half -> left half)

The result matches the version of octave.

If i try this with different files, its also working like in octave, means the maximum is moving from the center to left/right. That means that my idea that the output vector has to be "mirrored" seems to be right. Do you have any idea why?

A real (non-complex) time-domain input will always give a FFT that is duplicated (mirrored)

Some implementations will perform a shift to put the DC (frequency=0Hz) term at index 0. Otherwise you will find the DC term halfway through the array.


Discrete-time Fourier transform of $a^{|n|} u[n]$ I have a problem calculating the DTFT of this pair:
Could anyone tell me why the DTFT for $a^{|n|} u[n]$ is different from $a^{n} u[n]$'s?

There's an obvious typo in that DTFT pair, because - as you've probably noted yourself - we clearly have $\alpha^{|n|}u[n]=\alpha^nu[n]$.
The given expression for the DTFT corresponds to the two-sided time-domain sequence $\alpha^{|n|}$, if $|\alpha|<1$ holds.


Why does my stop-band filter attenuate other frequencies? Here is my code for a two-toned signal, where I use a stop-band to remove the higher tone and then plotted the before (in blue) and after (in red) in the frequency domain after convolving my signal with the filter coefficients.

If you put this code into Matlab you can clearly see that the higher frequency has been successfully removed by the filter, yet for some reason the amplitude of the lower frequency has been cut in half, and the more I increase the # of filter coefficients, the more it just flattens my whole curve, why does this occur?  And how can I prevent it so that the stop-band does not propagate outwards? Here is the image and the code:

fSampling = 8000;
tSampling = 1/fSampling;

t = 0:tSampling:0.005;
F0 = 1000;
F1 = 3000;

xt = sin(2*pi*F0*t) + sin(2*pi*F1*t);
ht = fir1(40,.25,'stop');
yt = conv(xt,ht);

fAxis = -4000:125:4000-125;

xF = fft(xt,64);
MagXF = fftshift(abs(xF));

plot(fAxis,MagXF);
hold on

yF = fft(yt,64);
MagYF = fftshift(abs(yF));

plot(fAxis,MagYF,'r')

Filters always have an inherent "roll-off" in their frequency response, because you can't practically realize a pass-band that's a perfect rectangular function. For a low-pass filter, the point where the magnitude of your frequency response drops to -3dB is called the pass-band and everything beyond is called the stop-band (technically, everything beyond the corner frequency, but we'll take the corner frequency to be the -3dB level). How fast your frequency response attenuates beyond the pass-band is dependent on the filter length.

If you look at the frequency response of your filter, ht, you'll see that it drops to -6dB at 1000 Hz:

So it makes sense that the power drops by 6dB after filtering, which you see in your figure as a halving in amplitude.

If you had looked at the documentation for the fir1 function that you used, you would've realized this too (emphasis mine):

  B = fir1(N,Wn) designs an N'th order lowpass FIR digital filter
      and returns the filter coefficients in length N+1 vector B.
      The cut-off frequency Wn must be between 0 < Wn < 1.0, with 1.0
      corresponding to half the sample rate.  The filter B is real and
      has linear phase.  The normalized gain of the filter at Wn is
      -6 dB.

Now to create sharper filters, with responses that get pretty close to a rectangle, you'll have to use IIR filters, which come with their own host of problems with stability, etc., but definitely an option. You can see my answer here for some ideas on implementing a discrete-form 2, second-order section IIR filter that gives very sharp corners. The example there is for a band-pass filter, but you can read up the docs for the functions used there and implement a low-pass version yourself.


Why is an RC circuit a first order system? I am going through the following tutorial: Time-Response Identification of a Resistor–Capacitor (RC) Circuit

Why do we call the RC circuit a first order system?

I understand the order of a system is based on the highest derivative involved in its governing equation. The input governing equation of an RC circuit is described by equation 1 in the tutorial:

$$e_i-iR-\frac{1}{C}\int i \;dt =0$$

I understand the equation. But, I'm confused because there is no derivative term in it and yet it is called a first order system.

The order of system gives the number of memory storing elements in the circuit.

What you have written above is a KVL equation of the circuit. But when you determine the order of the system you need to find the transfer function of the system. Laplace domain is normally used. The highest power of $s$ in  the denominator in a strictly proper transfer function gives the order of the system.

So in the $s$-domain,say you take the output across the capacitor,

$\dfrac{V_o(s)}{V_i(s)} = \dfrac{\frac{1}{sC}}{R+\frac{1}{sC}}$

simplifying we get

$\dfrac{V_o(s)}{V_i(s)}= H(s)=\dfrac{\frac{1}{RC}}{s+\frac{1}{RC}}$

So the highest power of $s$ is given by 1. So the order of the system is 1.


a neural network approach for FIR filter I am trying to write a code for a neural network to do the digital filtering on some signals. Is there any neural network model for digital filtering?

Yes. The FIR filter model you're used to is a series of Neurons with weighted inputs, and a linear activation function.
In other words, a standard FIR filter is a neural network.
I mean, it's called "CNN". The C is exactly the operation a filter does.


Processing Accelerometer Vibration Data With FFT For a project I am attempting to use an accelerometer to measure vibration in an RC aircraft and determine frequency from the result. To do this I am using an MPU-6000 accelerometer sampling @ 1000Hz. Before I actually deploy this system on an aircraft, I have taped the sensor to a speaker and used a frequency generator to generate a 100Hz tone for 30 seconds. This should be well within the accelerometer's range:

1000Hz / 2 = 500Hz maximum by Nyquist rule, much higher than 100Hz

My amplitude data in the time domain is bellow (generated w/ GNU Octave):

Zoomed in:

That looks like a ton of data to me but nothing seems to be obviously in error. Bellow is the data after performing a simple periodogram in GNU Octave with the command periodogram(DATA,[],[],1000):

A peak can be seen at 100Hz but there are all of those other spikes that I can't explain. I have tried to detrend my data using the GNU Octave command detrend(DATA,'linear') but that had no effect on removing the stray peaks.

I get the feeling that I should be doing some filtering but I am new to this and really don't know what I am doing. How can I remove all of the stray peaks in the periodogram? Thanks in advance!

Data file on Google Drive: http://goo.gl/w3Brol

EDIT:

Here is the periodogram for the unit sitting on the speaker with nothing playing for 30 seconds:

There doesn't appear to be any interference. The speaker was on, just not playing anything. This is making me think that my speaker is generating harmonics when I supply a sine wave. Could that be? As always, thanks for the help!

Silence data file on Google Drive: http://goo.gl/bNbnDA

EDIT #2:

Success! After isolating my speaker from the table, taping the device to the back of the speaker (I bet glue would help further), and lowering the volume I was able to get the result below! You guys were spot on, thanks for all your help!

After taking a closer look at your results, here is what I think.

Because you are playing back sinusoid through your speaker you are also observing it's properties. It has some significant harmonic distortions appearing as the harmonics of fundamental frequency: $200$, $300$, $400 \mathrm{Hz}$ in your case. That is either because:

You are driving your speaker at high volume
It is simply crap for accurate measurements (sorry...)

Obviously you could also get harmonic distortions coming from your sound card, but this is unlikely.

On the other hand what you see at $50 \mathrm{Hz}$ is the mains hum of AC voltage (you are definitely not in USA). This distortion also has it's harmonics at $100$, $150$, $200 \mathrm{Hz}$, etc. So you can see that they do overlap with your test signal.

Nevertheless if you look at higher frequencies, i.e. around $400 \mathrm{Hz}$, you will see two peaks next to each other. Those are two harmonics, one of your test signal and other one of your mains hum. Because there was a slight offset at low frequencies it diverges even more at higher frequencies.

In order to confirm the mains hum you should simply record silence. If it is present, then well - you have a slight problem. If you are using laptop, try to unplug the AC adaptor - it should help. Otherwise you are simply picking up everything with your power supply. You can also try external audio interface. On the contrary, if you are already using some kind of interface with XLR microphone or whatever - try to break the ground loop.

When it comes to speaker distortions try to use other one or maybe use it a slightly lower volume.


How to show that this filter is a HP filter i got a system with transfer function given by:

$$H(\omega)=1-e^{-j\omega}$$
I already plot it, and that seems to be a periodic function with $H(0)=0$, $H(\pi)=2$, , is that enough to show that this is a FIR bandpass filter. thank you very much

As was pointed out in the comments above, this is a highpass filter. There is no strict definition for what constitutes a lowpass, bandpass, or highpass filter. Instead, they possess certain qualities:

Lowpass filters attenuate high frequencies while passing low frequencies.
Bandpass filters attenuate low and high frequencies, while passing some band in the middle.
Highpass filters attenuate low frequencies while passing high frequencies.

Note the lack of precise terms in the descriptions above. In order to classify a filter into a category like those above, one evaluates the character of its magnitude transfer function. That requires a little manipulation:

$$
H(\omega) = 1 - e^{-j\omega}
$$

$$
\begin{align}
|H(\omega)| &= |1 - e^{-j\omega}| \\
&= \left|e^{-j\frac{\omega}{2}} \left(e^{j\frac{\omega}{2}} - e^{j\frac{\omega}{2}}\right)\right|
\end{align}
$$

Using Euler's identity, we can simplify the expression inside the parentheses:

$$
e^{j\frac{\omega}{2}} - e^{j\frac{\omega}{2}} = j2\sin\left(\frac{\omega}{2}\right)
$$

$$
|H(\omega)| = \left|e^{-j\frac{\omega}{2}} j2\sin\left(\frac{\omega}{2}\right)\right|
$$

We can drop the $j$ and $e^{-j\frac{\omega}{2}}$ terms here, because each of them have unity magnitude and therefore don't affect $|H(\omega)|$. Therefore:

$$
|H(\omega)| = \left|2\sin\left(\frac{\omega}{2}\right)\right|
$$

We can check this result by comparing it with a plot from MATLAB's Filter Visualization Tool (fvtool):

Which looks about right. This is plotted over the region $\omega \in [0, \pi]$, and shows that the filter has zero response at DC (attenuating low frequencies) and a magnitude response of $2$ at $\omega = \pi$ (passing high frequencies). That looks like a highpass filter.


Convolution reverb using a frequently updated frequency dependant impulse response I am implementing a software to auralize a virtual scene but am kind of a newbie to DSP. The focus should be the calculation of the IR which is frequency dependant. But this means that the renderer of the sound should be able to convolute reverb in a frequency dependant way (in realtime possibly). A IR has a set of Reflections. A reflection is a tuple of delay (double) and decay for each frequency band to a decay in dB (hash_map from freqband to double).

Now I want to convolute this IR with my input signal in the most efficient manner. What I intend to do is using the FFT to get the frequency spectrum and then apply a reverb function to this signal. The two main conceptual problems for me here are:

How do I apply delay in the frequency spectrum? (Is it even necessary to do this in the frequency spectrum?)
Is there a common method to apply such a reverb function to a signal? And if there is, do I need to have phase information modelled?

What I was thinking about was splitting the IR into several frequency responses dependant on the frequency band and then apply it to the frequency spectrum of the input signal. But then again this involves problem 1 and somehow includes phases because I would not know otherwise if I need to add or subtract.

There are basically two different methods to create reverb: you can store (simulated or measured) room impulses and apply it by convolution. Reverb impulse response are very, typically thousands of tabs. This makes direct convolution impractical and the method choice is overlap add or related FFT-framed algorithms. This will induce latency and and is difficult to update.

The alternative would be a parametric algorithm. Typical building blocks include tap-delay lines, schroeder allpass filters, comb filters and feedback delay networks. Most of these blocks can be tied to some meaningful room acoustical parameter: reverb time as a function of frequency, refection density, mean free path length, etc. They also can be updated quickly and continuously. However they have a tendency to sound metallic and artificial.

Building good sounding reverb is hard. Many a thesis has been written on the topic and commercial reverb algorithms and/or devices often take hundreds of hours of development + tuning. I recommend doing some research and then deciding how deep you want to go in.


Generate Complex White Gaussian Noise in MATLAB In the followinf code I am trying to generate a Complex Gaussian Noise:
n_3 = sqrt(0.1)*randn(1,K);
n_4 = sqrt(0.1)*randn(1,K);
beta_NLoS = (n_3+1i*n_4); % CN(0,0.1)

Does my code do as intended?

If you want a Circular Complex Gaussian Noise (Independent):

vComplexNoise = sqrt(noiseVar / 2) * (randn(1, numSamples) + (1i * randn(1, numSamples)))

For correlated noise you'll need to define the Co Variance Matrix and use Cholesky Decomposition.

Update

Following @Stanley Pawlukiewicz advise, run the following code:

numSamples = 100000;
noiseVar   = 4;

mA = sqrt(noiseVar / 2) * (randn(numSamples, 1) + (1i * randn(numSamples, 1)));

var(mA)

You should see result which is very close to noiseVar on the screen.


What subclass of nonlinear systems can be represented by Volterra series? Can every possible nonlinear time-invariant Volterra series be represented by a Volterra series?
If not, which ones can?
Also, are there any time-variant systems that can be represented by Volterra series? I don't see how from the equation, but want to be sure I'm not missing something.

The system must be time invariant and smooth in the functional derivative sense. That doesn't guarantee that the Volterra series converges (like with Taylor series, there are pathological counter examples), but almost all systems that have these properties have a convergent Volterra series.

The problem in practice is that the required expansion order for good approximations of your system is often far too great to provide for a useful implementation of the system. A condition for quick convergence is, that the functional derivatives are well bounded and decay quickly with order.

A little more detail:

Consider a non-linear system $f:V\to W$, where $W$ is an inner product space with induced norm and $V$ is a normed vector space. With this, look at a map

$$\begin{align}
W\times V&\to K_W \\
(w,x) &\mapsto \langle w,f(v+\lambda x)\rangle
\end{align}$$

This map is anti-linear in $w$ and non-linear in $x$. We treat $\lambda \in K_V$ and $v \in V$ as parameters. Further below, $v$ will be the center of the Volterra expansion and $\lambda$ will reveal the different expansion orders.

If we hold $w$ constant, we get a functional

$$
\begin{align}
F_w:V&\to K_W\\
x&\mapsto F_w[x]=\langle w,f(v+\lambda x)\rangle
\end{align}
$$

For bounded $x$, this functional can be expanded into a series of powers of $\lambda$ around $v$:
$$F[x]=\sum_{n=0}^\infty \frac{\lambda^n}{n!}F^{(n)}(x)$$

where $F^{(n)}(x)$ is the $n$-th (functional) derivative of $F$ around $v$ in direction $x$.

The assumption that all $F^{(n)}$ exist and the series converges for all finite norm choices of $w$ is what we mean, when we say: "the system $f$ is smooth".

Going from here to the full Volterra series is not too difficult. The key is understanding, how the functional derivatives relate to $w$ and $x$. Since $F$ is anti-linear in $w$, every single term in the series expansion must be, and therefore $F^{(n)}$ as the only factor with dependence on the choice of $w$. It is also clear by the option to absorb $\lambda$ into $x$, that the derivative $F^{(n)}$ must be of $n$-th order in $x$.

Hence we can write $F^{(n)}(x)=\langle w, G_n(x) \rangle$ with multilinear $\hat{G}_n:V^n \to W$ as well as an embedding $E_n:V \to V^n$ with $E_n(v)=(v,v,\dots,v)$ composed to $G_n=\hat{G}_n \circ E_n$.

In infinite dimensional vector spaces with uncountable bases, these multilinear maps are expressible through nested integrals over kernel functions. If you use the time basis for this expansion and make use time-invariance, you get nested convolution integrals and arrive at the Volterra series. I'll leave this mostly notational part to you. If you need help with it, I can also elaborate more.

As the last step, you use the properties of the inner product to draw $w$ out of the series expression and simply drop it, as the resulting equation needs to hold for all proper $w$.


Improving SNR using DSP techniques I am building an optical OOK (On off Key Ring) system without a carrier frequency. [However I have a guard time in between symbols, so a consecutive "1" messages will result in a pulse train as opposed to DC, see image].  Essentially, the presence of signal indicates one and lack of it indicates zero. I have a precise clock that is syncs the receiver with transmitter. The system operates with a low SNR and I like to improve SNR using DSP techniques.

I have a few questions:

I do selective sampling in my hardware, in other words, I do not continuously sample the channel but only sample when the probability of seeing the signal is maximum (i.e. this is a light pulse, I time the ADC such that ADC samples at the end of the pulse where I know the whole analog chain is stabilized). See image.

Naturally, this drawing doesn't show the noise but it is there. This is particularly a  low signal system and primary noise sources are shot noise, johnson noise and amplifiers internal noise. (optical system so there are no other interferer except Sun). My observation of noise indicates, it is similar in all frequencies. (At least what I see in Scope)

I now use a simple threshold comparison in software to determine if the data is one or zero.  Is there a better way? I have thought of some options but I like to hear from the experts.

So far I have considered the following options:

Do continuous ADC and try to integrate during rise time:Not entirely
sure on the benefit (There may be other benefits, I don't know).
Matched filter in software: Don't really understand the math but
based on what I read, a possibility
Sample during the guard time and subtract this from the signal ADC
value (This may    provide some further details but also not so sure,
guard time would be the noise measurement)
Change hardware to a synchronous decoder, costly, time consuming and
may not work well    since my data rate is fast and getting a
synchronous demodulator would mean expensive board since I have to
built a multi-MHz carrier frequency system.

Since you indicated that the power spectrum of your background noise is flat, I'll assume it is white. A major drawback with your current approach is that you are discarding a large amount of the signal power; even with the effect of front-end bandlimiting shown by in your diagram by the exponential-rise step response, a single ADC sample near the end of the rounded pulse provides a snapshot of the receiver input that is rather localized in time. You can take advantage of more of the signal power by sampling at a higher rate and applying a matched filter at the higher sample rate.

Theory:

You can look at this as a relatively simple problem in detection theory. In each symbol interval, your receiver needs to decide between two hypotheses:

$$
\begin{eqnarray*}
H_0 &:& signal \ is \ not \ present \\
H_1 &:& signal \ is \ present
\end{eqnarray*}
$$

This sort of problem is often solved using Bayesian decision rules, which attempt to make the optimum decision according to some specified measure of risk. This provides a framework where one can optimally make detection decisions based on a flexible set of criteria. For example, if there is a large penalty to your system for failing to detect the signal if it is in fact present (i.e. you pick $H_0$ when $H_1$ is true), then you can build that into your decision rule if needed.

For a detection problem such as yours, where you're trying to decide between zeros and ones at the receiver output, the penalty is typically assumed to be equal (outputting a zero when a one was transmitted, and vice versa, "hurt equally"). The Bayesian approach in that case reduces to a maximum-likelihood estimator (also described here): you pick the hypothesis that is most likely, given the observation that your receiver makes. That is, if the quantity that your receiver observes is $x$, then it would generate a decision based on the hypothesis that has the largest likelihood function value. For the binary decision case, the likelihood ratio can be used instead:

$$
\Lambda(x) = \frac{P(x\ |\ H_0 \ is \ true)}{P(x\ |\ H_1 \ is \ true)} = \frac{P(x\ |\ signal \ is \ not \ present)}{P(x\ |\ signal \ is \ present)}
$$

Using the above model, for each observation of the channel $x$, the optimal receiver would decide that the signal was not present (therefore outputting a zero) if the likelihood ratio $\Lambda(x)$ is greater than one (and therefore the signal was most likely to be not present based on the observation), and vice versa.

What remains is a model for the signal of interest and any other components in the receiver detection statistic $x$ that could affect its decisions. For a digital communications like this, it might be modeled as follows:

$$
\begin{eqnarray*}
H_0 &:& x = N \\
H_1 &:& x = s + N
\end{eqnarray*}
$$

where $n$ is a random variable taken from some distribution (often assumed to be zero-mean Gaussian) and $s$ is a deterministic component of the observation that is due to the signal that you're looking for. The distribution of the receiver observable $x$, therefore, varies depending upon whether hypothesis $H_0$ or $H_1$ is true. In order to evaluate the likelihood ratio, you need a model for what those distributions are. For the Gaussian case referenced above, the mathematics look like this:

$$
\Lambda(x) = \frac{P(x\ |\ H_0 \ is \ true)}{P(x\ |\ H_1 \ is \ true)} = \frac{P(x\ |\ x = N)}{P(x\ |\ x = s + N)}
$$

$$
\Lambda(x) = \frac{P(x\ |\ H_0 \ is \ true)}{P(x\ |\ H_1 \ is \ true)} = \frac{e^{\frac{-x^2}{2 \sigma^2}}}{e^{\frac{-(x - s)^2}{2 \sigma^2}}}
$$

where $\sigma^2$ is the variance of the Gaussian noise term. Note that the additive signal component only has the function of shifting the mean of the resulting Gaussian distribution of $x$. The log-likelihood ratio can be used to get rid of the exponentials:

$$
\ln(\Lambda(x)) = \ln\left(\frac{e^{\frac{-x^2}{2 \sigma^2}}}{e^{\frac{-(x - s)^2}{2 \sigma^2}}}\right) = \left(\frac{-x^2}{2 \sigma^2}\right) - \left(\frac{-(x - s)^2}{2 \sigma^2}\right)
$$

Recall that our decision rule picked $H_0$ if the likelihood ratio was greater than one. The equivalent log-likelihood decision rule is to pick $H_0$ if the log-likelihood is greater than zero. Some algebra shows that the decision rule reduces to:

$$
\begin{eqnarray*}
x < \frac{s}{2} \rightarrow choose\ H_0 \\
x > \frac{s}{2} \rightarrow choose\ H_1
\end{eqnarray*}
$$

Note that if $x = \frac{s}{2}$, then both hypotheses are equally likely, and you would need to just pick one; this isn't a practical concern for continuously-valued signals, however. So, given a known signal amplitude $s$, we can detect its presence against a background of Gaussian noise optimally by setting a threshold $T = \frac{s}{2}$; if the observed value $x$ is greater than $T$, we declare the signal present and emit a one, and vice versa.

Practice:

There are a few practical issues that creep into this simple, toy theoretical example. One: just mapping the scenario that you described into a deceptively simple-looking model might not seem straightforward. Secondly, it's very rare that you would know the amplitude $s$ of the signal that you're looking for, so threshold selection requires some thought.

As I referenced before, noise is often assumed to be Gaussian because the normal distribution is so easy to work with: the sum of a bunch of independent Gaussians is still Gaussian, and their mean and variances just add also. Also, the first- and second-order statistics of the distribution are enough to completely characterize them (given the mean and variance of a Gaussian distribution, you can write its pdf). So, hopefully that's a decent approximation at least for your application.

There are two ways to improve the performance of the detector given the model described above: you can increase $s$ (i.e. increase the signal power), making it stand out more against the noise. You could decrease $N$ (i.e. reduce the amount of noise), reducing the amount of interference that makes the presence of $s$ unclear. Or, equivalently, you can think of the signal to noise ratio instead. To see its importance, let's go back to the theory for a second. What is the probability of a bit error given our decision rule?

$$
\begin{eqnarray*}
P_e &=& P(choose \ H_0 \ |\ H_1\ true) P(H_1\ true) + P(choose \ H_1 \ |\ H_0\ true) P(H_0\ true) \\
&=& \frac{1}{2} P(x < \frac{s}{2} \ |\ x = s + N) + \frac{1}{2} P(x > \frac{s}{2} \ |\ x = N) \\
&=& \frac{1}{2} F_{x\ |\ x = s + N}\left(\frac{s}{2}\right) + \frac{1}{2} \left(1 - F_{x\ |\ x = N}\left(\frac{s}{2}\right)\right)
\end{eqnarray*}
$$

where $F_{x\ |\ x = s + N}(z)$ is the cumulative distribution function of the distribution of the observation $x$, given that $x = s + N$ (and likewise for the other function). Substituting in the cdf for the Gaussian distribution, we get:

$$
\begin{eqnarray*}
P_e &=& \frac{1}{2} \left(1 - Q\left(\frac{\frac{s}{2} - s}{\sigma}\right)\right) + \frac{1}{2}  Q\left(\frac{\frac{s}{2}}{\sigma}\right) \\
&=& \frac{1}{2}  + \frac{1}{2} \left(-Q\left(\frac{\frac{s}{2} - s}{\sigma}\right) + Q\left(\frac{\frac{s}{2}}{\sigma}\right)\right) \\
&=& \frac{1}{2} + \frac{1}{2} \left(-Q\left(\frac{-s}{2\sigma}\right) + Q\left(\frac{s}{2\sigma}\right)\right) \\
&=& \frac{1}{2} + \frac{1}{2} \left(-Q\left(\frac{-SNR}{2}\right) + Q\left(\frac{SNR}{2}\right)\right) \\
&=& Q\left(\frac{SNR}{2}\right)
\end{eqnarray*}
$$

where $Q(x)$ is the Q function:

$$
Q(x) = \frac{1}{\sqrt{2 \pi}} \int_{x}^{\infty} e^{\frac{-z^2}{2}} dz
$$

(i.e. the tail integral of the standard normal distribution's pdf, or $1$ minus the distribution's cdf) and $SNR$ is the signal-to-noise ratio $\frac{s}{\sigma}$. The above function is a strictly decreasing function of $SNR$; as you increase the ratio of the signal amplitude $s$ to the noise standard deviation $\sigma$, the probability of making a bit decision error decreases. So, it behooves you to do whatever you can to increase this ratio.

Remember our assumption that the noise was white and Gaussian? That can help us now. If the noise is white and Gaussian, then the noise components contained in each observation are jointly independent of one another. An important property of independent random variables is that when you sum them together, their means and variances sum. So, let's consider another simple case, where instead of taking one sample per symbol interval, you take two, then sum them together. I'll assume for simplicity that the pulse shape is rectangular (not an exponential rise), so the signal component $s$ in each observation $x_1$ and $x_2$ is the same. What is the difference in signal to noise ratio between just a single observation $x_1$ and the sum of two independent ones?

$$
SNR_1 = \frac{s}{\sigma}
$$

$$
SNR_2 = \frac{2s}{\sqrt{2 \sigma}} = \sqrt{2}SNR_1
$$

So, the signal to noise ratio in the combined observation is larger than using only a single sample (under the assumption of equal signal component and equal-variance white Gaussian noise in both samples that we took). This is a basic observation that points out the potential benefits of taking more than one sample per symbol interval and integrating them together (which, for a rectangular pulse, is a matched filter). In general, you want to cover the entire symbol interval with samples so that your receiver "ingests" as much of the transmitted energy for each symbol, thus maximizing the SNR in the combined output. The ratio of symbol energy to the background noise variance $\frac{E_s}{N_0}$ is often used as a figure of merit when evaluating digital communications system performance.

More rigorously, it can be shown that a matched filter has an impulse response that is identical in shape (that is, "matched", with the only subtle exception being that the impulse response is reversed in time) to the pulse shape that the receiver sees (so it weights more strongly samples that have larger signal components). That shape is a function of the transmitted pulse shape as well as any effects induced by the channel or receiver front end, such as bandlimiting or multipath.

To implement this sort of arrangement in practice, you would convolve the stream of samples taken by your ADC with the time-reversed expected pulse shape. This has the effect of calculating the cross-correlation between the pulse shape and the received signal for all possible time offsets. Your implementation is aided by the precise time synchronization that you have available, so you'll know exactly which matched filter output samples correspond to correct sampling instants. The filter outputs at those times are used as the detection statistic $x$ in the theoretical model above.

I referred to threshold selection before, which can be a complicated topic, and there are many different ways that you can choose one, depending upon your system's structure. Selecting a threshold for an on-off-keyed system is complicated by the likely-unknown signal amplitude $s$; other signal constellations, like antipodal signaling (e.g. binary phase shift keying, or BPSK) have a more obvious threshold choice (for BPSK, the best threshold is zero for equally-likely data).

One simple implementation of a threshold selector for OOK might calculate the mean of many observations. Assuming that zeros and ones are equally likely, the expected value of the resulting random variable is half of the signal amplitude, which is the threshold that you seek. Performing this operation over a sliding window can allow you to be somewhat adaptive to varying background conditions.

Note that this is only intended to be a high-level introduction to the issues inherent in digital communications with respect to detection theory. It can be a very complicated topic, with a lot of statistics involved; I tried to make it somewhat easy to understand while keeping true to the underlying theory. For a better explanation, go get a good textbook, like Sklar's.


Recreating matlab second order filter in python I want to create some custom filters in Matlab filterDesigner tool and then apply those filters to a large set of audio files. I want to use python to do the application of the filters to the audios.
I created a Butterworth filter in Matlab and exported it as an ASCII file. In there I get two datas. SOS Matrix and Scale Values. Then I read the file in python and got the SOS matrix and filtered my signal using the result as follows (using scipy.sosfilt)
from scipy import signal

x_sos = signal.sosfilt(sos, x)  # x is my signal

The result I get has the desired shape except for the amplitude being different (I applied the same filter in Matlab and compared the output from the python code). I understand that scale values needs to come into this equation. But how do I incorporate that? I don't understand this answer. This question has answers explaining the scale values but still I have no idea how I should use that in the python code

EDIT: (After getting the answer) To anyone who is interested in the final code, here it is. And there's an explanation again in this thread
from scipy import signal
import numpy as np

x_sos = signal.sosfilt(sos, x) * np.prod(scale)  # x is my signal

The scale values should be all $1$ except for the first Unless that's $1$ too. Multiply your output signal with that first scale value OR multiply the first three numbers in the first three row of you SOS matrix with that number before filtering.
EDIT
Turns out the filter designer makes this more convoluted than it needs to be. Just multiply all the scale values and multiply output with the product.


Construct Postfilter for overall Minimum or Linear Phase Filter I have a given FIR filter impulse response $h[k]$ for a linear system $y[k] = x[k] \star h[k]$. $x$ is some input, $y$ is the output of the system.

My goal is to construct a postfilter $g[k]$, so that we have $y[k] = x[k] \star h[k] \star g[k]$ with $h_{\rm overall}[k] = h[k] \star g[k]$ being minimum phase.

Though $h$ can vary, it will in general always have most of its zeros on the unit circle.

Question now: How can we construct $g[k]$?
My original approach is quite simple:

Let $a_i$ be zero of $h$ outside the unit circle, then form $G_i(z) = \displaystyle\frac{z - 1/a^\ast}{z-a}$
Form the overall $G$ as $G(z) = \displaystyle\prod_i G_i(z)$

In theory this should work, right?

I implemented this approach in MATLAB and the results are completely useless. My assumption would be that there are (if my approach is correct) are of numerical nature. So my questions are:

Is my approach described above correct?
How can I implement this in Matlab the best way?

What you're trying to do is impossible to implement (i.e., if we require that the post-filter is causal and stable). If the original filter $h[n]$ is not already a minimum-phase filter, it has zeros outside the unit circle. These zeros need to be compensated for by the post-filter $g[n]$, i.e., the filter $g[n]$ must have poles outside the unit circle, and, consequently, it cannot be causal and stable.

On the other hand, constructing a linear phase filter using a post-filter is possible. Just make sure that for each zero $z_{0,i}$ of the filter $h[n]$, the post-filter $g[n]$ has a zero at $1/z^*_{0,i}$, i.e., for each zero you need another one mirrored across  the unit circle.


DSP or signal/image/data processing jokes Some other StackExchange/StackOverflow sites are entertained with a certain level of humor or fun.  What is your favorite “data analysis” cartoon? is particularly worth  mentioning (IMHO)  this xkcd cartoon on causality and correlation (DSP folks know what I mean):

So Is humor allowed in answers? I hope so, se for instance:

StackOverflow: What is your best programmer joke?
MathOverflow: Do good maths jokes exist
SE.Stats (Cross Validated): Statistics Jokes
SE.Stats (Cross Validated): What is your favorite “data analysis” cartoon?
SE.Maths Questions about math jokes

Humor in signal processing or image processing does exist, albeit shallow, and narrow-spread. I have encountered it first (I had zero humor before) at the ICASSP 2000 DSP humor exhibition in Istanbul, Turkey. The image below comes from Humor in DSP by EURASIP:

So I believe that we deserve some humor and joke tags and actual lightness.
 Go ahead!

An airplane is leaving Warsaw (the capital of Poland), and gets caught in a terrible winter storm. The plane rolls, pitches and yaws. The crew is expecting the plane to crash or break up at any time. But one DSP student suddenly stands up and asks some passengers on the right aisle to move across the plane. Passengers are reluctant but the student insists. “Hurry, our fate depends on your move.” Most passengers quickly leap to the other side, and suddenly the plane stops shaking. The pilot, who had a DSP degree too,  enjoys calming everybody:

  "Thank you for the Poles who have moved to the left-half of the plane for stability."

Additional references:

The Routh Stability Criterion
Warsaw, Poland
Stability Revisited

Another visual pun (dedicated to @MarcusMuller) is this chocolet decomposition:


Deconvolution of a 1D Signal with Known Kernel (Square Wave) I have a signal measured from a radiation detector in a narrow beam of radiation. The peaks I get are quasi-gaussian in shape, see attached picture.

The signal is not a function of time, rather a function of distance. The x-axis is is mm and the y-axis is in arbitrary detector response.

The detector used to measure this signal had a finite width, which is contributing to a broadening of the peaks. What I want to do is deconvolve a square wave of width equal to that of the detector from the signal, hence removing some of the broadening effect.

I am hoping to do this in matlab, however I am having trouble using the deconv function due to each data set being two vectors, x and y, and the fact that each data set is a function of linear distance not time.

Any ideas on how to go about this?

General Solution
We have a Deconvolution problem with known operator.
One way to define the objective function is:
$$ \arg \min_{x} \frac{1}{2} {\left\| h \ast x - y \right\|}_{2}^{2} = \arg \min_{x} f \left( x \right) $$
The are 3 method to solve this:

Use Gradient Descent by deriving the gradient of $ f \left( \cdot \right) $.
Write the problem in Matirx Form and either use direct solver (One could use Toeplitz System Solver) of this Least Squares problem or again use iterative solver in Matrix Form.
Solve this in Fourier Frequency Domain.

I previously solved similar problems using the Matrix Form (See my answer to Deconvolution of 1D Signals Blurred by Gaussian Kernel) so this time we'll solve it in the convolution form.
The Gradient
Clearly the derivative of Inner Product is related to the adjoint of the operator.
In the case of Convolution the Adjopint is the Correlation Operator.
Let's define the Convolution Operation $ \ast $ with MATLAB Code.
So $ y = h \ast x $ will be in MATLAB:
vY = conv(vX, vH, 'valid');

Defining $ f \left( \cdot \right) $ in MATLAB would be:
hObjFun = @(vX) 0.5 * sum((conv(vX, vH, 'valid') - vY) .^ 2);

With the Derivative being:
vG = conv2((conv2(vX, vH, 'valid') - vY), vH(end:-1:1), 'full');

Where when we convolve with the flipped version of a vector it means we're basically doing correlation.
Solution in MATLAB
So, one must pay attention that since we use valid convolution the size of $ y $ and $ x $ are different (Depends on the size of $ h $).
Usually we chose $ h $ to have odd number of elements to have its radius to be defined.
Then the MATLAB solution would be:
vX = zeros(size(vY, 1) + (2 * kernelRadius), 1);
vObjVal(1) = hObjFun(vX);

for ii = 1:numIterations
    vG = conv2((conv2(vX, vH, 'valid') - vY), vH(end:-1:1), 'full');
    vX = vX - (stepSize * vG);

    vObjVal(ii + 1) = hObjFun(vX);
end

Specific Solution
In your case you asked for $ h $ to be Box Blur.
So I did the above with this model. Let's go through results.
First, let's see if the solution we got to is a real solution.
In order to examine that, we will convolve the result of the optimization with $ h $ and compare result to the input data. We expect it to be similar to the input.

As can be seen above it is, indeed, a perfect reconstruction of the input data.
We can see this by the objective value as well:

Yet the result is not good:

So, why is the result so bad?
There are 2 options:

The SNR is not good enough to solve the inverse problem as is. It requires additional regularization like Wiener / Tikhonov (They are the same). I implemented it in the code so you have Wiener Filter to play with and understand its working.
The model of Box Blur doesn't fit.

So in my code I implemented Wiener / Tikhonov Regularization (See paramLambda in my code).
The full code is available on my StackExchange Signal Processing Q51460 GitHub Repository.


Phase wrapping in Bode plot I'm trying to model the transfer function of a linear system (a MOEMS mirror) based on a datasheet plot of its response. I'm developing the model in the python-control model and I think that I've got the model pretty close, except that the Bode plot appears to wrap the phase response to keep it within $\pm 180^\circ$. My question is threefold:

Is there a material difference between my model and what's in the datasheet with respect to the frequency wrapping?
Is there a simple way in the python-control library to unwrap the plot
Are there any other suggestions to improve my model?

Here is the response plots from the datasheet:

Here is my model:

There is not really a difference between wrapped and unwrapped phase, it is mainly more visually pleasing when it is continues. The only way this could occur in an actual transfer function is when you would have two complex conjugate unstable poles and two complex conjugate stable zeroes with the same complex and real part, but the real part with opposite sign and the real part should be very small.

One way you could unwrap the phase is by splitting the transfer function up into multiple terms. Since multiplying two transfer functions, means adding their magnitude in dB and phase in degrees of radians. So you could first calculate the bode plot of the two complex conjugate poles (which should stay within 180°, unless it already wraps -180° to 180°) and then add the bode of the real pole to it. Or in this case you could look when the phase is bitter than 0°, in which case you subtract 360°.

In order to improve you model I would take a slightly better at the pole location, since their phase does not seem to drop instantly, so a small real part. There also might be a small delay ($e^{-Ts}$), since the phase of the data drops further than -270°.


Recovering DTFT from Z-transform The relationship between the Z-transform and DTFT can be expressed like:
$$ H(e^{j \omega}) = H(z)|_{z = e^{j \omega}}$$
Graphically, evaluating the Z-transform on the unit circle is shown as sweeping $\omega$ as $e^{j \omega}$ is plotted, represented as a phasor rotating around the unit circle.
My question is, isn't there a problem in that the argument of an exponential is supposed to be unitless? Any complex number (on the unit circle) can be represented as $e^{j \phi}$ where $\phi$ is the number's phase and is unitless (radians). But $\omega$ has units of frequency, so what's happening here? Is there implicitly a time unit of 1?
I've seen some texts cancel the units by normalizing the frequency, so instead of following $e^{j \omega}$ you follow $e^{j 2\pi (f / f_s)}$, where $f_s$ is a sampling frequency. How would assuming an arbitrary $f_s$ affect the calculated frequency response?

The $\omega$ in the frequency response of a discrete-time system $H(e^{j\omega})$ is indeed unitless. The frequency response $H(e^{j\omega})$ is periodic with period $2\pi$. If the discrete signal is obtained via sampling with sampling rate $f_s$, then the relation between $\omega$, the actual frequency $f$, and the sampling frequency $f_s$ is
$$\omega=2\pi\frac{f}{f_s}\tag{1}$$
Consequently, as a function of $f$, $H(e^{j2\pi f/f_s})$ is periodic with the sampling frequency $f_s$.


What preprocessing technique can help me to filter the background text noises using matlab...the image is shown below? Can anyone help me to improve the degradation of the following image....i am working on OCR system in Amharic language and what is the best image filtering technique that can help me to separate text data from background noise?

This problem is called (ink) bleeding. A popular approach is to separate or segment the document image using graphical models. Relevant papers include:

Document Ink Bleed-Through Removal with Two Hidden Markov Random Fields and a Single Observation Field
User-Assisted Ink-Bleed Reduction
Bleed-through removal in degraded documents
Directed Assistance for Ink-Bleed Reduction in Old Documents
Ink-Bleed Reduction Using Functional Minimization


Creating a noise model with a fixed PSNR My goal is to create a noise model so that the PSNR of the noisy image compared to the original is some pre-defined constant $\alpha$. This means the mean squared-error will also be a constant.

My model is as follows:

The additive noise is a random normally distributed variable with $\mu = 0$ and unknown $\sigma$.

$K(i,j) = I(i,j) + f(x; 0, \sigma)$

So,

$$ MSE = \frac{1}{mn}\sum_{i = 0}^{m - 1}\sum_{j = 0}^{n - 1}\ [K(i,\ j) - I(i,\ j)]^2 = \frac{1}{mn}\sum_{i = 0}^{m - 1}\sum_{j = 0}^{n - 1}\ [\ f(x;\ 0,\  \sigma)]^2$$

where $$f(x;\ 0,\ \sigma) = \frac{1}{\sqrt{2\pi}\sigma}e^{\frac{-x^2}{2\sigma^2}}$$
i.e. the Gaussian pdf with unknown variance.

What I'm confused about is how to solve for sigma when the values of the pdf are supposed to be independent samples from the distribution. If I'm understanding correctly, this means that the MSE itself becomes a random variable and cannot be a constant.

If we were to make it a constant, wouldn't the values of $f(x)$ for each pixel have to already be known for the equation to be solvable? In which case the distribution would need to be already known to randomly sample the distribution to produce the values of $f(x)$...

Anyhow, if anyone has advice on how to construct noisy images with a pre-defined PSNR, that would be appreciated.

From wikipedia:

$$
    PSNR = 20 \cdot \log_{10}(MAX_I) - 10 \cdot \log_{10}(MSE)
$$

So, you could:

To identify the max value of your image: $MAX_I$
Generate noise samples for $\sigma = 1$
Solve for $MSE$ in the equation above.
Scale all you generated noise samples by the same factor to achieve that $MSE$ value.


PCA to reduce dimensionality to 99% variance I'm attempting to use PCA to reduce the dimensionality of a dataset I have. I want to explain 99% of the variance in the dataset, and I think I've been able to determine that, but I'm unsure what I have to do to my original dataset to reduce the dimensions. Here's what I've got so far (MATLAB code):

rng 'default'
M = 1000; % Number of observations
N = 500; % Number of variables observed
X = rand(M,N);

% Remove mean.
X = (X - repmat(mean(X),[size(X,1) 1]));

% Determine the eigenvector (V) and eigenvalues (D) of the covariance of
% the matrix X.
[V, D] = eig(cov(X));

% Calculate variance of the vectors (largest to smallest).
var_vec = flipud(diag(D));

% Calculate the cumulative percentage of the variances.
percentages = cumsum(var_vec) / sum(diag(D));

% Find the 99% variance
PC99 = find(percentages >= 0.99, 1);

PC99 tells me the that components 1 to 464 (out of 500) contain 99% of the variance, which I believe to be correct, but I'm now unsure on how to manipulate my original data in order to reduce the dimension from 1000x500 to 1000x464.

Thanks for any help.

If what you want is PC compression, a more useful form of PCA would be achieved via Singular Value Decomposition (which is, in many cases, more accurate and faster than an eigendecomposition).

Assuming that X is a "tall" matrix with dimensions [n, p], n>p;
Then:

 [U, S, V] = svd(X);

eig_vals = diag(S.^2); %recall that eigenvalues of cov(X) are the singular values of X, squared.
percentages = cumsum(eig_vals)./sum(eig_vals); % cumulative variance.

% Find the 99% variance
PC99 = find(percentages >= 0.99, 1);

X_reduced = U(:,PC99)*S(PC99,PC99)*V(PC99,PC99)'; % compressed version of X

If X is a wide matrix [n < p], you have to adjust the above expression a bit (all of V will be used, but only part of U, unlike here where V is truncated as well). We leave it as an exercise to the interested reader!

I will leave wiser people to comment on the judiciousness of PCA for compression. In many cases there are much more efficient schemes, such as wavelet decomposition followed by "thresholding" (zeroing out the coefficients  that fall below a certain threshold). But since your question assumes that you know what you want, so does the answer.


Moving from deterministic signals to stochastic signals in s-domain (Power Spectral Density) Assume we have the following system (coming from control systems theory, hence in s-domain)
$ Y(s) = H_A (s) \cdot A(s) - H_B (s) \cdot B(s) $
I now wish to consider $a(t)$ and $b(t)$ as white noise of unit variance, and I'm interested in the Power spectral density of $y(t)$ (rather the RMS of y(t) derived via the integral of the PSD of $y(t)$, but regardless).
Intuition tells me, that I should get something along the lines of
$ |Y(j\omega)|^2 = |H_A (j\omega)|^2 \cdot 1 + |H_B (j\omega)|^2 \cdot 1 $
But I cannot show how. Especially the switch from subtraction to addition stumps me.

You have to look at the autocorrelation function of $y(t)$:
$$R_y(\tau)=E\{y(t)y(t+\tau)\}\tag{1}$$
with
$$y(t)=(h_A\star a)(t) - (h_B\star b)(t)\tag{2}$$
where $\star$ denotes convolution. If you write out $(2)$ with integrals and plug it into $(1)$ then, with the given assumptions on $a(t)$ and $b(t)$, you'll see that the mixed terms with the negative sign are zero, and the two remaining terms with a positive sign become
$$R_y(\tau)=\sigma_a^2r_{h_A}(\tau)+\sigma_b^2r_{h_B}(\tau)\tag{3}$$
where $\sigma_a^2$ and $\sigma_b^2$ are the variances of the processes $a(t)$ and $b(t)$, respectively, and $r_{h_A}(\tau)$ and $r_{h_B}(\tau)$ are the deterministic autocorrelation functions of $h_A(t)$ and $h_B(t)$, respectively. Finally, taking the Fourier transform of $(3)$ gives
$$S_y(\omega)=\sigma_a^2|H_A(j\omega)|^2+\sigma_b^2|H_B(j\omega)|^2\tag{4}$$
for the power spectrum of $y(t)$.


How to get more accurate results with OpenCV I would like to ask about this function  in OpenCV:

cvTermCriteria( CV_TERMCRIT_ITER | CV_TERMCRIT_EPS, 20, .3 )

How do I get more accurate results - do I need to increase the value of epsilon or decrease it?

This termination criteria tells the algorithm to stop when it has either done
20 iterations or when epsilon is greater than .3.  You can play with these parameters for speed vs accuracy but these values work pretty well in many situations.

So if you want to improve accuracy then you will get more accurate value with more iterations but the time taken for doing such iterations will be improved , I mean speed will be reduced.


transfer function of wiener filter I am trying to understand the function of a Wiener filter. I get that the Wiener filter minimizes the mean square error between the estimated random process and the desired process.

Does this mean that when you have a signal with low and high frequencies and the noise is only in in the high frequencies that the resulting signal only contains low frequencies? What would the transfer function look like in this case?

In general, an ideal Wiener filter (i.e. a non-causal filter with an infinitely long impulse response) has the following frequency response

$$W(\omega)=\frac{S_{dx}(\omega)}{S_x(\omega)}\tag{1}$$

where $S_{dx}(\omega)$ is the cross-power spectral density of the desired signal $d[n]$ and the input of the Wiener filter $x[n]$, and $S_x(\omega)$ is the power spectral density (PSD) of the input process $x[n]$. If the input is just a noisy version of the desired signal, and if the desired signal and the noise are uncorrelated, then we have

$$S_x(\omega)=S_d(\omega)+S_n(\omega)\tag{2}$$

where $S_n(\omega)$ is the PSD of the additive noise, and $S_d(\omega)$ is the PSD of the desired signal. Furthermore, we have

$$S_{dx}(\omega)=S_d(\omega)\tag{3}$$

because the noise and the desired signal are uncorrelated. Plugging (2) and (3) into (1) gives

$$W(\omega)=\frac{S_d(\omega)}{S_d(\omega)+S_n(\omega)}\tag{4}$$

This expression is very easily interpreted: for frequencies where the noise is much stronger than the signal (i.e. $S_n(\omega)\gg S_d(\omega)$ the Wiener filter strongly attenuates the input signal, whereas for frequencies with a very good signal-to-noise ratio (SNR) (i.e. $S_d(\omega)\gg S_n(\omega)$) the filter passes the input with a gain of approximately $1$.

So your guess is correct: if there's a lot of noise in the high frequencies (relative to the desired signal), then these frequencies are attenuated, and the low frequencies, where the SNR is high, will pass the filter with a gain of approximately $1$.


Recognition of household objects ( image processing, pattern recognition ) using wide lens? Camera setup: We have setup some stereo cameras in a living apartment. That is, the indoor environment is monitored. With the stereo cameras, wide lens (3.5mm) are used to cover a big volume.

The height from the floor is around 2.8 meters. The objects (such as mug, bottle, telephone) are at least 3 meters away. For instance, an object (75mmx102mm), which located 3 meters away from the camera, is represented by 15x20px in the camera image. Thus, the images are getting smaller in the far field.

I have around 30 different objects to be recognized.

I do not use the depth information, because it is not so accurate. I just use RGB values from a single camera. Image resolution is 1360 x 1024 pixels.

Approaches: 1. Point detectors/descriptors, matches (some models in the database, and check the match one by one) 2. Bag of visual worlds + SVM classification (5 object categories)

I had experiences with Haar-cascade but I never tried for my current issue.

What methods/approaches should I try to investigate?

Thank you in advance,

I share your assumption about depth being useless here.

The approach #1 based on point detectors seems also useless, because there is probably a difference in scale between your reference (learning) images and the representation of objects in the real pictures that is so big that point detection becomes useless. Information is not present at the same level.

So, you're left with approach #2, which sounds reasonable since it's also more-or-less the state-of-the-art in object detection benchmarks.
You may improve the performance of your system by using a training set of pictures that have the size of the objects in the real (test) images.

Also, you will use only the descriptor part of the existing keypoints and apply it in sliding windows. Since oyu have a very poor resolution on objects, I believe that SIFT-like descriptors (e.g., HoG) will have poor performance, while binary descriptors (LBPs, BRISK, FREAK, ORB) will probably behave better (they are less dependant on pixel accurate resolution).

Also, if your objects are very different, you may try a simple image correlation (using again learning images downscaled to the actual size of the objects in the test pictures).


Showing frequency and amplitude after doing an FFT on a signal I have executed an FFT on a .wav file of a person saying "ahhh" which shows the frequency and amplitude of the audio file.  When I run it on a .wav file it shows the frequency 0 Hz as having the highest amplitude.  I'm not sure where the problem is within my code and how I can fix it.

How does 0 Hz have the max amplitude?  The reason I'm trying to fix this is that I want to find the frequency with the highest amplitude in the wave file then run some calculation against it but running calculations against 0hz doesn't make any sense to me especially when I have to do multiplication with it.

When I look at the array_sort variable it says the frequency is 0hz and 0.117341 as it's amplitude.

Here's the wave file and .m file I'm using
http://dl.dropbox.com/u/6576402/questions/fft/11262012_44100ahh.wav

http://dl.dropbox.com/u/6576402/questions/fft/test_rtfftphase_question.m

clear all, clc,clf,tic
dirpathtmp=strcat('/tmp/');
[vp_sig_orig, fs_rate, nbitsraw] = wavread(strcat(dirpathtmp,'/nov26/11262012_44100ahh.wav')); %must be mono
fs=fs_rate;

t_rebuilt=linspace(0,2*pi,fs); %creates same size time for new signal as orginal signal good for error checking

vp_sig_len=length(vp_sig_orig); %get sample rate from vp fs_rate needs to be an even number?

% Use next highest power of 2 greater than or equal to length(x) to calculate FFT.
nfft= 2^(nextpow2(length(vp_sig_orig)));
% Take fft, padding with zeros so that length(fftx) is equal to nfft
fftx = fft(vp_sig_orig,nfft);
sigfft= fft(vp_sig_orig);
sigifft=ifft(sigfft);
sigphase = unwrap(angle(sigfft')); %get phase of orginal signal
% Calculate the number of unique points
NumUniquePts = ceil((nfft+1)/2);
% FFT is symmetric, throw away second half
fftx = fftx(1:NumUniquePts);
% Take the magnitude of fft of x and scale the fft so that it is not a function of the length of x
mx = abs(fftx)/length(vp_sig_orig); %replaced for testing from stackexchage

if rem(nfft, 2) % odd nfft excludes Nyquist point
  mx(2:end) = mx(2:end)*2;
else
  mx(2:end -1) = mx(2:end -1)*2;
end

amp=mx;
ampinv=abs(amp-max(amp));

% This is an evenly spaced frequency vector with NumUniquePts points.
freq_vect = (0:NumUniquePts-1)*vp_sig_len/nfft;
freq=freq_vect';

%get phase of new signal
phase = unwrap(angle(fftx)); %get phase of orginal signal
array=[freq,amp];
array_sort=sortrows(array,-2); %sort by largest amplitude first

I'm using Octave 3.2.4 on Linux Ubuntu 64bit.

A strong amplitude response at 0 Hz simply means that you have a very strong DC offset. In other words, it just means that the mean of your signal is not 0.

If this is the only problem you have, then all you really need to do, is remove the mean of your signal. In other words:

vp_sig_orig = vp_sig_orig - mean(vp_sig_orig);


Why the curvature of straight line in curvature scale space is not zero? I have read the paper "Shape similarity retrieval under affine transforms", and use c++ to implement it.The step as follow;

Calculate the first and second Gaussian Derivatives

Convolving each coordinate of the curve with a Gaussian function,and calculate the first and second Derivatives of x,y

Compute curvature of the curve

The paper main method as follow

Question:I compute the curvature of a rectangle as sigma is 0.66667. But the result what i got is strange. The curvature of four lines of rectangle are different and not zero? The curvature are 1.9, -6.3, -5.6, 1.5. Then, i roatate rectangle,the curvature of every point is different. Why? Can somebody help me ?Thank you very much

My guess is you're calculating the derivatives wrong, or you're applying the whole thing to the wrong input data. I can show you the expected output for each step graphically, so you can compare it to what you get. I've implemented this in Mathematica; I'll post the code at the end, but I don't know if it'll help you much - Mathematica has built in functions for e.g. convolutions, that you would have to write yourself in C++.

First the input data for a rectangle would look like this:

That's 150 points along a rectangular path, using linear interpolation between the corner points.

Next, the Gaussian and it's derivatives (using $\sigma=3$ for illustration):

If you convolve these with the X and Y coordinates separately, you get this:

The first row is just for illustration - you don't actually need the convolution with the Gaussian itself, only the two derivatives.

At the sides of the rectangle, the first derivatives (second row) are constant and the second derivatives (third row) are 0. (You might want to write the results you get to a text file and look at them e.g. using Excel to check this.)

Then, using the curvature formula you quoted, you get this:

Along the straight sides of the rectangle, the curvature is 0 as expected. At the corners, the curvature of a rectangle is infinite; however, what we're calculating here isn't the real curvature, but the curvature of a Gaussian-smoothed rectangle, with filter size $\sigma=3$:

To check this result, I'll draw a circle with radius $r=\frac{1}{\kappa }$ at each point:

Mathematica code:

(* interpolate between corner points *)
corners = {{-2, -1}, {2, -1}, {2, 1}, {-2, 1}, {-2, -1}};

dist = Rescale@Prepend[Accumulate[Norm /@ Differences[corners]], 0];

int = Table[
   Interpolation[Transpose[{dist, corners[[All, i]]}],
    InterpolationOrder -> 1], {i, 1, 2}];

pts = Outer[#1[#2] &, int, Range[0., 1 - 10^-9, 1/150]]\[Transpose];

cornerIdx = Most[dist*Length[pts] + 1];

Row[{
  ListLinePlot[pts, AspectRatio -> Automatic, Mesh -> All,
   ImageSize -> 400, PlotLabel -> "Points",
   Epilog -> {Red, PointSize[Large], Point[pts[[cornerIdx]]]}],
  Spacer[50],
  ListLinePlot[Transpose[pts], PlotRange -> All,
   GridLines -> {cornerIdx, {}}, PlotLabel -> "X/Y Coordinates",
   PlotLegends -> {"X", "Y"}, ImageSize -> 400]
  }]

(* Gaussian & derivatives *)
\[Sigma] = 2;

g0 = Table[
   1/(Sqrt[2 \[Pi]] \[Sigma])*Exp[-x^2/(2 \[Sigma]^2)], {x, -10., 10}];
g1 = Table[-x/(Sqrt[2 \[Pi]] \[Sigma]^3)*
    Exp[-x^2/(2 \[Sigma]^2)], {x, -10., 10}];
g2 = Table[(-\[Sigma]^2 + x^2)/(Sqrt[2 \[Pi]] \[Sigma]^5)*
    Exp[-x^2/(2 \[Sigma]^2)], {x, -10., 10}];

ListLinePlot[{g0, g1, g2}, PlotRange -> All,
 PlotLegends -> {"Gaussian", "1st derivative of Gaussian",
   "2st derivative of Gaussian"}]

(* Convolve x and y with gaussian&derivatives *)
derivatives =
  Outer[ListConvolve[#1, #2, Round[Length[g1]/2] + 1] &, {g0, g1, g2},
    Transpose[pts], 1];

Grid[MapThread[
  ListLinePlot[#1, PlotRange -> All, PlotLabel -> #2,
    GridLines -> {cornerIdx, {}}] & ,
     {derivatives,
    {
    {"X (smoothed)", "Y (smoothed)"},
    {Overscript[X, \[Bullet]],
     Overscript[Y, \[Bullet]]}, {Overscript[X, \[Bullet]\[Bullet]],
     Overscript[Y, \[Bullet]\[Bullet]]}}}, 2]]

{{xs, ys}, {dx1, dy1}, {dx2, dy2}} = derivatives;

(* calculate curvature *)
\[Kappa] = dx1*dy2 - dx2*dy1/((dx1^2 + dy1^2)^(3/2));
ListLinePlot[\[Kappa], PlotRange -> All,
 PlotLabel -> "Curvature \[Kappa]", GridLines -> {cornerIdx, {}}]

ListLinePlot[Transpose[{xs, ys}], Mesh -> All]

(* osculating circle animation *)
animationIdx =
  Select[Range[Length[xs]],
   Abs[\[Kappa][[#]]] > 10^-2 || (Mod[#, 5] == 0) &];

Monitor[frames = Table[Rasterize@Row[
      {
       ListLinePlot[Transpose[{xs, ys}], Mesh -> All,
        ImageSize -> 400,
        PlotLabel -> "Osculating circle r=1/\[Kappa]",
        Epilog -> Module[{r, pt, n},
          (
           pt = {xs[[i]], ys[[i]]};
           r = Clip[1/\[Kappa][[i]], {-1, 1}*10^4];
           n = Normalize[{dy1[[i]], -dx1[[i]]}];
           {Red, PointSize[Large], Point[pt], Circle[pt - n*r, r]}
           )]
        ],
       ListLinePlot[\[Kappa], PlotRange -> All,
        PlotLabel -> "Curvature \[Kappa]",
        GridLines -> {cornerIdx, {}}, ImageSize -> 400,
        Epilog -> {Red, Line[{{i, 0}, {i, 10}}]}]
       }], {i, animationIdx}], i];


How to determine the type of a digital filter given its expression? I'm trying to understand a variant of the filter presented in http://www.stockspotter.com/files/PredictiveIndicators.pdf in the "Code Listing 1.". Its expression is given as follow:
$$y_n = c_0 x_n + c_1  y_{n-1} + c_{2}  y_{n-2}$$ with:
\begin{align}
c_0&=(1-c_1-c_2) \\
c_1&=2 a \cos\left(2\frac{\pi\sqrt{2}}{T_s}\right) \\
c_2&=-a^2
\end{align}
and:
\begin{align}
a&=\exp\left({-\frac{\pi\sqrt{2}}{T_s}}\right) \\
T_s&=10
\end{align}
I'm wondering how these coefficients were computed. I'm told it might be a 2-pole Butterworth filter, but my attempts at finding the coefficients for such a filter failed (and my searches online give me many different results, sometimes involving $x_{n-1}$ and $x_{n-2}$ whereas they are not used in this expression.
Does this ``SuperSmoother'' filter has another name in the literature ? How could I describe the effect of such a filter (beside being a low-pass filter)?
The paper is cited in the Linux Kernel's packet rate estimation smoothing.

For the records, here is the response of this filter using the notebook provided in the answers, compared with a moving average and an exponential moving average:

As someone remarked on IRC, this filter looks a lot like the "multiple feedback digital low-pass filter" from Chamberlin 1987 book "Musical Applications of Microprocessors", with $F=\frac{\sqrt{2}}{16}$ and $Q=1$:

Seeing that the paper cites the author of the paper as inventor of the "SuperSmoother" filter, and this filter was (supposedly) good for this specific use case, there's no indication this filter is based on anything but the author's inventive force (his fantasy). He does mention it's a "converted analog filter made from capacitors and resistors", and you'd often apply  Bilinear Transform to do a continuous/discrete time transform.
So, you won't find that "SuperSmoother" filter in literature. It seems to be a purpose-specific invention, not derived methodically using mathematical approaches.
Notice that you're very much on scientifically thin ice with the whole article.
There's no reason you'd want to base your analysis smoothing filter on an analog filter design; the author also doesn't explain why you'd want that. That's something that someone not really understanding how to design filters would do if they have but an analog filter to go on – and building analog filters that are as good for any particular discrete analysis purpose as a filter directly designed for the digital purpose in digital design is um, impossible; the optimization goals are in discrete time, not continuous time.
Also, nope, Butterworth is certainly not a minimum-lag filter, it's a maximum-flatness filtering approach.
He makes it sound like minimal-lag filter design is a new thing that he invented. That's nonsense, minimum-phase filters are literally entry level concepts in any discrete system lecture I've had material of.
The whole paper is very handwaving and has premises that are plain wrong¹. I'd put it into the category of pseudoscience – which means that while I really applaud you looking for meaning and scientific sense in there, you won't.

¹ for example

We also know that the
market is fractal; a daily interval chart looks just like a weekly, monthly, or intraday chart.

No, that's nonsense, and that can be easily seen by the trivial autocorrelation. Or by the fact that there is an average weak-periodic component, but a month isn't "naturally" split into 7 equally sized units, and you won't find that in any autocorrelation. Really, this description is magic/superstition.
The paper claims to be based on scientific studies, but fails to cite even a single one – this should really make you think whether the author isn't trying to sell you that he's got a skill he really doesn't have.


Has anyone used Siglib? Answer a question on the FFTr function? I'm hoping someone - anyone can answer a question on the Siglib library call to perform an FFT on data.

There's no examples anywhere I can find that provide a real application using data from an ADC, so that's where my wheels are spinning.
The examples from Numerix have only canned data, generated by the function SDA_SignalGenerate. Specifically, in the FFT example, SDA_SignalGenerate is passed a value of $0.15625$ for frequency. Can anyone provide some insight as to what this represents? The FFT_LENGTH of the FFT function is $512$, but I do not see what the SDA_SignalGenerate frequency of $0.15625$ means, or what it is normalized to.

Thanks for any help anyone can provide.

This is the frequency normalized by the sampling frequency, i.e. in your case the frequency in Hz is

$$f=0.15625\cdot f_s$$

where $f_s$ is the sampling frequency in Hz. This can also be found in the SigLib reference manual (p. 251).


LTI filtering for wide-sense stationary process Why is it that if $U[n]$ is wide-sense stationary and it is convolved with $h[n]$ to produce $W[n]$, the autocorrelation becomes $R_{WW}[n] = R_{UU}[n]*h[n]*h[-n]$?

I know that in general $R_{WW}[n_{1},n_{2}]=R_{UU}[n_{1},n_{2}]*h[n_{1}]*h[n_{2}]$ and that wide-sense stationary means $m_{U}[n] = m_{U}$ along with $R_{UU}[n_{1},n_{2}]=R_{UU}[n_{1}-n_{2},0]$, but I can't get to the above relation from these facts.

It's basically just about using the definitions and doing the math:

$$\begin{align}R_{WW}[n]&=E\{W^*[k]W[k+n]\}\\
&=E\left\{\sum_mU^*[k-m]h^*[m]\sum_lU[k+n-l]h[l]\right\}\\
&=\sum_m\sum_lh^*[m]h[l]E\left\{U^*[k-m]U[k+n-l]\right\}\\
&=\sum_m\sum_lh^*[m]h[l]R_{UU}[n+m-l]\\
&=\sum_mh^*[m]\sum_lh[l]R_{UU}[n+m-l]\\
&=\sum_mh^*[m]\left(h\star R_{UU}\right)[n+m]\\
&=\sum_mh^*[-m]\left(h\star R_{UU}\right)[n-m]\\
&=h^*[-n]\star h[n]\star R_{UU}[n]\end{align}$$

where $*$ denotes complex conjugation, and $\star$ denotes convolution.