The TI DAC8532 uses a Kelvin resistor-string (string of resistors that divide the reference voltage, with individual switches from the nodes to the output), so it is guaranteed to be monotonic. What you see does not sound normal, you must be writing zeros to it over the serial interface somehow. How did you implement the serial protocol?

The low-pass filter you put on the output is the reconstruction filter. It is used to reduce imaged frequency content (repeated frequency content above the Nyquist frequency). It's not always needed. The DAC has a built-in zero-order hold element that attenuates the images somewhat. If you are driving a system with low-pass frequency response, like a DC motor with a heavy mass attached, then you don't need the reconstruction filter.


consecutive filtered blocks with filtfilt

Both the SciPy and MATLAB implementations of filtfilt use a mirrored part of the beginning and end of the signal (roughly with a length equal to the impulse response of the filter) to reduce the effect of initial conditions on the filtered signal (i.e. to reduce ringing the signal is artificially lengthened by concatenating mirrored parts of the signal).

As these mirrored signal parts are probably not equal to the data in the previous and subsequent data blocks, concatenation of the filtered blocks will not result in a continuous signal.

A possible workaround might be to select longer data blocks from the original signal and discard a part from the beginning and end after application of the filter but before concatenation. Note that this might require a FIR filter to work properly although I would expect acceptable results for a IIR filter if the data block can be made long enough.


Reduce the Number of Intensity Levels of a Grayscale Image in MATLAB

I think by number of levels you want the image full scale grey-scale to be divided piece-wise into given number of levels.
For example: -

If number of levels = 2, then you want only two grey-scales in your image i.e. (0 and 128) or (128 and 255) depending upon if you are using floor or ceil within the range
If number of levels = 4, then you want 4 different levels in your image i.e. (0, 64, 124, 192)

Solution -
This operation can be done in a single line

new_img = ceil(img_gray./step)*step;

Here I have divide the entire range into desired number of parts (levels) with size as step and used ceil function to restrict the result to the upper bounds (you can use floor for lower bounds)

Full Code -

img_color = imread('peppers.png');
img_gray = rgb2gray(img_color);
imshow(img_gray);

[rows, cols] = size(img_gray);
noOfDesiredIntensityLevels = 2;
step = ceil(255/(noOfDesiredIntensityLevels - 1));

new_img = ceil(img_gray./step)*step;

figure
imshow(new_img);

% optional code to show the levels
allSteps = 0;
currStep = 0;
while(currStep < 255)
    currStep = currStep + step;
    allSteps = [allSteps currStep];
end

allSteps

Edit 1 - I have included few extra lines which shows you the different levels (steps) in an array. Its just optional but nice to know the exact values of the levels.


Do I have to flip my kernel when performing an FFT-based convolution?

Flipping & dragging is an animative method used in graphical computation of convolution in time (or space) domain. It is the result of an argument manipulation in $h[n-k]$ (or $x[n-k]$) signal animated as a function of n, but drawn on an axis of k, in the convolution sum: $$y[n] = x[n]*h[n]=\sum {x[k]h[n-k] }$$ On the other hand in frequency domain DFT (FFT) based convolution implementations you don't need it as you simply multiply the two DFTs.


Intuitive definition of ergodicity for random signal

When a process or system in which the average, over time, of the physical quantities that describe it coincides with the average calculated over a set of possible states of the process or of the system itself.
I would use an example like the following:
-You would have the similar results after flipping 100 coins as flipping a coin 100 times.
I hope it helps.


it's in textbooks.  you'll find that, for a sharp filter of given sharpness, that IIR filters will cost less than FIR.

given it's IIR, you'll find that the Direct Form I and Direct Form II will have less computational cost than the State-Variable Filter, Lattice, Normalized Ladder, Gold-Rader, or Zölzer, or Harris-Brooking forms (and there are a zillion other forms, i s'pose, but i can't see them being cheaper to implement than the Direct Forms).  and then you can implement many forms (like DF1, DF2, or Lattice) as one-big direct filter section or as a cascade of second-order sections.  but there are numerical tradeoffs doing it the cheap (from a computation time POV) way.

given it's FIR, there are at least 3 different ways of implementing an FIR, all with different cost issues.  if the FIR filter is totally general (but still causal), then the choice is between the simple direct method (the "transversal" FIR form) and what we like to call Fast Convolution (implementing a long FIR using the FFT).  for some FIRs of a greatly restricted class (that is FIRs that look like IIRs), there is a method that is called Truncated IIR (TIIR) that is about as cheap as an IIR, but it's mathematically FIR.

good question, no short answer that is also comprehensive.  need to get O&S and read it.


N orthogonal basis vectors can completely describe a point in N dimensional space.  Same for N * N orthogonal vectors and N * N samples describing an N*N dimensional space.

A DCT matrix is constructed from an orthonormal set of basis vectors.


Speech Compression - In LPC how does the linear predictive filter work on a general level?

LPC voice coders (starting with the old LPC10 standard, which seems to be the one you refer to here) are based on the source-filter model of speech production. Speech can be characterized by the following properties:

The raw sound emitted by the larynx (through vibration of the vocal folds, or just air flowing through it, the vocal folds being opened).
The transfer function of the filter achieved by the articulatory system, further filtering this raw sound.

Early LPC coders (LPC10) adopt the following model of those two steps:

The larynx emits either white noise, characterized by an amplitude $\sigma$ ; or a periodic impulsion train, characterized by an amplitude $\sigma$ and a frequency $f_0$
The transfer function of the articulatory system is of the form $\frac{1}{1 - \sum_k a_k z^{-k}}$, and is thus entirely characterized by the coefficients $a_k$

The principle of early LPC coders (such as LPC10) is thus to estimate these parameters from chunks of the incoming audio signal ; transmit these over the network ; and have a sound generator reproduce a sound from these parameters on the receiver. Observe that in this process, nothing of the original audio samples is actually transmitted. To make a musical analogy, it's like listening to a piano performance, transcribing it, sending over the score, and getting someone to play it on the other end... The result on the other end will be close to the original performance, but only a representation has been sent over.

  does this mean that we keep a running average of last p samples at
  both the encoder and decoder?

No, this is not how it works. On the encoder side, we run a process known as AR (autoregressive) estimation, to estimate the set of coefficients of the AR filter that matches best the spectral envelope of the input signal. This is an attempt at recovering the coefficients of the filtering that has been performed by the articulatory system.

These coefficients are sent over the network (along with pitch, the voice/unvoiced flag, and loudness). The decoder uses those coefficients to filter a synthetic excitation signal, which can be either white noise (unvoiced frame) or a comb of periodic impulsions (voiced frame). These coefficients are used in the following way on the decoder to recover the output signal $y(n)$:

$$y(n) = excitation(n) - \sum_k a_k y(n - k)$$

Observe that since this is an all-pole, IIR filter, the samples linearly combined are previous samples produced by the filter.

  So that at the encoder we only transmit
  data that corresponds to the difference between this average and
  actual signal?

There is no "averaging" and there is no transmission of a difference signal. The decoder is a sound synthesizer with some parameters ; and the encoder searches for the set of parameters for this synthesizer that most closely matches the input signal.

  Why is it an a linear combination of these previous samples?

Other options would be possible indeed here, but the advantage of using an auto-regressive model are the following:

In this model we use an all-pole IIR filter to simulate the articulatory system. This is a good model, since it can capture strong spectral dips and peaks with a small number of coefficients ; and the articulatory system is indeed capable of emphasizing / attenuating narrow bands of frequencies (See formants). If we had used an all-zeros FIR filter, we would have needed many more coefficients to accurately capture the kind of filter responses that the articulatory system can achieve.
On the encoder end, the problem of of estimating the filter coefficients is computationally efficient - it can be done using the Levinson-Durbin recursion on the first few samples of the autocorrelation. More complex linear models (ARMA) or non-linear models are more computationally expensive or intractable.
On the decoder end, the synthesis is very simple indeed - we just need enough memory to keep track of the n previous samples emitted by the decoder.

  Surely an AVERAGE of these previous samples would be a better indication of the next sample?

This is not true. For example, assume the input signal is a sine wave. The prediction $\hat{y}(n) = (2-\omega^2) y(n - 1) - y(n - 2)$ has zero error, while $\hat{y}(n) = \frac{y(n - 1) + y(n - 2)}{2}$ is wrong most of the time (in particular, it does not capture the fact that the sine wave is strictly increasing half of the time, strictly decreasing the rest of the time, so $y(n)$ shouldn't be "between" the previous values). Speech signals are not sinewaves of course - but you can squeeze much more modelling power out of a model with $p$ parameters (order p AR), than out of a model with zero (averaging).

Also, it's worth keeping in mind that there's a mathematical definition of "better" (minimizing the expected value of the squared of the error), which yields a procedure for finding the optimal value of the coefficients.


The library function log of zero results in a NaN (not a real number, -inf).  So the log of a square with a base of zero can't be plotted.


There is a package from octave. I am using debian jessie and this version of octave-video producing errors and can not be installed. Also, there is no debian package with name octave-video.


What is the difference between the sampling frequency of signal and sampling frequency of filter

I sadly don't have the reputation to comment on Fat32's answer, but let me try to answer directly instead: the frequency response of a digital filter is always relative to the processing sampling rate.
So if you design a filter with cutoff at 100Hz and 1kHz sampling rate, then you are really designing for a normalized cutoff frequency $f/f_s$ of 1/10 and whatever the processing rate (ie. the sampling rate of the processed signal), the cutoff will always be 1/10 times the processing rate. If you design a similar filter with cutoff at 200Hz and 2kHz sampling rate, then you will get the same filter.
Now, if you want to filter a signal sampled at 1kHz and you want the cutoff to be 100Hz, then you need to design a filter with a normalized cutoff frequency of 1/10. If you use the same filter on a signal sampled at 2kHz, then the cutoff will be at 200Hz. If you want to keep the cutoff at 100Hz for the 2kHz signal as well, then you need a different filter, this time for a (normalized) frequency of 1/20 instead.
The filter coefficients don't carry any information about the sampling rate, but rather everything is always relative to the sampling rate of the signal you are filtering and this applies to all digital filters, whatever the design method and whether FIR and IIR.


Increasing SNR and Dynamic Range using Oversampling

OP clarified that the question in the comments as follows:

If we ignore any modulation for now and assume that we are receiving
pure tones plus the band limited noise and we try to improve the SNR
in post processing how much improvement can we expect by oversampling
and is there a limit to it? My original question was about this
aspect.

First consider the case of sampled white noise (no band-limiting of the sampled signal other than anti-alias filtering prior to the A/D converter), the improvement in SNR due to oversampling a pure tone would be the processing gain in a matched filter such as a correlator to $e^{j\omega n}$ given as $\sum x[n]e^{-j\omega n}$:
$$PG = 10\log_{10}(N) \text{ dB}$$
This is because as we sum the samples in the correlation process, the signal component of each sample which is correlated from sample to sample will increase at $20\log(N)$ dB while the noise components will increase at $10\log(N)$ dB: The total signal magnitude would be $x_1+x_2+x_3 +\ldots$ while the total noise rms magnitude would be $\sqrt{x_1^2+x_2^2 + x_3^2 + \ldots}$
But with band-limiting, as we increase the sampling rate beyond the Nyquist frequency for the bandwidth of the signal, there will be no further increase in SNR: since the signal is bandlimited (no longer white), adjacent samples will be correlated both in signal and in noise, and as we sum the samples in an attempt to increase SNR, both the signal and noise components of the signal will increase at the same rate ($20\log_{10}(N)$).
Assuming a constant noise out to $f_b$ and then rolling off after that, we would increase SNR after correlation up to the point where $f_s/2 = f_b$ (for a real signal) Increasing the sampling rate beyond that will not increase SNR. Further we have a practical limit for long averaging situations from ADC spurious free dynamic range (SFDR), phase noise from clocks and local oscillators, and similar $1/f$ noise sources where the signal is no longer stationary and increasing the averaging duration further will begin to degrade SNR.  The Allan Deviation (ADEV) is an excellent statistical tool for determining the optimum duration to average given we have no other constraints on the time interval. (For more on ADEV for this application see What determines the accuracy of the phase result in a DFT bin? )

If the input signal has a fixed analog bandwidth, and you are properly not limited by the quantization noise, and front-end filtering has been sufficient to not have aliasing to occur, then oversampling will not have any effect on the SNR in band.
The idea of oversampling is to spread the quantization noise out so that it is no longer swamping out the noise floor of your signal which should be the limiting factor to SNR, and to allow for appropriate filter design to eliminate noise out of band from folding in at a lower sampling rate. These are the two most prominent reasons to oversample: increasing the effective resolution by spreading out the same total quantization noise over the full sampling bandwidth, and increasing the total frequency space in which we can do digital radio processing of signals outside our primary band of interest, for filter rejection and multi-carrier operations. For more details on the first and it's limitations please see:  What are advantages of having higher sampling rate of a signal?
It is assumed that the proper receiver design will filter out all noise that is outside of your signal bandwidth (the matched filter), so given you are sampling sufficiently high according to Nyquist (2x the bandwidth, plus some margin for realistic filter realization) then there is no reason to sample any higher once the signal is properly filtered (and this only increases power dissipation, increased resources, etc).
The strategy should be to sample sufficiently high at the front-end of the digital receiver for the filter designs involved and in consideration with the analog filtering that may be done for anti-aliasing, then once out of band noise, interference and other channels are filtered out, reduce the sampling rate to be as small as possible (typically 2 samples/symbol).
This may add a further intuitive explanation: Consider a band-limited random noise process with the sampled time domain signal shown below along with a histogram of those samples along the right axis, and then immediately below that is the spectrum. As long as we sample that spectrum more than twice the bandwidth shown, no matter how much faster we sample we will still get the same histogram, from which we get the standard deviation of that noise.
If we had additional noise energy outside the signal's bandwidth, then oversampling can of course help us to eliminate that interference through filtering where it could otherwise fold into band with a lower sampling rate-- but this isn't improving SNR by oversampling, this would be implementing proper receiver design by filtering out interference (which that task alone may require a higher sampling rate).


Perform frequency analysis on grouped pulses

The first step is to low pass filter the signal with 100Hz cutoff and then resample everything at 200Hz. This doesn't change the spectrum in the 0-100Hz range. Therefore, it's not possible to gain anything of use from the faster signals.

The 0-100Hz spectrum will be the same as if there were a single 45-450μs pulse occuring at 10Hz. There will peaks at multiples of 10Hz (harmonics) and the maximum will be the 10Hz peak itself.


Ideal low pass filter output at given sampling frequency

That $\frac{1}{T_s}$ factor in your sampled signal expression is $\frac{20}{\pi}$.

$$\omega_s = 40 \ rad/sec$$
$$2\pi f_s = 40$$
$$f_s = \frac{40}{2\pi} = \frac{20}{\pi}$$

The sampled signal is basically given by :
$$X_{sampled}(f) = f_s \sum^{\infty}_{k=-\infty} X(f - kf_s)$$
Meaning : Sampling a signal with impulse train at sampling frequency $f_s$ gives you, images of $X(f)$ scaled by $f_s$ and centered at $kf_s$.

You can check this answer to understand why : Sampling with Impulse Train explained pictorially

The original signal is at frequency $\frac{15}{\pi}$ and being sampled at sampling frequency $\frac{20}{\pi}$. And hence the frequency components in the sampled signal will be : $\pm k\frac{20}{\pi} \pm \frac{15}{\pi}$.

Notice that for $k=0$, you get the original signal and for $k=\pm 1$, you have an alias term $\cos{(2\pi (\frac{20}{\pi} - \frac{15}{\pi})t)} = cos(10t)$ which will get inside the lowpass filter cutoff frequency.

And, that will give you the final result as :
$$f_s (\cos{(2\pi \frac{15}{\pi}t)} + \cos{(2\pi \frac{5}{\pi}t)})$$
$$ = \frac{20}{\pi} (\cos(30t) + \cos(10t))$$


Am I understanding the phase angle from FFT correctly?

If you're interested in the time delay between the two signals, it's probably most useful to compute the cross-correlation function of the two signals and find its maximum. The location of the maximum will give you the time lag.


How does companding compress dynamic range?

Companding is used when you want a signal to noise ratio (SNR) that's largely independent of the signal level.
Quantization noise is constant, so if the signal level varies, so does the SNR.
Let's consider an audio signal that's quantized at 8 bit. The quantization noise sits at -53dBFS. If you have a loud signal, i.e. -10dBFS you have a healthy 43 dB SNR. However, for a soft signal at -40dBFS, the SNR is very poor at 13 dB.
If you compand you can get an SNR of maybe 30dB for both cases, which is much better. The audibility of the quantization noise is NOT a function of the absolute noise level but the noise level relative to the signal. You can tolerate more noise, if the signal itself is louder since the louder signal masks more of the noise.


When listening in to an AM signals of various frequencies, how do we exactly tune in?

Ok, I think here is a good point to introduce equivalent baseband, because you implicitly are already using it!

So, what does your cosine 5Hz audio signal look like, if you were to set $f_c=0$, in spectrum?

Exactly, one dirac at +5 Hz, and one at -5 Hz! Hence, when you mix that up to $f_c\ne 0$, you get symmetric sidebands. That works for any real signals (ie. for any audio signal) – and even when you add them. That's where your two sidebands come from – positive and (hermitian) symmetrical negative components of any real-valued signal.

So, in baseband, your sum of 5 and 30 Hz cosines have four spectral components – at -30, -5, +5 and +30 Hz, and mixed up to the carrier frequency you get the same discrete spectral components, but added to $\pm f_c$. (your figure also shows a component at $f_c$, but that's only there if you got a DC component in your baseband signal, and/or you're not suppressing the carrier)

Now, there's very different receiver architectures for AM-modulated audio. None of them "just looks every 1/(audio period)". "Looking every so often" presumes you're building something digital. Most AM receivers (most AM being pretty obsolete by now) are not digital. The simplest detector methods really are just rectifying diodes and a low pass filter to get rid of the RF content – I'll leave googling for "AM diode detector" up to you; there's a plethora of good articles on that out there. This is all continuous-time, so there's no "looking that often" there – it all happens by electronically processing the continuous signal.

Now, assuming we're really aiming for digital here:

That 1/(audio period) wouldn't even make sense for audio – it breaks Nyquist; you need at least twice the audio bandwidth as sampling rate to be able to reconstruct the signal.

What one can build is simply a mixer with $f_c$, which will first (continuous-time!) multiply with a harmonic of that frequency (effectively, a complex sinusoid), and then sample (that's the act of looking only every so often) that, and then you get a digital signal that's basically nothing but your original audio sum.

Note that you've been asking about "how to deal with the fact that you add two sines of different frequencies", but that question is totally independent from the AM aspect of that: you just can. Under the Fourier Transform, i.e. in the spectrum, these are orthogonal, i.e. you can perfectly separate them, both in a computer, in a circuit, or with your ear. That's why you're hopefully able to listen to sounds that aren't made of a single tone. In fact, in reality, single-tone sounds are extremely rare.


DSP Hardware for Detecting an Audio "Click"

Probably the easiest (and most accurate) way to do this is the matched filter approach suggested above.  Basically, you can record each click independently and use those as matched filters against the microphone signal.  If the clicks are always the same, click A will ring up very strong on the A matched filter.  Same with B and C.

To implement the matched filter, you can get a development kit with a processor or DSP chip on it.  It will take some coding to implement the solution (probably C).  In college, I worked on the "TMS320C6713 DSP Starter Kit (DSK)" but you can try Arduino or Beagleboard.  Your system will probably not be highly complex.  As for the microphone, you can get a relatively simple one that connects to the 3.5mm jack on the development kit.  (http://www.amazon.com/Pyle-Pro-PDMIC58-Professional-Handheld-Microphone/dp/B003GEBGA0/ref=sr_1_4?s=electronics&ie=UTF8&qid=1440204270&sr=1-4&keywords=microphone+3.5mm)

A matched filter works by multiplying the incoming signal by the expected stream of data.

Example:
If click A lasts 1mS and your sampling rate is 44.1kHz, you would have a 44 sample matched filter (the click may be much shorter).  You can store this in the code and point multiply your incoming signal by this matched filter.

MATLAB example code:
https://stackoverflow.com/a/19471009/3920284

The maximum value you get when you multiply by your input stream will correspond to the time when your click happened.  You can set a threshold like the example does.

One thing that may change over time is your clicks' sound: as the hardware ages, the clicks will sound different since wearing out may occur.


How to Remove the Periodic Oscillations from a Signal

Their are some problems in the code and same algorithmic issues:
Code:

You are estimating the wrong frequency: f2=1/((365*24*60)/2); should be f2=1/((365*24*60))/2;

Algorithm:
Basically you are using a fourier series (http://en.wikipedia.org/wiki/Fourier_series) with those terms corresponding to the known frequencies. But you should be aware of the limitations of that approach:

You have to estimate the dc part (frequency 0), because your signals' dc value is not zero!
In theory this approach works only for full periods of each frequency. That's why f1 is estimated quite good, but f2 is not. To work arround that you need either more data, because the more periods are in the data the smaller will be the estimation error. And you should use a windowing function, to get rid of some side effects by cutting the signal at the beginnig and the end.

The improved code could look like this:

clear
doWindow = true;
scaleFreq = 1;
%% Creating the data

g=2.5;
t=linspace(0,g,g*365*24*60);              %n-years of 1min data
o=rand(size(t));                         %random noise
%we create two sinusoids
p1=2*sin(2*pi*t*scaleFreq);
p2=1.5*cos(1*pi*t*scaleFreq);
p=p1+p2;
x=p+o;                                     %and add noise to them

%% The method

N=length(x);            %no of data
n=(1:N);
T=g*365*24*60;          %time that the data covers
f1=scaleFreq*1/(365*24*60);       %frequency of 1 yr oscillations
f2=scaleFreq*1/((365*24*60))/2;   %frequency of 1/2 year oscillations

a1=f1*T;
a2=f2*T;

if(doWindow)
    window = hanning(numel(x)).';
else
    window = ones(size(x));
end
xWin = x.*window;%windowing
xMean = mean(xWin);%calculate dc
xWin = xWin-xMean;%subtract dc

buildFunkCos1 = cos((2*pi*a1*n)/N);
buildFunkCos2 = cos((2*pi*a2*n)/N);
buildFunkSin1 = sin((2*pi*a1*n)/N);
buildFunkSin2 = sin((2*pi*a2*n)/N);

A1=(1/N)*sum(xWin.*buildFunkCos1)
A2=(1/N)*sum(xWin.*buildFunkCos2)
B1=(1/N)*sum(xWin.*buildFunkSin1)
B2=(1/N)*sum(xWin.*buildFunkSin2)

if(doWindow)%do correction due to windowing
    A1 = 2*A1;
    A2 = 2*A2;
    B1 = 2*B1;
    B2 = 2*B2;
end

C1=2*sqrt(A1^2+B1^2)%amplitude
C2=2*sqrt(A2^2+B2^2)

fi1=atan(B1/A1)%phase shift
if(A1<0)
    fi1 = fi1 + pi;
end
fi2=atan(B2/A2)
if(A2<0)
    fi2 = fi2 + pi;
end

%the reconstructed set of data for these two frequencies
v1=C1*cos(2*pi*t*scaleFreq-fi1);
v2=(C2*cos(1*pi*t*scaleFreq-fi2));
v=v1+v2+xMean;%sum estimates (freq 1,freq 2, dc)
r=x-v;    %I wasn't sure what was meant by removing, so i just subtracted

%% Visualization
figure(1)                                  %now lets plot the data
plot(t,[p1;p2;p]);
legend('p1','p2','p');
title('single sine signals and sum');

figure(2)                           %and, original signal and noise separately
hold off
plot(t,[v1;v2;v]);
legend('v1','v2','v');
title('estimated sine signals and sum');

figure(3)                    %and once again data (blue), reconstructed signal
hold off
plot(t,x,'b-')            %(green), and the final result afer subtraction
hold on                   %(red)
plot(t,v1+v2,'g--')
plot(t,r,'r:')
xlabel('t[min]')
ylabel('w/e')
legend('measured','reconstructed','result');

figure(4)
hold off
plot(t,[buildFunkCos1;buildFunkCos2;buildFunkSin1;buildFunkSin2]),
legend('buildFunkCos1','buildFunkCos2','buildFunkSin1','buildFunkSin2');

Here you can switch the windowing and also test what happens if more periods are in the signal. Try, scaleFrequ=10 und doWindow=true and the result should be what you are looking for.


Finding the linear prediction coefficients for a sampled sinusoid

In the absence of more details, I assume that

the phase $\Theta$ is a random variable and
$\Theta$ is uniformly distributed on the interval $[0,2\pi]$, so that its probability density function is $\frac{1}{2\pi}$ on that interval.

This is a common choice of model for such problems, and I will show that its guarantees that the stochastic process $x$ is wide-sense stationary (WSS).

\begin{equation}
\begin{split}
\mathsf{E}[x[n]] &=~ \mathsf{E}[2\sin(0.25\pi n + \Theta)]\\
&=~ 2\mathsf{E}[\sin(0.25\pi n + \Theta)]\\
&=~ 2\int_{0}^{2\pi}\sin(0.25\pi n + \theta)\frac{1}{2\pi}d\theta\\
&=~ \frac{1}{\pi}\int_{0}^{2\pi}\sin(0.25\pi n + \theta)d\theta
\end{split}
\end{equation}
The integral is the integral of sine over a $2\pi$-wide interval, so it is equal to 0. Hence $\mathsf{E}[x[n]] = 0$, independent of $n$.

\begin{equation}
\begin{split}
\mathsf{Var}[x[n]] &=~ \mathsf{E}\left[\left(x[n] - \mathsf{E}[x[n]]\right)^2\right]\\
&=~ \mathsf{E}\left[\left(x[n] - 0\right)^2\right]\\
&=~ \mathsf{E}\left[x[n]^2\right]\\
&=~ \mathsf{E}\left[4\sin^2(0.25\pi n + \Theta)\right]\\
&=~ 4\int_{0}^{2\pi}\sin^2(0.25\pi n + \theta)\frac{1}{2\pi}d\theta\\
&=~ \frac{2}{\pi}\int_{0}^{2\pi}\sin^2(0.25\pi n + \theta)d\theta\\
&=~ \frac{2}{\pi}\times\pi ~=~ 2,
\end{split}
\end{equation}
which is independent of $n$. We see that $x$ is a WSS process.

This is important because we rely on expected values of the form
\begin{equation}
\begin{split}
\mathsf{E}[x[n]x[m]] &=~ \mathsf{E}[4\sin(0.25\pi n + \Theta)\sin(0.25\pi m + \Theta)]\\
&=~ \frac{4}{2\pi}\int_{0}^{2\pi}\sin(0.25\pi n + \theta)\sin(0.25\pi m + \theta)d\theta
\end{split}
\end{equation}
I leave it as an exercise to show that the value of this integral depends only on the difference $n-m$ and not on $n$ and $m$ separately.

In particular, we can define $R[i] = \mathsf{E}[x[n]x[n-i]]$ and know that $R[i]$ is really a function of $i$ alone and has no dependence on $n$.

For second-order linear predictive coding (LPC), we estimate $x[n]$ from the previous two values:
\begin{equation}
\widehat{x}[n] = a_1x[n-1] + a_2x[n-2].
\end{equation}
The error is then
\begin{equation}
\begin{split}
e[n] &=~ x[n] - \widehat{x}[n]\\
 &=~ x[n] - a_1x[n-1] - a_2x[n-2].
\end{split}
\end{equation}
The coefficients $a_1$ and $a_2$ must be chosen to minimize mean square error, which is $\mathsf{E}[e[n]^2]$.

To minimize mean square error, the error must be stochastically orthogonal to the random variables that make up the estimate:

\begin{eqnarray}
\mathsf{E}[e[n]x[n-1]] &=& 0,\\
\mathsf{E}[e[n]x[n-2]] &=& 0.
\end{eqnarray}
These will give us linear equations whose solution reveals the proper choice of $a_1$ and $a_2$.

\begin{equation}
\begin{split}
\mathsf{E}[e[n]x[n-1]] &=~ \mathsf{E}\left[\left(x[n] - a_1x[n-1] - a_2x[n-2]\right)x[n-1]\right)]\\
&=~ \mathsf{E}[x[n]x[n-1]] - a_1\mathsf{E}[x[n-1]^2] - a_2\mathsf{E}[x[n-2]x[n-1]]\\
&=~ R[1] - a_1R[0] - a_2R[1]
\end{split}
\end{equation}
Since we want this quantity to be equal to zero, we have the linear equation
\begin{equation}
R[0]a_1 + R[1]a_2 = R[1].
\end{equation}

\begin{equation}
\begin{split}
\mathsf{E}[e[n]x[n-2]] &=~ \mathsf{E}\left[\left(x[n] - a_1x[n-1] - a_2x[n-2]\right)x[n-2]\right)]\\
&=~ \mathsf{E}[x[n]x[n-2]] - a_1\mathsf{E}[x[n-1]x[n-2]] - a_2\mathsf{E}[x[n-2]^2]\\
&=~ R[2] - a_1R[1] - a_2R[0]
\end{split}
\end{equation}
Since we also want this quantity to be zero, we have the linear equation
\begin{equation}
R[1]a_1 + R[0]a_2 = R[2].
\end{equation}

We now have the system of equations
\begin{equation}
\left(\begin{array}{cc}
R[0] & R[1]\\
R[1] & R[0]
\end{array}\right)
\left(\begin{array}{c}a_1\\a_2\end{array}\right)
=
\left(\begin{array}{c}R[1]\\R[2]\end{array}\right),
\end{equation}
whose solution is
\begin{equation}
\left(\begin{array}{c}a_1\\a_2\end{array}\right)
= \frac{1}{R[0]^2 - R[1]^2}
\left(\begin{array}{c}
R[0] & -R[1]\\-R[1] & R[0]
\end{array}\right)
\left(\begin{array}{c}R[1]\\R[2]\end{array}\right)
\end{equation}

I leave it as an exercise to compute the correlations $R[0]$, $R[1]$, and $R[2]$.

An aside: In higher-order LPC we take advantange of the very nice structure of the matrix in the linear equations (it's a Toeplitz matrix). That structure allows one to use the so-called Levinson-Durbin algorithm to solve it numerically. While the matrix is real symmetric and thus is numerically amenable to numerous algorithms, Levinson-Durbin is the one used in practice because of its speed and because it yields other useful quantities as it runs. These reflection coefficients (RCs) that it yields along the way, are proxies for the linear prediction coefficients (LPCs).

I should note here that there are sometimes differences in sign. See, for example, the help page for MATLAB's lpc function. But once the signs are chosen, there is consistency afterward. Check the definition of LPC in your example before starting.

The filter to compute the error is $x[n] - a_1x[n-1] - a_2x[n-2]$, so the impulse response is $\mathbf{h} = (h_0, h_1, h_2) = (1, -a_1, -a_2)$. The corresponding frequency response is
\begin{equation}
H(e^{i\omega}) = 1 - a_1e^{-i\omega} - a_2e^{-2i\omega}.
\end{equation}

In some applications (such as vocoders), the error ($x[n] - \widehat{x}[n]$) is assumed to be a noise process. If the LPCs have been received, then a realization of the noise process is generated, and an approximation of the original signal is created by using the inverse filter, which is an IIR filter with $z$-transform
\begin{equation}
\frac{1}{H(z)} = \frac{1}{1 - a_1z^{-1} - \cdots - a_kz^{-k}}.
\end{equation}


GCC-PHAT (Generalized cross correlation MATLAB)

I remember I was facing a similar problem some time back. I am pretty sure it has to do with the DC and nyquist values when performing the equalisation and then the inverse. The first element in the array which is becoming 1.0 is the DC value.

I have this code which I’m pretty sure works for GCC-PHAT:

function [r,tau] = gccphat( x, y, fs )
M = max(numel(x),numel(y));

%%Transform both vectors
% X = fft(x,2^nextpow2(2*M-1));
% Y = fft(y,2^nextpow2(2*M-1));
%
% % Compute cross-correlation
%
% R = X.*conj(Y);
% c = ifft(R./abs(R));

%%
N = 2*M-1;
Nfft = 2^nextpow2(N);

R = bsxfun(@times, ...
        fft(y,Nfft), ...
        conj(fft(x,Nfft)));
rtmp = fftshift( ...
        ifft(exp(1i*angle(R))) ,1);
r = rtmp(Nfft/2+1-(M-1)/2:Nfft/2+1+(M-1)/2,:);

lags = (-(N-1)/2:(N-1)/2).';
lags = lags/fs;
[~,idx] = max(abs(r));
tau = N/(2*fs)+lags(idx);

end


Based on Schuyler Smith's presentation on the paper, I think they calculated signed Gaussian derivative "gradients" at 4 orientations 0°, 45°, 90°, 135°, and then got from that 8 unsigned gradients at orientations 0°, 45°, 90°, 135°, 180°, 225°, 270°, 315°.

If you have gradients $\frac{\partial}{\partial x}$ and $\frac{\partial}{\partial y}$ along coordinates $x$ and $y$ at orientations 0° and 90° (so that they are orthogonal, that is, perpendicular), you can calculate the gradient $\frac{\partial}{\partial x_\phi}$ with respect to a coordinate $x_\phi$ at any given orientation $\phi$ by:

$$\frac{\partial}{\partial x_\phi} = \cos(\phi)\frac{\partial}{\partial x} + \sin(\phi)\frac{\partial}{\partial y}.\tag{1}$$

The same way you can get the output of a Gaussian derivative filter at any orientation by taking a $\cos(\phi)$ and $\sin(\phi)$ weighted sum of the outputs of the $x$ and $y$ coordinate orientation Gaussian derivative filters.


might look at:

B. F. La Scala and R. R. Bitmead, "Design of an extended Kalman filter > frequency tracker," Acoustics, Speech, and Signal
Processing, 1994.

  Abstract: This paper discusses the problem of designing a frequency tracker using the extended Kalman filter (EKF). The design
  trade-off between balancing noise rejection and tracking at a maximal
  slew rate is considered, as are the effect of varying design noise
  covariance values. The performance penalties for over- and
  under-design of noise covariances are examined and theoretically
  supported design guidelines are suggested.

There's also the book :

Quinn, Barry G., and Edward James Hannan. The estimation and tracking of frequency. Vol. 9. Cambridge University Press, 2001.

As Peter noted, in SONAR, tracking low SNR tones is done but they typically run off a DFT.  If you want to go this route, you would have to do some experimentation on bin widths and short term averaging.  Fundamentally, you need to filter the signal down to a band that results in a positive SNR.  The band can contain more than a few bins, but the signal has to display a reliable peak when it occupies the bin.  A much more expensive approach would be a bank of chirps.

I would also suggest, for your remark about process and measure noise being correlated.  One of the reasons I like the book is that they derive the KF directly from the case where measurement and process noise is correlated.

Kailath, Thomas, Ali H. Sayed, and Babak Hassibi. Linear estimation. Vol. 1. Upper Saddle River, NJ: Prentice Hall, 2000.


Please help me understand this paper about Discrete Wavelet Transform!

Basically, the method decomposes a signal into one approximation $A_7$ and seven details $D_7$,..., $D_1$. Since the wavelet is orthogonal, the energy of the signal is equal to the sum of the energy of each subband. To me, equation (1) is not correct, one should not sum approximations over levels. Approximation coefficients $A_7$ are discarded, they kind of remove the low frequency or baseline parts. They keep the four most energetic detail subbands, which they order by decreasing energy. Then they set the others to zero, and reconstruct the signal by inverse wavelet transform. So no level is missing, they are just erased. This turns out to be a very crude form of denoising, as they finally only keep the central frequencies, in the interval $[1/2^7\; 1/2^3]$. This could be a misuse of wavelets, as a simple bandpass filter could have been better.


FIR filter (from zeros to filter coefficients)

First of all a discrete filter having a zero somewhere means it's Z-transform (which is a polynomial) acquires the value 0 there.

If you remember, there is a theorem in calculus saying: for a polynomial to have real coefficients any complex roots which are not on the real line must come in complex conjugate pairs. So if $z_0 = e^{j\omega /6}$ then $z_1 = e^{-j\omega / 6}$ and that gives : $$(z - z_0)(z - z_1)C = (z - e^{j\omega /6})(z - e^{-j\omega /6})C$$

Then to find the constant $C$ solve the equation with the absolute value of the polynomial equals 1 for $z = e^{j\pi}$.


HINTS:

Just multiply the terms to get $H(z)$ in the form $$H(z)=h[0]+h[1]z^{-1}+h[2]z^{-2}$$ If you know the four types then you should immediately see which type it is.
Your frequency response looks good. The term in parentheses is real-valued. Does it change sign? You can directly write down magnitude and phase from that representation (just watch out for a negative sign, and what that means for the phase).


Avoid saturated pixels in cross-correlation computation

I have finally found an aswer in a technical report "Image Alignment and Stitching: A Tutorial" from Richard Szeliski.

A weighted SSD (Sum of Suqared Differences) function is presented on page 13:

$$E_{WSSD}(u)=\sum_{i}w_{0}(x)w_{1}(x_{i}+u)\left[I_{1}(x_{i}+u)-I_{0}(x_{i})\right]^{2}$$

where $w_{0},w_{1}$ are spatially-varying weight functions (for template and the image). $u$ is a template offset.

The result may be divided by the overlap area to avoid bias towards smaller overlap.

This however works for SSD metric, not NCC. But it may suffice once images are normalized beforehand (using contrast/brightness adjustment computed using linear regression over corresponding pixel intensities).


A few observations

Your pictures looks odd. "Wave 1" doesn't look like an original audio file but something that has already been encoded/decoded and aggressively lowpass filtered in the process. Your reconstructed audio (Wave 2) actually as MORE high frequency energy, which I suspect is simply just quantization noise from the other bands. The expected spectrum of music is roughly pink (3dB/octave) up to about 10 kHz and then a steeper decline (maybe 10dB octave) to 20 kHz and above. There shouldn't be a step function in the spectrum.
The standard method of recovering for lost high frequency spectrum is Spectral Replication. That's been already implemented in some of the more modern codec standards (e.g. AAC) although not in MP3. See for example http://www.stockis.se/AES112th_SBR.pdf
If you want to do any commercial work in this area, I highly recommend studying the patent situation carefully.  A massive amount of work has already been done in this field and since the license revenue is quite substantial, the work has been well protected by patents and there is incentive to enforce them.


From BFMatcher constructor documentation:

  NORM_HAMMING should be used with ORB, BRISK and BRIEF, NORM_HAMMING2 should be used with ORB when WTA_K==3 or 4 (see ORB::ORB constructor description)

And ORB constructor documentation:

  WTA_K - [...] Other possible values are 3 and 4. [...] Such output will occupy 2 bits, and therefore it will need a special variant of Hamming distance, denoted as NORM_HAMMING2 (2 bits per bin).

In other words the right choice is the one that matches your ORB configuration (NORM_HAMMING when WTA_K==2, and NORM_HAMMING2 when WTA_K==3 or WTA_K==4), otherwise you'll wind up with incorrect norms and correspondingly poorer/inconsistent matching performance.


What is the slope of a transfer function?

Okay. I think i got it after some clarification from Dan. Thanks Dan. !

I think what i was confused with was the step size in my frequency vector.. i had thought to have the final result in db/mhz.. the steps should be in 1 Mhz.. but i now realize that the freq vector could be in whatever steps i like but if its in Hz then i need to have my final answer normalized by 1e6 to get from  db/Hz --> db/Mhz

So this is the code i finally used..

fs_ADC=200e6;
N=2^11;
step_size=fs_ADC / N ;
f = 0: step_size: fs_ADC;
H = fir1(f);
H_db=20*log10(abs(H));
slope=1e6*diff(H_db)./diff(f);  % db/Hz--> db/mhz
plot(f(1:end-1)/1e6,slope)

Figure result till 30 MHz is shown below:


How does a matched filter know about input signal beforehand?

A matched filter is used when the shape of the "desired" part of the received signal is known. Under the additive white Gaussian noise (AWGN) channel, a matched filter is the optimal detection scheme in the maximum-likelihood sense. This is a commonly-modeled case, so matched filtering is a typical signal-processing operation.

Based on your question, I think you're misunderstanding what a matched filter is: it doesn't require full knowledge of the exact transmitted signal. Instead, it helps to sift out the transmitted signal from the noise that is received with it, allowing the receiver to make more accurate decisions on what was actually transmitted.

Consider a baseband communication system that transmits a rectangular pulse with amplitude $1$ to communicate a zero and $-1$ to communicate a one. The appropriate matched filter in this case would be a rectangular pulse that matches the shape that is used by the receiver.
When the train of transmitted rectangular pulses is passed through the receiver's rectangular matched filter, the result is a train of triangular pulses (the autocorrelation of a rectangular pulse is a triangle pulse).
The receiver then makes symbol decisions based on the amplitudes of the triangular peaks that it observes. At this point, the SNR is maximized, as the matched filter tends to pass the signal of interest while averaging out any additive noise.

This is a simple, but representative example. The matched filter doesn't carry any of the transmitted information in itself (as you noted, that would be nonsensical); instead, it is used as an optimal means of extracting the transmitted information from the observed waveform.

I would also call your attention to this very good answer to a previous question about matched filters. It contains many of the diagrams that would make the above discussion easier to understand but I didn't have time to produce. It's a great resource for additional details.


algorithm for feature extraction from audio signals

For speaker recognition, the features you should probably start with are MFCC.

There are libraries offering MFCC extraction modules, such as YAAFE, aubio (C/C++), the MIR toolbox or Dan Ellis' implementation (Matlab) - and of course speech recognition frameworks (Sphinx, HTK) provides MFCC extraction tools.


How to calculate time domain SNR using known sequence

The peak of the cross correlation should be the transmit signal power times the channel's attenuation.
Realizing that, it's really just
\begin{align}
\text{SNR} &= \frac{P_\text{signal}}{P_\text{noise}}\\
&= \frac{P_\text{signal}}{P_\text{received} - P_\text{signal}}\\
&= \frac{P_\text{tx}\cdot a_\text{channel}}{P_\text{received} - P_\text{tx}\cdot a_\text{channel}}\\
&= \frac{P_{crosscorr,max}}{P_\text{received}-P_{crosscorr,max}},
\end{align}
which is
=Rxy(tm)^2 / (Py-Rxy(tm)^2)

in your notation.
So, either that book is wrong or you're missing something there.


Bandwidth is of course occupied by the signal itself. The modulation design of the signal determines the bandwidth. That's why pulse shaping techniques are adopted  by many different kind of modulation techniques.

However, our bandwidth we occupy has a DIRECT relation with our maximum bit rate. In Shannon's theorem:

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

where

C = channel capacity (bits/s),
B = channel bandwidth (Hz),
S = signal strength (watts),
N = noise power (watts).

This theorem tells us that we CANNOT transmit more bits per second through that channel in good condition.

Hence, our signal determines our bandwidth, and our bandwidth determines our maximum achievable bit rate.

It depends on some number of things like Inter Symbol Interference and Nyquist Theorem. It's relatively a long topic to explain each element. Please take a look at the terms I have mentioned (Nyquist and Inter-Symbol Interference) and if there is a spot if you have a hard time to understand, please don't hesitate to ask.

As always, Have a nice day.


Is SSIM a measure of mutual information?

...is the quality measure SSIM between two images (one base-line, one distorted) directly correlated to the mutual information between the two [?]

The short answer to this is "yes" but it tells you nothing. Because the question is a little bit ill posed.

You have two signals, let them be signals or images, it doesn't matter, you have two observations. And you want to assess how much they are alike.

What does alike mean? You might have an idea about this already, if it is sound do they sound similar? If it is images, do they look similar?

So, things look alike when they co-vary...Sound pressure over time covaries between two sound clips, brightness (reflected light flux) over space covaries between two images.

One of the first metrics of assessing how much alike two signals (images) are that someone might encounter is cross correlation .

If you employ (normalised) cross correlation, on your baseline image (i.e calculate its autocorrelation), you will observe the highest value at the centre of the spatial cross correlation. A value of 1.0.

Now, notice here, we decoupled the concept of co-variance from the way by which we assess it. Co-variance is this notion that things might be going up and down at the same time/place. Cross correlation is a quantification of that notion.

And, there are others. Cross correlation, mutual information, transfer entropy, granger causality, partial correlation, ... and many others.

ALL of these are correlated with each other. Their ranges, sign and even interpretation might be different, but if you adjust for all of these, in general, all of these metrics go in one direction when the observations co-vary and in another when they don't co-vary.

Let's turn to the SSIM, the SSIM is the product of three factors: luminance, contrast and structure with the following expressions:

(Image captured on wikipedia)

Where the $\mu$'s denote means and the $\sigma$'s denote variances with the $\sigma_{xy}$ in particular denoting covariance. By the way, you might want to notice the similarity of the $l,c,s$ components with the covariance and even the correlation coefficient.

Luminance doesn't even count, if I remove the mean from two images, the luminance component goes to zero. It is a little bit more difficult to zero the Structure but I can do it, if the two images do not covary with each other, a condition that is met if the two images are pure white noise. This will leave the contrast with some value that would have to be equal to the separate variances.

Long story short: SSIM contains a covariance term, covariance and mutual information are quantifications of the same concept, therefore, SSIM and mutual information are correlated.

BUT, this tells you nothing. The point is how are these two correlated and the interpretation. So, there is a different mapping between cross correlation and mutual information because mutual information has that $\log$ factor in it and, if it is mutual information then would it assume memory of some length?

So now we are on the next level, WHAT are you trying to achieve? To then move on to WHICH aspect of your signal do you need to quantify (and how) in order to get there?

Hope this helps.


Spatial spectrum of EEG data: non-uniform DFT?

QUESTION 1: Do you think an nonuniform DFT makes sense for such an electrode grid?

The answer to that specific question is yes.

  QUESTION 2: The reason I would like to used the nonuniform DFT is that I would like to extract the most parsimonious frequency representation. ...

With this, I suppose that you are after the most parsimonious frequency representation in the spatial domain. In other words, the answer to the question, "How many waves and at what frequencies and amplitudes compose the "surfface" of the electric field potential that is detected by the electrodes at the scalp". If that is accurate (?), it will be a useful reference later on.

  ... For the non-uniform case it is unclear to me, how many non-interpolated coefficients one can extract and what these coefficients are. Any idea?

In the case of a type 2 NUDFT, it is still $\frac{N}{2}$ discrete frequencies.

In the non-iterative cases, there is an implicit interpolation step anyway in order to be able to then use the Fast Fourier Transform, whether this is done via oversampling or "gridding". As is mentioned in "Using NFFT 3 – a software library for various nonequispaced fast Fourier transforms" (page 15):

"The chief idea of the NFFT algorithm is to use standard FFTs in combination with an approximation scheme that is based on a window function $\phi$. This function needs to be mutually well localized in time/spatial and frequency domain."

  QUESTION 3: I was checking the NFFT implementation, however I am confused by the fact that the grid points must be constrained to the [-1/2 1/2] grid. Is it enough if I scale down my grid to this range?

I cannot find that specific reference but this probably refers to normalised frequency. If you notice (1) from here, the only thing that changes between the uniform and non-uniform Discrete Fourier Transform (DFT) is the substitution of the time variable with a sequence. In other words, the discrete sinusoids of the DFT (the $k$s) still "run" at the frequencies they would run anyway but you take extra care to correlate the acquired signal samples at the right phase. That is, take their relative timestamps into account, rather than simply putting them one next to the other as if they were acquired at equi distant sample intervals.

  QUESTION 4: The electrodes are actually placed on the scalp surface. The grid below is a 2-D representation of the 3D scalp surface. Is there a way to compute a nonuniform DFT directly on the 3D grid?

Yes, with the non-uniform spherical version of the FFT. The brain is not a sphere though.

Now, let me suggest a different approach to achieve what you are trying to do.

The reason I asked you "...how do you envisage this "spatial spectrum"..." is because given a $N \times M$ matrix $X$ of acquired data over $N$ samples and $M$ electrodes, there are $N$ different "spatial spectra" that you can obtain. This is because, time and space are two different variables and obtaining the spatial spectrum means that you obtain the two dimensional Fourier Transform over all stations ("theses") at the same time instance. So, you freeze time and you obtain the spatial spectrum of *the electric potential measured across the scalp".

The next natural thought here is to say "that's not a big deal I will integrate over time and possibly over frequency as well and therefore will be extracting the spatial spectrum of the strongest wave activity across the scalp". Which leads us to things like filtering at a brain wave rhythm and then evaluating the distribution of those waves across the scalp and then, by that integration step, we essentially "freeze" activity over a unit of time and that allows us to obtain a spatial spectrum.

Eventually, you would be getting closer and closer to something that a technique like LORETA gives you for free.

LORETA estimates a small number of sources that are responsible for the electric activity across the scalp and projects that activity back to the scalp to form a potential map. So, essentially, if you were to use LORETA to estimate the distribution of the field you also get a nicely interpolated potential from the small number of sources that the low number of electrodes can estimate and you can then do a two dimensional FFT on that "reconstruction".

For more information you might want to have a closer look at this.

Hope this helps.


I would try a median filter.

Let your original signal be $f[n]$.
Median filter $f[n]$ using $N$ pixels, where $N > 2 \times S + 1$, where $S$ is the maximum number of samples in the spike. The resulting signal, lets call it $g[n]$ should have all the spikes removed.
Find the absolute of the difference between the two signals, $h[n] = |f[n] - g[n]|$. This signal represents the spikes.
Count the number of positive transitions in $h[n]$ that are above a threshold. This is the number of spikes.

Example:

Original signal $f[n]$

Median filtered signal $g[n]$

Absolute difference $h[n]$


Your reasoning sounds correct (with unit considerations), you can see a discussion on both continuous and discrete aspects of the Fourier transform in What are the units of my data after an FFT?

What is unclear to me is why  you expect the PSD to be in unit $\frac{m^2}{Hz}$? Powers are powers: if a quantity has unit $\textrm{u}$, its square has unit $\textrm{u}^2$. However, when dealing with physics, the quantities you mention are often somehow periodic, and evaluated via averaging over one period, or asymptotically. So basically you had a normalization factor, proportional to the inverse of a time duration, or $\textrm{s}^{-1}$. Hence, possibly, your $\textrm{Hz}$ factor.


Is there an asymmetric crypto analog to direct sequence spread spectrum?

(This does not answer your question which is a good one, so I may move this to a comment below your question)

Warning! Typical Direct Sequence Spread Spectrum that uses linear spreading codes (based on linear feedback shift registers implemented from irreducible and primitive polynomials) are NOT a good cryptographic technique. See this post where the approach to break such codes is detailed:

Minimum Output Samples needed to crack a "Gold Code" Generator (Dual LFSR)

This is because although the signal will indeed be "buried in the noise" at certain reception points, it is very rare to not be able to get closer to the transmitter (either spatially or with highly directional antennas) such that a positive SNR can be established enough to get the relatively small number of consecutive sequences to crack the code.

Use of non-linear codes together with data layer encryption are recommended. I wanted to mention this due to the common misconception that DSSS alone is a good communication security technique.


Estimation of occupied frequency bins (location of non-zero fourier coefficients)

Is there a method to identify location of non-zero Fourier coefficients of a signal (just locations, not values) with minimum computational cost (less that computations of FFT)?

So, first of all, computation of FFT:

Remember, the $10N$-point FFT, which you need for the resolution you demand:

  Suppose my signal is a vector of length $N$ associated with some white noise.Suppose the demanded frequency accuracy is
  $$\frac{f_\text{sample}}{10N}$$

has complexity of $\mathcal O\left(10N\log(10N)\right)$.

You've got a sparsity of ca. $\frac1{20}$:

  comment clarify I give an example, imagin we have signal with lengh of 1000 or 2000. We take its DFT and there are less than 50 coefficient non-zero.

which means that it's not that sparse, and all algorithms exploiting sparsity by throwing linear algebra at the problem can't really benefit much from the situation; a $K$/$L$-rank matrix/matrix operation typically has best-case complexity $K\cdot L$, so considering $K=N$  this can only work out¹ in the favour of the sparse operation if $L<10\log(10N)$.  Now you gave $N\approx 1000$, and let's chose a very modest logarithm base here – $\log_{10}$ – so, you'd end up with the condition $L < 10\log_{10} 10000 = 40$, and as you can see, with $L\approx 50$, this is not generally the case – this can really only work out if your sparsity-exploiting "tone finder" is very, very, very efficient.

But, you'll say, there's superefficient parametric estimators such as ESPRIT, that can get a much higher resolution with less samples!

Yes, that's what I'd have suggested in your situation, in general, too, but you've got to realize that these subspace/Eigenvalue based approaches work on a autocorrelation signal model – i.e. they need an estimate of the auto-correlation-matrix, which you typically get by calculating the dyadic product of the $N$-sample vector – a operation with complexity $\frac 12 N^2$, i.e. much much more than the complexity of the FFT.

So, to answer your question: No, there's to my knowledge no estimator with lower computational cost, but there's more mathematically efficient estimators than the DFT; however, you might actually get a $\frac1{10}$ resolution with less than $N$ samples, but really, you need to have a rank of your matrix estimate sufficiently higher than the number of tones you're looking for, so I'm afraid this won't change that much.

¹ I'm well aware that's not really how Landau symbols work. However, FFT implementations are really highly optimized, and so is LAPack, and thus, my gut feeling says I can really use the complexity limits as "constants".


Is it valid to increase amplitude (and presumably FFT quality) by simply scaling the data?

This can be a valid approach. You're observing a very practical issue that arises often when using fixed-point (i.e. integer) arithmetic (although it can happen in floating-point also). When the numeric format that you're using to perform calculations doesn't have enough precision to express the full range of values that can arise from your calculations, some form of rounding is required (e.g. truncation, round-to-nearest, and so on). This is often modeled as an additive quantization error to your signal.

However, for some combinations of algorithm and rounding scheme, when the magnitude of the input signal is very low, it's possible to get what you observed: a large number of zero outputs. Basically, somewhere in the sequence of operations, the intermediate results are becoming small enough to not break the threshold required to quantize to a nonzero level. The value is then rounded to zero, which can often propagate forward to the output. The result is, as you noted, an algorithm that generates a lot of output zeros.

So can you get around this by scaling the data up? Sometimes (there are very few techniques that work all the time!). If your input signal is bounded in magnitude to a value below full-scale of the numeric format (16-bit signed integers run from -32768 to +32767), then you could scale the input signal up to more fully use the range available to it. This can help to mitigate the effects of roundoff error, as the magnitude of any roundoff error becomes smaller in comparison to the signal of interest. So, in the case where all of your outputs are getting rounded to zeros internally to the algorithm, this may help.

When can such a technique hurt you? Depending upon the structure of the algorithm's calculations, scaling the input signal up might expose you to numeric overflows. Also, if the signal contains background noise or interference that is larger in magnitude than the algorithm's roundoff error, then the quality of what you get at the output is going to be typically limited by the environment, not error introduced in the computation.


Your algorithm, as described, is not correct.  You are simply counting the number of negative numbers, not the number of zero crossings.  For the particular data set that you use, B1, the two numbers are the same, but in general they won't be.

The way to do it is to multiply each number by the number after it, and if that number is negative, then you have a zero crossing.

zeroCrossings = sum((data(1:end-1) .* data(2:end)) < 0);

You also need to add some logic to handle samples that are exactly zero, but hopefully you get the idea at this point.


Obtaining the magnitude of the frequency response by plugging $e^{jω}$ into the z-domain transform function?

It's  derived using Euler's equation http://en.wikipedia.org/wiki/Euler's_formula. You start with $$H(z) = (\frac{1}{1 - \frac{1}{5}z^{-1}}) (\frac{1}{1 + \frac{1}{2}z^{-1}})$$

Then plug in $e^{-jwT}$ for $z^{-1}$ and you get

$$H(z) = (\frac{1}{1 - \frac{1}{5}e^{-j\frac{\pi }{2}}}) (\frac{1}{1 + \frac{1}{2}e^{-j\frac{\pi }{2}}})$$

Now according to Euler we simply have $e^{-j\frac{\pi }{2}} = -j$ so the whole thing simplifies to

$$H(z) = (\frac{1}{1 - \frac{1}{5}\cdot (-j)}) (\frac{1}{1 + \frac{1}{2}\cdot (-j)})$$.

Utilizing $(-j) \cdot j = 1$ we can multiple each part of the fraction with j/j and get
$$H(z) = (\frac{j}{j - \frac{1}{5}}) (\frac{j}{j + \frac{1}{2}})$$

Then we can take the magnitude as

$$\left | H(z) \right | = \frac{\left | j^{2} \right |}{\left | j-\frac{1}{5} \right |\cdot \left | j+\frac{1}{2} \right |}$$

This isn't exactly the same formula that you have given but since magnitude j and j^2 are the same, the final result comes out to be the same.


You'll not experience aliasing – in the end, the IDFT result per definition only contains frequency components that can be represented by the original frequency domain.

If you will, just stop worrying about aliasing: the DFT and inverse DFT are just base change matrices in $\mathbb C^N$, and are invertible.

Aliasing can only happen when something is non-invertible. This is not the case here; so, no problem.

What you might see is leakage, but you're already taking care of that by windowing.


Take the log of your extremities, find your points as if it were on a linear scale beween these extremities, take the exponential of the result.

Example:

The log of your extremities are $\log 100$ and $\log 4000$. Your 10 points on a linear scale between these values would be: $\log 100 + k \frac{\log 4000 - \log 100}{9}$, $k \in \{0, 1, \ldots, 9\}$.

Taking the exponential: $e^{\log 100 + k \frac{\log 4000 - \log 100}{9}}$. Notice that the first and last values are 100 and 4000 as expected.

This works well with linspace in programming languages which have this primitive, for example with python/numpy:

numpy.exp(numpy.linspace(numpy.log(100), numpy.log(4000), 10))

The second part of your question is not clear. To generate more points by preserving the same spacing, continue increasing $k$ above 9.


Two basic questions related to complex Fourier series

First, or second, with complex notation and the binomial theorem,

$$\sin^5 x=\left(\frac{e^{ix}-e^{-ix}}{2i}\right)^5 = \frac{e^{i5x}-5e^{i3x}+10e^{ix}-10e^{-ix}+5e^{-3ix}-e^{-i5x}}{32i} = \frac{\sin(5x)-5\sin(3x)+10\sin x}{16}$$

On the other hand,

\begin{align}
\sin(5x) &= Im(e^{i5x}) = Im([\cos x+i\sin x]^5) \\
&= 5\cos^4x\sin x - 10\cos^2x\sin^3x + \sin^5x\\
&= 5(1-2\sin^2 x+\sin^4 x)\sin x-10(1-\sin^2x)\sin^3x+ \sin^5x\\
&= 5\sin x - 20\sin^3x +16\sin ^5x
\end{align}

so you got the coefficients right, but the notation reversed.

In the second task, is the star connected to convolution, cepstrum transforms or similar?

No, it is just $g(x)=\exp(α x)=e^{α x}$ defined on the interval $[-\pi, \pi]$. Then

\begin{align}
c_k
&=\frac1{2\pi}\int_{-\pi}^\pi g(x)e^{-ikx}dx
 =\frac1{2\pi}\int_{-\pi}^\pi e^{(α-ik)x}dx
\\
&=\frac1{2\pi}\left[\frac{e^{(α-ik)x}}{α-ik}\right]_{x=-\pi}^{x=\pi}\\
\end{align}


If it were a low-pass filter, I would say you need to turn the resonance (or $Q$) down enough that both poles are real.

What's the order of your bessel filter?

But a bandpass filter is necessarily resonant.  i might think that a BPF must have poles that are complex.  (now that i think about it, a BPF does not necessarily have poles that are complex conjugate.)

if the BPF has gain exceeding 1 (or 0 dB) for its resonant frequency, then then it should be expected that the output will be larger than the input at that resonant frequency.  but if the gain is less than 1 for all frequencies, i dunno why it would overshoot at that frequency.

but there are zeros for BPF also.  one at DC and another at Nyquist.  so maybe you have to turn down the $Q$ to $\frac12$ so that no pole is "ringing".


Variational Regularization Method in Image Processing

This is an example of the Fidelity Term and Prior Term model.
In many Inverse Problems we assume some model on the additive noise. This part is modeled by the Fidelity Term ($ \mathcal{D} \left(A \boldsymbol{f}, \tilde{\boldsymbol{g}} \right) $ in your example). For Gaussian Noise it is given by Least Squares Term:
$$ \frac{1}{2} {\left\| A \boldsymbol{f} - \tilde{\boldsymbol{g}} \right\|}_{2}^{2} $$
The prior terms is usually something else we know about the data in question. It can be something about its distribution, about its derivatives or in the case of Variational Models something about its response to some operator (For example, Derivative).
In 3.12 above they use the Tikhonov Model which I wrote about in your other question.


Firstly, for a response to be classified as the response to an impulse. The input first must qualify as an impulse to the system. For ex: a quick clap could be the impulse input to measure the sound impulse response of the room.

In your case, hitting the hammer at point E must be a "good enough" impulse input for all the materials and lengths.

The issue however is complicated in this case as locations A and B are simply outputs when viewed inline with systems theory for the rope system you constructed. Therefore, the excitation at input should be an "effective impulse" at both locations A and B for you to effectively denote the response of both sensors at A and B to be the response of an impulse.

So the essential question here is to first classify what qualifies as an impulse for rope at E and the suseuqent uniformity of excitation across the length.

If the excitation  at E is not uniform at A and B then you would have to define two effective impulses for your rope system. One for the response at A and other at B.

Once you know the impulse response. Then any input can be written as sum of weighted impulses and if the system is linear and time invariant then the output is simply the sum of delayed responses as below for a discrete time signal

$$ x(n) = \sum_{k=-\infty}^{k=\infty} x(k)\delta(n - k)$$

Let the impulse response to this $\delta(n)$ input be $h(n)$ then the output is simply,

$$y(n) = \sum_{k=-\infty}^{k=\infty} x(k)h(n-k)$$

This is possible both because the system is linear as well as time invariant, so if the impulse response $\delta(n)$ is $ h(n)$ then the response for $\delta(n-k)$ is $h(n-k)$


Why do envelope extraction before pitch detection?

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

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

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


This is (trivially) possible because all Fourier transform like operations are ideally invertible. So what you ask for is actually easily realized by:

$X(e^{j\Omega}) = DTFT\{ IDFT\{X(k)\} \}$

In formulae this is

$X(e^{j\Omega}) = \sum_{l=0}^{N-1} \sum_{k=0}^{N-1} X(k)e^{2 \pi j kl/N} e^{-j\Omega l}$

which can be a bit simplified to

$X(e^{j\Omega}) = \sum_{k=0}^{N-1} X(k) \sum_{l=0}^{N-1} e^{jl (2 \pi k/N-\Omega)}$

or, finally

$X(e^{j\Omega}) = \sum_{k=0}^{N-1} X(k) \frac{e^{2 \pi j k- jN\Omega}-1}{e^{2 \pi j k/N-\Omega}-1}$.

EDIT: The question on whether you need all samples X(k) lies in the behavior of the interpolation function $I(x) = \frac{e^{j N x}-1}{e^{j x}-1}$. This is because the above result is then equal to

$X(e^{j\Omega}) = \sum_{k=0}^{N-1} X(k) \, I(2 \pi k/N - \Omega)$.

I(x) has some "sinc like" behavior, so the question on how many X(k) you need, is really when I(x) has decayed to a value "small enough".


You might be interested in the following paper:

Nehab, D; Maximo, A; Lima, R; Hoppe, H: GPU-Efficient Recursive Filtering and Summed-Area Tables. ACM Trans. Graph. 30(6):176 (December 2011). doi://10.1145/2024156.2024210
http://w3.impa.br/~diego/publications/NehEtAl11.pdf

The basic idea is to recognize that linear filters are associative (at least in theory) and thus you can

divide your input into chunks.
every thread initializes its iir filter state with 0s and processes its chunk.
every thread stores the final state of its iir filter for its chunk.
now (because of associativity) you can figure out what the beginning filter state is of every chunk.
so then you run the filter over the chunk again, this time correctly initialized and get the correct output.

As @jonsca pointed out, you were a bit vague about your use case.  The above all assumes that you can do batch processing on your input.  If you need some kind of real-time filtering none of what I said is relevant.


Bit-reversal equivalence on IFFT (radix-2 Cooley-Tukey)

To answer my own question : no you can't use the exact same algorithm.

In the case here of a decimation-in-time, radix-2 Cooley Tuckey, in-place FFT/IFFT algorithm, we have the choice between two implementations which are exactly equivalent :

The above uses bit-reversed inputs and produce natural order results.

The above uses natural order input and produces bit-reversed outputs.

You can see that the operations are the same : for each butterfly operators, same input data, same twiddle factors. It's just the address of that data (here, the number of lines from the top) that changes.

The same applies for IFFT. When applying the bit-reversal after the IFFT, I would have needed to change the IFFT algo so that it changes the addressing scheme to match the second algo.

To go a bit further on the topic : "ALGEBRAIC FORMULATION OF THE
FAST FOURIER TRANSFORM" from Tran-Thong. (http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=6323745)


Plotting convoluted and original signals on the same graph

Looks like you need to do a couple of things:

Use conv instead of convol
Change the order of the parameters of convol like this: conv( win_response,siggram_n(:, i))
Tell conv that it should return an array the same size as the first parameter conv( win_response,siggram_n(:, i),"same")

If you don't use "same" then you get a convolution result back that is the length of the first array plus twice the length of the second array.  Convol has the parameters in the opposite order, and doesn't have a "same" option.

Also, change the name of your variable conv to something else.  It is bad practice to use keywords as variable names.  In this case, you have to change it else you can't use the conv method properly.

See the documentation for conv


Yes, this is common. Common non-linearities (such as op-amp clipping) preserve the symmetry of the upper and lower section of the waveform; and as such, they do not introduce even harmonics - just odd harmonics.


How good is the Leech lattice for AWGN coding or MSE quantization of Gaussian sources?

Tables 7.1 and 7.2 from Lattice Coding of Signals and Networks by Zamir are essentially this.

Also the figures in this paper: https://arxiv.org/pdf/1103.0171.pdf


So, you've got the audio signal in digital on your PC.

Analyze that digitally; you don't need any measurement equipment to do that.

On the contrary: Measuring properties of the analog audio signal is highly undesirable, because you'd be including the sound card of your laptop in your measurement, but you probably don't care about that even remotely. If you actually wanted to know the properties of your sound card, then you'd need to separately measure that – but you'd do that with a "clean" test signal generated by your PC, not something received over a noise system.


When grabbing the first 10 samples of $h_1[n]$ you are essentially multiplying the whole sequence $h_1[n]$ with the rectangular window

$$
  w[n] =
\begin{cases}
  1 & 0 \leq n < 10 \\
  0 & \mbox{otherwise}
\end{cases}
$$

whose DFT is given by:

$$
\begin{align}
  W\left[k\right]
    &= \sum_{n=0}^{9} e^{-j 2\pi k n / N} \\
    &= \frac{1 - e^{-j 2\pi k \cdot 10/15}}{1 - e^{-j 2\pi k / 15}}
\end{align}
$$

By the circular convolution theorem:

  Multiplication in the discrete-time domain becomes circular convolution in the discrete-frequency domain. Circular convolution in the discrete-time domain becomes multiplication in the discrete-frequency domain.

In other words,

$$
  \mbox{DFT}\left\{h_1[n] \cdot w[n]\right\}
  = \frac{1}{N} H_1\left[k\right] \otimes W\left[k\right]
$$

where $\otimes$ is the circular convolution operator.
This can be seen with the following matlab script (you can play around with h1 to see what happens with other waveforms):

n = [0:14]; % discrete-time index
k = [0:14];
h1 = sin(pi*n/14); % some arbitrary function

% window function
w = [ones(1,10) zeros(1,5)];

% time domain operation
h2  = h1.*w;
% convert to frequency domain
H2t = fft(h2);

% frequency domain operation
W   = fft(w);
H1  = fft(h1);
H2f = cconv(W,H1) / 15;

% Compare H2t & H2f
figure(1);
subplot(2,1,1);
hold off; plot(k, abs(H2t), 'b');
hold on;  plot(k, abs(H2t), 'bx');
hold on;  plot(k, abs(H2f), 'ro');
title('Magnitude');
subplot(2,1,2);
hold off; plot(k, angle(H2t), 'b');
hold on;  plot(k, angle(H2t), 'bx');
hold on;  plot(k, angle(H2f), 'ro');
title('Phase');

Note that if $h_1[n]$ is already no longer than 10 samples (i.e. if those truncated samples are already zero), then windowing in the time-domain will have no effect and similarly the circular convolution with $W[k]$ in the frequency domain will happen to produce $H_1[k] = \frac{1}{N} H_1[k] \otimes W[k]$ (but only as a special case).


Verdict: I conclude scipy's higher leftmost peak is due to both pywt's wavelets' lower amplitude at lower scales and scipy's wavelets' stronger correlation with lower frequencies at lower scales. Overall, they both correlate approximately evenly with irrelevant frequencies (scipy being a little worse), but pywt "lucks out" with its bad correlations being offset with smaller wavelet amplitudes at lower scales, thus the smaller peak.
I'm not exactly sure why scipy's direct recomputation is worse than pywt's resampling; I'd say resolution problem with very limited number of samples early on, and resampling integrated working out better (it's also slightly longer), but haven't investigated explicitly.
For additional assessment on low scales pywt vs scipy, see bottom section. The rest of the answer serves to justify the verdict, but doesn't add much further; the most informative is Fig. 2. Some code is provided, which can be reused to test other transforms.

A more naive kind of norm points to the answer; zero-mean + scale coeffs between $[-1,1]$ per scale:

For the top plot, a $[0,1]$-rescale is equivalent, so it shows the normalized correlation strength with the wavelet irrespective of energy (actually not quite; caveat1). This suggests the difference to be rooted in wavelet amplitudes; let's compare. (-- coeffs code3)

Do they correlate differently? -- For all plots, blue = scipy, orange = pywt.
Scale-standardized: divide each scale vector by own maximum (highest absolute correlation among all timeshifts). Plot absolute value of each scale vector, so timeshifts along horizontal. A greater value at a given timeshift indicates greater correlation of the plotted scale with relative to all other timeshifts at that scale.

Above shows scipy's early scales having comparable correlations at each scale for both halves of the signal, whereas for pywt the right half is dominant.
Does this imply scipy correlates more strongly with lower frequencies at lower scales than pywt (i.e. for same wavelet amplitude)? Yes; we'd get identical scale-standardization if we forced wavelet amplitudes the same. See zoomed amplitude-equalized plot in Fig 2; note how at indices 2 and 4 pywt exceeds scipy, and observe the above gif.

Timeshift-standardized: divide each timeshift vector by own maximum (highest absolute correlation among all scales). Plot absolute value of each scale vector, so
timeshifts along horizontal. A greater value at a given timeshift indicates greater correlation of the plotted scale relative to all other scales at the timeshift.

Above shows, for left half, about same correlation strength for all 20 scales relative to strongest-correlating scale. Right-half fluctuates much more for both, but scipy keeps itself multiples above pywt. This implies most of the energy discrepancy at these scales is due to the right-half.
But why due to right-half? And why is left-half relatively fixed? Figure 1; toward higher scales (lower frequencies, left-half), scipy and pywt agree in both wavelet amplitude and shape.

Wavelet amplitudes comparison
Instead of looking at max amplitude, I define a measure of "mean amplitude": mean of absolute value of tail-trimmed wavelet, where "tail" = any absval 1e7 times less than peak amplitude (instead of strictly zero which is rarer). This is to unbias the mean for wavelets with long tails: (-- code2)
left  = [i for i, v in enumerate(x)       if abs(v) > x.max() / 1e7][0]
right = [i for i, v in enumerate(x[::-1]) if abs(v) > x.max() / 1e7][0]
meanamp = np.mean(np.abs(x[left:-(right + 1)]))

At scales[0], the ratio is 2.58 - when squaring for the energies, that's a x6.64 larger coefficient for scipy due to wavelet amplitudes difference alone. However, the ratio of coefficients' energies at that scale is 8.33, so there's another factor - and this factor is, surely, scipy's stronger correlation with low frequency part of the signal, as shown earlier.
To see how much exactly the amplitude accounts for, we can normalize by it; scale coefficients by ratio of (mean) amplitudes:

We see for most lower scales scipy's energies greater even with amplitude accounted for, leaving only correlation strength as the source of difference. Also interesting to see pywt correlate much more strongly at the first true mode (~scale 50).

Caveat1: Fig. 0 is misleading
Notice scipy jump above pywt at the right edge, yet, if we look at logscale mean amplitudes, pywt's is both greater, and correlates stronger as per the amplitude-normed plot. What's the deal? Recall how scale-standardized plot was defined; this is same, but now for coefficients - scale vectors are rescaled relative to themselves, rather than between pywt and scipy.
A greater value in Fig 0 results if most points aren't far below the maximum (as all points are divided by max) across timeshifts, or that the wavelet at the scale correlates more consistently for each timeshift. So it's more of a consistency measure; here's the last scale, shown for every timeshift:

So... pywt's better at lower scales? -- Uncertain, but pywt wins in two important test cases:
Higher max frequency

Noisy, combined frequency signal


Playing back an audio file at an arbitrary rate with upsampling and downsampling

Some solutions that are ordered by complexity, simplest first.

The simplest solution is the one you have already implemented- pick the nearest value.
The next simplest is linear interpolation.  You use the two nearest values and use a weighted average based on how far the point is from each.  For instance, if the point you want is .3 samples past $x_1$ and .7 samples before $x_2$, the value would be $y = x_1*.7 + x_2*.3$.  This is still pretty simple and should have much better sound performance than your current solution.
The next simplest is quadratic interpolation, or, more generally, polynomial interpolation.  Polynomial interpolation is also sometimes referred to as spline interpolation.
The most complicated solution is Lanczos resampling.  This is where you use a windowed sinc filter to low-pass filter the data.  You pick the particular set of taps to use from the filter based on the particular delay (i.e. at what point you are between the samples).


Why is my DFT/FFT always 90° out of phase?

The Fourier Transform (and the DFT is the discrete form of this with a finite length sequence) is simply a correlation of the input signal with each frequency bin, of the continuous form:

$$X(\omega) = \int_{t=0}^Tx(t)e^{-j\omega t}dt$$

(If we consider a finite time duration, however I will stay with the continuous time version of the Fourier Transform as it will simplify the explanation).

Here we see that for any freqency $\omega$ we are doing a correlation of $x(t)$ with $e^{j\omega t}$  (correlation is a complex conjugate multiplication followed by integration, or summation when discrete).

As a correlation we see that the result will be maximum when x(t) is fully correlated to $e^{j\omega t}$ over the given time interval, and similarly 0 when it is fully uncorrelated (orthogonal).

Now make note of Euler's identity given as:

$$e^{-j \omega t} = cos(\omega t) - j sin(\omega t)$$

And replacing that in our first formula we can expand the Fourier Transform given into two correlations:

$$X(\omega) = \int_{t=0}^T x(t)cos(\omega t)dt - j\int_{t=0}^T x(t)sin(\omega t)dt$$

From this we see that when $x(t) =sin(\omega t)$, (and duration T covers one full cycle or an integer number of cycles), the second term will be fully correlated, while the first term will be zero, and vice versa when $x(t) =cos(\omega t)$.  Any correlation to the second term results in an imaginary result IF x(t) is real.  Any phase shift on x(t), as in $x(t) = sin(\omega t + \phi)$ where $\phi$ is anything except $\pm N\pi/2$ for N any integer, will result in a complex result with $X(\omega)$ having real and imaginary components.

Note the discrete form of this, the DFT, is as given below where the same result can be explained, but may not be as easy to see for some that are less familiar with it.

$$X[k] = \Sigma_{n=0}^{N-1}x[n]e^{-jn\omega_o \frac{k}{N}}$$


Step-by-Step calculation of stationary wavelet transform gives shifted results?

Why normalized filters

Last question first, "Why $1/\sqrt2$": Because it makes the (Euclidean) norm of the filter $1$, so that the whole wavelet filter bank operation (if done right, that is, the decimated version) is orthogonal/isometric. It is a design choice, there is nothing wrong in staying with the integer values and correct the combined factor during reconstruction.

The book example is not really the standard wavelet transform

However, the book you cite is doing the wavelet transform undecimated. Which is a valid mathematical transform, but leaving out the decimation step, the down- and up-sampling, will lead to data bloat. By applying the two filters, the amount of data representing the signal is doubled, so some of that is redundant. This redundancy can be reduced by transmitting (or keeping for further processing) only every second element of these sequences.

Undecimated or otherwise redundant transforms give freedom in how to represent the signal, that is, in some non-linear post-processing one can set some of the data to zero and compensate (more-or-less) lossless-ly by changing other coefficients. However here I think this example is just a teaching aid.

Decimated version of the examples

So taking of cD1 and cA1 only every second entry gives cD1(2:2:end)=[0 0 0 0] and cA1(2:2:end)=[160 160 0 0], together they encode the same signal, which now has more zeros and is thus better suited for compression. In the reconstruction, zeros will be added to fill the gaps, and then the convolution with the reversed/biorthogonal filters is applied. So

[1 -1]*[0 0 0 0 0 0 0 0] = [0 0 0 0 0 0 0 0 0]

and

[1 1]*[160 0 160 0 0 0 0 0] = [160 160 160 160 0 0 0 0 0]

so that in the sum the result is twice the original signal, $(\sqrt2)^2$, with a zero added at the end.

If one shifts the whole process by one place, the reduced signals are [ -80 0 80 0] and [80 160 80 0] so that in the reconstruction one has to add

[1 -1]*[-80 0 0 0 80 0 0 0] = [-80 80 0 0 80 -80 0 0 0]

and

[1 1]*[80 0 160 0 80 0 0 0] = [80 80 160 160 80 80 0 0 0]

to give [0 160 160 160 160 0 0 0 0]. Here we get, additional to the factor 2, also the shift, or an additional zero at the start of the sequence. However, by using the first, third etc. elements in down-sampling, we used implicitly an extension of the sequence S to the left by zero. This additional zero element is what is obtained now explicitly in the reconstruction.

In general wavelet transforms, the wavelet filters have lengths larger than 2, leading to more "interesting" boundary effects that include a shift by about the length of the filter.


There a few variations on the definitions of "energy signal" and "power signal". For example, Sklar in his textook "Digital Communications" defines "energy signal" as a signal with energy $E$ such that $0 < E < \infty$, and a "power signal" as a signal with power $P$ such that $0 < P <\infty$.

However, others (see for example this answer) define "energy signal" as a signal with $E < \infty$.

In your case, $x[n]=0$ for all $n$; then this signal has zero energy and zero power. Whether it is classified as an energy or a power signal depends on which specific definition you adopt. For instance, if you go with Sklar's, then this signal is neither an energy signal nor a power signal.


Audio localization accuracy with two microphones

Each ear have a Frequency Response $H_L(f)$ and $H_R(f)$ and transmission delays $\tau_L$ and $\tau_R$ according the 3D angle from where the sound source is located relative to each ear.

Hence, for a fixed magnitude and frequency sound, $s(t)$, you have the response $H_L(f)s(t-\tau_L)$ and $H_R(f)s(t-\tau_R)$.

Actually, you cannot solve this for every position, or equivalently, you need a head movement (or ear movement in animals) for retrieving the proper angle. Besides, there exist some cone of confusion which means there is actually not a well defined 3D angle, mainly you only can locate 2D cones instead of location.

I will not solve it for you, but from that information you can trivially determine a location considering head movement onto some estimated direction, or an enough rich movement trajectory.


For question 1: apply the definition of time invariant: find the output as normal; find the output with the same input but delayed by $T$
$$
y_1(t) = \frac{dx(t)}{dt}\\
y_2(t) = \frac{dx(t-T)}{dt}\\
$$
Does $y_2(t) = y_1(t-T)$ ?

For question 2: Find a definition of stability and apply it. For example, for a system to be BIBO stable it needs to have

$$
\int_{-\infty}^{+\infty} \left|h(t)\right| dt < \infty
$$

If $h(t)$ is periodic and non-trivial (zero), can that be true?


Frequency of a single sample in a digital signal and aliasing

A single sample being non-zero is known as an impulse.  The effect in the discrete frequency domain depends upon where it is within a the frame.  The DFT will give you the coefficients required to reproduce whatever finite signal is withing the sample frame.  It is misleading, even erroneous, to say the signal is composed of those frequencies.  This can only be said about periodic signal for which the sampling frame is a whole number of periods.  Recall that the DFT and inverse DFT are basically the same, except for a sign change in the exponent.  When you have a pure tone with a whole integer frequency in the sample frame, the result is a single spike in the DFT at the bin corresponding to that frequency.  So when you have a single spike, or impulse, in the signal domain, the result will be a DFT with the bin values being the same as a pure complex tone.

"Alias" is a different concept entirely.  It has to do with the sampling of tones higher than the Nyquist frequency appearing to be tones of a lower frequency.  The easiest way to envision what alias frequency will look like is to imagine the bins of the DFT arranged along the unit circle.  The Nyquist bin will occur at $\pi$ radians.  Any alias frequency has to appear somewhere on the circle.  It may be beyond $\pi$, or even beyond multiple times around the circle.  The frequency it will appear to be (for real valued signals) will be the frequency associated with the bin number of the first half of the circle that has the same projection onto the horizontal axis.

I would recommend that you read my blog articles (link on my profile page) from the start to get a perspective on this material.

Hope this helps.

Ced


This is not a proof, but maybe an intuition why this conjecture can be true for the points in general position. From the properties of rank we know that:
$$
\mathrm{rank}(F) = \mathrm{rank}(F^\top) = 2.
$$
Hence:
$$
\begin{align}
\mathrm{rank}(F+ F^\top) &\leq \mathrm{rank}(F)+\mathrm{rank}(F^\top)\leq 4.
\end{align}
$$
Because $(F+F^\top)\in \mathbb{R}^{3\times 3}$, $\mathrm{rank}(F+ F^\top)\leq 3$. To have any rank deficiency here (e.g. $\mathrm{rank}(F+ F^\top)< 3$), we should have an additional rank constraint on $F$, which we donot seem to have. The only additional constraint is that $F$ has 7 degrees of freedom and that $\det(F)=0$. These non-linear constraints are not sufficient to reduce the rank.


You can simulate an ISI (inter symbol interference) effect, by passing the message signal pulses through a dispersive communication channel model filter.

The dispersive channel will spread the pulses and make them overlap, hence creating ISI.

You can simulate a dispersive channel via nonlinear phase filtering or multi path effects.


Although it looks like it's not being actively maintained, you might check out "aubio" to see if it will meet your needs : http://aubio.org/

  aubio is a tool designed for the extraction of annotations from audio signals. Its features include segmenting a sound file before each of its attacks, performing pitch detection, tapping the beat and producing midi streams from live audio.


How to correct phase of Diagonal Volterra kernels obtained with exponential sweep sine?

A potential pitfall in non-linear modeling is using different signal levels. The sweep test needs to be carried out at a sufficiently high signal level to cover the maximum possible real-world input levels, but not much higher as this will require more kernels for the same effective quality. Usually real-world signals to the system are much of the time at lower level than their peak level, with the 1st order linear kernel dominating.

To me your kernels look just fine. I don't know how common it is to see that a higher-order kernel (3rd order in your plots) is more prominent than the 1st order kernel. The saw-like phase comes from the large initial delay, which gives a steep linear dependency of phase on frequency. Wrapped phase of $\pi$ is the same as a wrapped phase of $-\pi$, and you are seeing wrapped phase, which gives the sawtooth pattern. If you remove the initial delay in the inspection of the results, you should be able to unwrap the phase more easily for inspection.

Most systems you would model by the technique will have right-sided kernels. The systems are causal and cannot store energy for very long after they have received an impulse as input. That would require a delay line. Instead, much of the energy passes through almost instantly as the initial spike, and the residual energy is released as the tail.

I suggest you test with a system you have complete control over, for example a memoryless non-linearity that is the sum of the $N$ first Chebyshev polynomials of the first kind, one for each of the $N$ harmonics you want to generate, including the 0th and 1st harmonics in the count. Then any phase shifts due to errors in the algorithms/software should be clearly visible. If the algorithms work correctly, each harmonic impulse response should be a delayed sinc, and the the coefficients of the memoryless nonlinearity polynomial should be resolved as the amplitudes of the scaled impulse-like kernels.

What I have done in my own implementation is to generate the sweep and its fade in and fade out in the frequency domain by programming the timing into the group delay, and amplitude into the magnitudes, rather than using a direct time domain equation for the sweep. Either way should work and will have its own problems. The frequency-domain approach gives amplitude oscillations near the fades. In the frequency domain it is easy to correct the non-whiteness of the sweep in the inverse filter. (I can't publish my implementation in detail at the moment.)

If there are phase problems, I'd try $\cos$ instead of $\sin$ in the time-domain sweep formula, and fade the sweep in and out. The fades will take care of eliminating transients at the ends. The initial and final phases of the sweeped sinusoid could be made such that the value of the sweep is zero both in the beginning and in the end, but with the fades in and out that is no longer needed for transient reduction, because the fades gently bring the value to zero in any case.

I don't see a problem with the sweep going all the way to the Nyquist frequency, if the non-linearity is analog and the analog-to-digital converter (ADC) has an anti-alias low-pass filter at the Nyquist frequency. This way the harmonic impulse responses and consequently the Hammerstein kernels have only non-aliased frequencies up to Nyquist frequency. It is sufficient information to model the non-linear system in the sub- Nyquist frequency band if the system is compatible with the Hammerstein model. Only the implementation has to be careful not to introduce aliasing. In any case, you'd typically want the linear kernel to have frequency response information up to or at least near the Nyquist frequency, which dictates the sweep end frequency.


Is there a term called "Fourier Beamforming"?

I'm guessing it's used in the following context.

Assuming you have a physical source with some of an amplitude distribution in a plane: you can calculate the polar pattern of that source simply as the Fourier Transform of the amplitude distribution as a function of space (not time!).

Fourier Beamforming would be the inverse process: you start with the desired polar pattern and do an in inverse Fourier Transform: that gives you the desired in amplitudes as a function of space, i.e. it tells where to put sources or sensors and what weights to apply.


Maybe you find this useful:

N=500;                          %signal length
t=0:N-1;
s=20*randn(size(t));
x=s;
p=200*[ones(1,5) -ones(1,5)];   %square artifact template
x(:,101:110)=x(:,101:110)+p;    %add artifact to signal at two random points
x(:,401:410)=x(:,401:410)+p;
figure;plot(x)
threshold = 1e5;                %peak detection threshold
c=xcorr(p,x);                   %cross-correlation
[pk l]=findpeaks(c,'MinPeakHeight',threshold);
l=N-l+1;                        %find the actual lag
y=x;
for j=1:numel(l)                %remove artifact
y(:,l(j):(l(j)+length(p)-1))=y(:,l(j):(l(j)+length(p)-1))-p;
end
hold on;
plot(y,'r--')
legend('corrupted signal','recovered signal')


How to understand Kalman gain intuitively?

I found a good way of thinking intuitively of Kalman Gain $K$. If you write $K$ this way

$\displaystyle \quad\ \bf{K_k} = \bf{P_k^-\, H_k^{\rm T} (H_k P_k^-\, H_k^{\rm T} + R_k)^{-1}}
  = \bf{\frac {P_k^-\, H_k^{\rm T}}{H_k P_k^-\, H_k^{\rm T} + R_k}}$

you will realize that the relative magnitudes of matrices ($R_k$) and ($P_k$) control a relation between the filter's use of predicted state estimate ($x_{k}⁻$) and measurement ($ỹ_k$).

$\displaystyle \quad\
  \lim\limits_{\bf{R_k \to 0}} \bf{{P_k^-\, H_k^{\rm T}} \over\
  {H_k P_k^-\, H_k^{\rm T} + R_k}}\
    = \bf{H_k^{-1}}$

$\displaystyle \quad\
  \lim\limits_{\bf{P_k \to 0}} \bf{{P_k^-\, H_k^{\rm T}} \over\
  {H_k P_k^-\, H_k^{\rm T} + R_k}}\
   = \bf 0$

Substituting the first limit into the measurement update equation

$\displaystyle \quad\
  \bf{\hat x_k} = \bf{x_k^-} + \bf{K_k}(\bf{\tilde y_k}-\bf{H_k}\bf{x_k^-})$

suggests that when the magnitude of $R$ is small, meaning that the measurements are accurate, the state estimate depends mostly on the measurements.

When the state is known accurately, then $H P^⁻ H^T$ is small compared to $R$, and the filter mostly ignores the measurements relying instead on the prediction derived from the previous state ($x_k⁻$).


How can I comprehend the complex coefficient of complex Fourier series?

if you think of $e^{j n \omega_0 t}$ as a "vector" (complex numbers are not quite the same as 2-dimensional vectors, but they look the same, and are both describable with two real numbers) with length 1 and starting out with an angle of 0 (pointing straight to the right, along the positive real axis, and this vector is rotating counter-clockwise at a rate of $n \omega_0$ radians per unit time.  this vector is what you see looking down the spiral in the graphic referred to in the question, looking along the time axis.

then multiplying by $C_n$ changes the length of the vector to $|C_n|$ and the initial angle of the vector is the angle of $C_n$ or $\arg\{C_n\}$.

that's how i would visualize it.


The signal in question is a rectified sine-wave. This signal has a non-zero average, i.e. a DC component. It also has AC components at f, 3f, 5f, etc. where f is the frequency of the non-rectified signal.