The first thing to remember is that the FFT sink does measure actual PSDs, it performs an FFT, which may be scaled however, and displays its mag-square. So when you say "2 Volts p2p", this is incorrect, as it implies certain things, such as Ohm's law, which don't work that way in DSP.

Now, for your result. First, the window you've chosen (probably Blackman-Harris) will attenuate the main lobe height by a couple of dB. If you had chosen a boxcar window (rectangular), you would have noticed the peak would be at -6 dB. Half the power! And now, things fall into place: In a real-valued FFT (as you are running), the power is distributed between the two peaks that are generated at plus and minus the sinusoid's frequency. And, sure enough, if you run the whole thing in complex mode using a rectangular window, you end up at 0 dB.
If you had chosen the QT Frequency Sink, that defaults to showing the entire Nyquist zone and you would have seen the 2nd sinusoid.

Now, is GNU Radio broken for chosen such normalizations? Not really, since they are always arbitrary. But it is always important to realize that the frequency sinks are not true power spectrum analyzers.


The idea is the same as the 1D case. Start with the regular old continuous time fourier transform definition (the integral) and calculate the Fourier transform of $e^{j \omega_0 t} f(t)$. Similarly, calculate the inverse Fourier transform (using the integral) for $e^{j \alpha \omega} F(\omega)$.

Alternatively, do this with the discrete time fourier transform (so the transform is a sum and the inverse transform is an integral).


Differences between the mean filter, trend(mean) removal, and baseline correction filter?

A mean filter is one precise type of linear filters that replaces a value by a weighted combination averaging the neighboring values. A constant offset value, or an affine trend, is "quite" invariant by an averaging filter. Therefore, it may estimates a constant or linear baseline well.
Trend removal might be less well-defined. For the mean, there are many alternatives, like the exponential weighted moving average (EWMA), which is an IIR filter.
Baseline is even less  well-defined. It can dwell between simple instrumental offets and complicated "references" that may include noie.
The last one is the fuzziest to me. But there is a notion of oscillation, that could be incorporated into a filtering model with knowledge of the signal. Tracking and informed filtering (eg Kalman) can be useful for that purpose


Band Limited Table Look Up Oscillator C++

Please note that you cannot implement a band-limited oscillator with a pure (stateless) function like:

Sample s = saw(x);

The reason is that the function must be aware of the rate at which x change to generate a signal with the correct bandwidth

For example, let us assume your sample rate is 48kHz. If you call your function in this order:

saw(0.00000000);
saw(0.00000001);
saw(0.00000002);

You can see that the argument is slowly changing to the point that aliasing is not an issue, so your function could return something like x / M_PI - 1 (assuming input is in [0, 2 pi[ ; and output in [-1, 1]) without any lookup table. But now, if you call it this way:

saw(0.0);
saw(3.0);
saw(6.0);

We are getting close to Nyquist frequency and the output must then be heavily band-limited - actually closer to sin.

So a band-limited sawtooth or square cannot be generated by mere waveshaping of the naive saw. You need your function to be aware of the target frequency so that it can make the right decision.

Another related misconception that can be found in your code is that you assume that you would need only one table. There is just not one single "band-limited sawtooth" waveform that would work for all frequencies. For low frequencies, there is very little different between the band-limited saw and the naive saw. For high frequencies, the band-limited saw is effectively a sine. So you have to define how much memory you want to dedicate to your task, and generate many band-limited tables corresponding to various frequencies. At run-time, you decide which one to use depending on the frequency of the signal. If this is for a musical application in which you don't want to hear the transition from one table to the next as notes rise in pitch, you'll need to add one further layer of interpolation on top of that to crossfade from one table to the next. A typical implementation would use one table by octave or per 8 semitones.


There is no in-build MATLAB function for DHT, but you can implement it very easily using the fft function:

function H=dht(x)
    X = fft(x);
    H = real(X) - imag(X);
end


Minimal number of additions in convolution

One method which attracted my attention recently is Discrete Time Random Sampling.

Which is an approximation method useful to reduce both multiplications and additions in filtering.

Filter and $N$-length signal by a $L$-length filter can be performed in $\mathcal{O}(LM)$  instead of  $\mathcal{O}(LN)$. Where $M\ll N$.
Which is done by picking samples at random in a way described in the paper.

See also Frequency-Shaped Randomized Sampling for a more sophisticated approach which allows you to shape the spectrum of the error.


I was going to suggest the machine learning route but a quick search turned up OSS so you might not have to re-invent the wheel: Identify speakers with sndpeek

Do the speakers talk over each other (which complicates things) or is it mostly one after another?


Jumping between bins of STFTs will lead to the discontinuities seen in your phase estimates.  Interpolating both frequency and phase from a sequence of STFTs is possible, but requires non-standard techniques, such as fftshifts and 2D interpolators.

If you know the frequency range of your desired slow oscillation, then a bandpass filter in conjunction with a PLL follower might provide a better looking phase result using more standard DSP techniques.


Extracting image metadata programatically

Open a few files in a hex editor and check, if the header containing the information is always the same size, i.e. if the image data always begins at the same adress. If this is so, you can just retrieve the header with python ("normal" open, not binary). Then, you have the complete header as a string and can create any data structure you want, using any xml package you want.


Can the Cramér-Rao Lower Bound be constant?

Prologue

The CRLB of an unbiased estimator $\hat \theta$ for the paramater $\theta$ is

$$
\text{var}(\hat \theta) \ge - \frac1 {\text E_x\left\{ \frac{\partial^2 (\ln f(x;\theta))}{\partial \theta^2} \right\}}
$$

with $f(x;\theta)$ being the likelihood function of $\theta$ under the observation $x$. It is the inverse of the so-called Fisher information.

Answer

  Is it possible that the CRLB remains constant for all SNR?

Yes. For that, the expectation (over all $x$, hence the index)

$$\text E_x\left\{ \frac{\partial^2 (\ln f(x;\theta))}{\partial \theta^2} \right\} $$

just needs to be constant.

The "easiest" explanation would be if noise just didn't have anything to do with the likelihood – ie. noise is just not actually observed by the estimator.

But that's kind of least exciting result – "well, if I add a lot of noise to something and then observe something completely unrelated, the noise has no effect" isn't that surprising.

A lot more complicated, is the option that the Fisher information is actually constant, although $f(x;\theta)$ actually depends on the noise.

The expectation operator $E_x$ is definable as

$$
\text E_x \{g(x)\} = \frac1{\mu (\xi)} \int\limits_{x\in\xi} g(x)\,d\mu\text,
$$

i.e. as the integral over all possible $x$ divided by the measure $\mu$ of the set $\xi$ from which $x$ can be chosen.

Assuming smoothness of $x$ and compactness of $\xi\subset \mathbb R^1$ (just to make all the derivates easier to handle and to avoid cascading integrals), this specializes to:

$$\begin{align}
\text E_x \{g(x)\} &= \frac1{x_{max}-x_{min}} \int\limits_{x = x_{min}}^{x_{max}} g(x)\,dx\\
\text E_x\left\{ \frac{\partial^2}{\partial \theta^2}{ (\ln f(x;\theta))} \right\}  \cdot\left({x_{max}-x_{min}}\right) &= \int\limits_{x = x_{min}}^{x_{max}} \frac{\partial^2}{\partial \theta^2}{ (\ln f(x;\theta))} \,dx\\
&= \left[ \frac{\partial}{\partial \theta}{ (\ln f(x;\theta))}  \right]_{x = x_{min}}^{x_{max}}\\
&= \left[ \frac1{f(x;\theta)}  \frac{\partial f(x;\theta)}{\partial \theta}  \right]_{x = x_{min}}^{x_{max}}\\
&= \frac1{f(x_{max};\theta)}  \frac{\partial f(x;\theta)}{\partial \theta}_{x=x_{max}}
- \frac1{f(x_{min};\theta)}  \frac{\partial f(x;\theta)}{\partial \theta}_{x=x_{min}}\\
&= c,&c\text{ const.}\\
\end{align}$$


How to decipher this Microphone sensitivity Chart for Choosing a SPL meter

Microphones are typically specified in dB re 1 V/Pa. So am microphone with a sensitivity -44 dB (like this one http://www.panasonic.com/industrial/electronic-components/parametric-search.aspx?src=/www-cgi/jvcr13pz.cgi?E%2BSS%2B2%2BABA5010%2B0%2B4%2BAM) would produce 6.3mV RMS at 1 Pa of sound pressure. 1 Pa also corresponds to 94 dBSPL. So you can easily convert dB SPL in dB Pa by subtracting 94.

The datasheet for our microphone looks strange. I have never seen unit of 1.0V/0.1Pa before. This is potentially a print error.

The weighting curves do indeed mimic the frequency response of the human auditory system at different overall levels. However they are all 0 dB at 1 kHz, so if you calibrate at 1 kHz the weighting doesn't matter.


what does the output of ifft in ofdm look like?

You are correct that the output of the IFFT will normally have both a real and imaginary component.  That is fine and normal.  When the signal is modulated from baseband to the carrier frequency the real part is typically modulated by a sine wave and the imaginary part by a cosine (or vice versa).  This allows both parts to coexist without interfering with each other, since sines and cosines have zero cross-correlation.

Regarding how the numbers are represented digitally, they are represented the exact same way that real numbers are represented.  The only difference is that there is an extra set of numbers for the imaginary part of the signal.


Conditions for analytical signal reconstruction with the Hilbert Transform (HT)

Let me use the notation from Bedrosian's research report "The analytic signal representation of modulated waveforms" in a slightly simplified version:

$$s(t)=c(t)m(t)\tag{1}$$

where $s(t)$ is the modulated signal, $c(t)$ is a (not necessarily sinusoidal) carrier function, and $m(t)$ is the message signal. The carrier is assumed to be a relatively narrow-band band pass signal, and the message signal has a low pass character.

Let $c_a(t)$ be the analytic signal corresponding to $c(t)$. By definition, $c_a(t)$ has only positive frequency components. The analytic signal $s_a(t)$ corresponding to $s(t)$ can then be written as

$$s_a(t)=c_a(t)m(t)\tag{2}$$

if and only if the spectra of $c(t)$ of $m(t)$ do not overlap. Because only then does $s_a(t)$ as defined in $(2)$ contain positive frequency components only. If the spectra overlap, then $s_a(t)$ is not analytic because it will have negative frequency components. This is illustrated by the figure below, taken from Bedrosian's report:

The condition that the spectra of $m(t)$ and $c(t)$ do not overlap is equivalent to the requirement that $s(t)$ be a true band pass signal, i.e., $s(t)$ has no DC component.

I'm sure that this is what is alluded to by the sentence you quoted. I do agree that this is not obvious from the chosen formulation.


If there was some code outlining this, I would be interested

I don't have the expertise to describe pYin precisely, but I can point you to the original code for that paper https://code.soundsoftware.ac.uk/projects/pyin/repository
There's also a more recent implementation which made it into the librosa library since this question was asked. The code and pull request discussing the implementation are here: https://github.com/librosa/librosa/pull/1063/files#diff-f4245b657ccbb33a7bf56a5f7754585dfbbb80d2e545e840f678aceb548aa206R735


Next steps for analyzing this PWM algorithm and choosing a low pass filtering technique?

The pulse generation algorithm sounds like it is more like Pulse Density Modulation (PDM), not Pulse Width Modulation (PWM), because PWM would pack the ones into the beginning of counting, and leave the zeroes in the end.

Indeed some values of N and M cause the output pattern to have a certain repetitive bit stream, for example the 100/200 means the output is 50% 0 and 50% 1, so the output is a bitstream of 1,0,1,0... which is easy to filter. If you set N and M to be exactly ratio of 1/3, the output pattern will also be repetitive, 0,0,1,0,0,1... So, in light of this this also explains the spectrogram of N=50 M=200, which is 25%, 0,0,0,1,0,0,0,1.

These patterns that have a short repetition rate show up as pure tone peaks in the spectrogram. Choosing the values so that there is no obvious short repeating patterns basically dithers the bitstream and widens the spectrum to have multiple lower peaks.

Depending on what is the usage of the output will determine how to filter it. For many cases a simple RC filter is enough, as long as the PDM bitstream output rate is high enough for the RC filter cutoff frequency.


First express $x(t)$ in piece-wise form:

$$
x(t) =
\begin{cases}
0 ~&,  \text{ for}  & t  < -2 \\
0.5(t+2) &, \text{ for}  &-2 < t < 0 \\
-t+1 &, \text{ for}  & 0 < t < 1 \\
0 &, \text{ for }   &1 < t \\
\end{cases}
$$

Now, express $x_1(t)=x(2t)$ in a similar form:
$$
x_1(t) = x(2t) =
\begin{cases}
0 ~&,  \text{ for}  & 2t  < -2 \\
0.5(2t+2) &, \text{ for}  &-2 < 2t < 0 \\
-2t+1 &, \text{ for}  & 0 < 2t < 1 \\
0 &, \text{ for }   &1 < 2t \\
\end{cases}
$$

Recognize how the formula and intervals are adjusted for $x_1(t)$
$$
x_1(t) = x(2t) =
\begin{cases}
0 ~&,  \text{ for}  & t  < -1 \\
t+1 &, \text{ for}  &-1 < t < 0 \\
-2t+1 &, \text{ for}  & 0 < t < 0.5 \\
0 &, \text{ for }   &0.5 < t \\
\end{cases}
$$

Now, express $x_2(t)=x(2t+2)$ in a similar form:
$$
x_2(t) = x(2t+2) =
\begin{cases}
0 ~&,  \text{ for}  & 2t+2  < -2 \\
0.5(2t+2+2) &, \text{ for}  &-2 < 2t+2 < 0 \\
-(2t+2)+1 &, \text{ for}  & 0 < 2t+2 < 1 \\
0 &, \text{ for }   &1 < 2t+2 \\
\end{cases}
$$

Recognize how the formula and intervals are adjusted for $x_2(t)$ as well
$$
x_2(t) = x(2t+2) =
\begin{cases}
0 ~&,  \text{ for}  & t < -2 \\
t+2 &, \text{ for}  &-2 < t < -1 \\
-2t-1 &, \text{ for}  & -1 < t < -.50 \\
0 &, \text{ for }   &-0.5 < t \\
\end{cases}
$$

Now add those two piecewise defined signals $x_1(t)$ and $x_2(t)$, considering the intervals carefully:

$$ y(t) = x_1(t) + x_2(t) = x(2t) + x(2t + 2) =
\begin{cases}
0 ~&,  \text{ for}  & t < -2 \\
t+2 &, \text{ for}  &-2 < t < -1 \\
-t&, \text{ for}  & -1 < t < -0.5 \\
t+1 &, \text{ for }   &-0.5 < t < 0 \\
-2t+1 &, \text{ for}  & 0 < t < 0.5 \\
0 &, \text{ for }   &0.5 < t \\
\end{cases}
$$

Expressing this using unit steps (alone) is not possible, you need the ramp functions. Then you can provide a representation using scaled and shifted ramps.

A MATLAB plot of the resulting signal is provided below ;


using Short-Time Energy and zero crossing rate

Audio processing FFT of chunks of audio
This will give you the short term energy as a function of time and frequency.

The zero crossing in itself is easy to code (the example is conceptual, but very inefficient)

x0 = x(1);
for i = 2:length(x)
  x1 = x(i);
  if (x0 > 0 && x1 < 0) || (x0 < 0 && x1 > 0)
     % you have a zero crossing
     doSomethingWithAZeroCrossing(i);
  end
  x0 = x1;
end

However, zero crossings are fairly sensitive to benign and typically irrelevant changes in the audio signal, so zero crossings have only limited usefulness. In particular

even small amounts of noise can generate lots of extra zero
crossings and "jitter" the existing ones
Fluctuations in DC offset or very low frequency noise (like an Air conditioner or a heating fan picked up my a microphone) move the zero crossings all over the place


Is there a penis detection demo similar to face detection?

Is there a script / tutorial / demo for penis detection?
[...]
Fairly serious quesion, future of internet memes is at stake

Yes, there is.
Common Pattern Recognition techniques will be able to spot one even with what would be considered today "traditional" approach (i.e. without "Deep Learning").
There already is a sub-category of the subject with people working on "Adult Content Detection". Here are a few relevant publications.
This is an example for detecting child pornography which contains an overview and references to other techniques along with what the authors are proposing.
There is an increasingly large amount of tutorials out there on how to use a number of popular libraries to build an "Image Classifier", some of them not even requiring writing code. So, getting to a "script" that does what you are after is not a huge amount of work.
Hope this helps


How to perform an FFT on a signal with a sampling rate of 44100?

A 44100 Hz sample rate just means that you get 1 sample approximately every 22.67 microseconds. It also means you get N samples in N*22.67 microseconds, and you're free to choose N. If you choose N=44100, you get 44100 samples in one second. If you choose N=1024, you get 1024 samples in 23.22 milliseconds.


DFT of 2d real signal and Hermitian symmetry

Yes, multidimensional DFT's have complex conjugate symmetry in each dimension

I think you might be misinterpreting the FFTW manual which is referring to internal storage formats.

Edit:

  why the transformed vector in FFTW is diminished only by half (only
  one dimension values are stored on N/2+1 complex numbers)? Why isn't
  each direction diminished?

From the FFTW manual:

  The multi-dimensional transforms of FFTW [...] compute simply the separable product of the given 1d transform along each dimension of the array

So for example if you provide a 256x256 array of real numbers the result will be 256 1 dimensional DFT's. As they mention at the bottom of that documentation this is not a multi-dimensional DFT.

see http://www.fftw.org/fftw3_doc/Complex-Multi_002dDimensional-DFTs.html


Does a simple photograph contain more information than a complex painting?

It depends how you define the term "information" or "entropy".

The conventional definition of entropy of an image is to think the image as a two-dimensional matrix of pixels and
$$H = - \sum_k p_k \log_2(p_k),$$
where $p_k$ is the probability, which is calculated from histogram, associated with gray level $k$.

This kind of entropy is correct if we ignore the correlation between pixels. For example the two images have the same entropy by this definition.

It is not true if the correlation between pixels is considered. For example if the color first pixel in top-left has probability $p_k$, the next pixel surely has the same color and its color does not have the same probability $p_k$.

We human being, with you as an example, use this kind of correlation to perceive the images. This kind of correlation are called "details" and, therefore, we/you expect that rich-detail images should have more information/entropy than simple ones. This is the reason why you have found it counter-intuitive.

PS:

I tried calculating the entropy of the two images you have posted, but they are not different "by one or two orders of magnitude" !!!!

"Along the River During the Qingming Festival" entropy about 7

"The sky" entropy about 6

They must not be the same files of the professor.


Are there any quick way to tell if a kernel is a high-pass filter?

Proper way: do a Fourier transform and look at the magnitude of the result
Quick and dirty: take the sum of all all elements. If it's zero or small compared to the absolute sum, it has some high pass-ness to it.


Image cleaning - How to remove line art from a colored image?

I have a couple of solutions, but I have not test either of them for you image.

I. Use simple neighrest neighber interpolation, i.e.

sweep through whole pixels,
Check if the ith pixel is black (threshold it, or get its HSV transform,
whatever)
if it's black, substitute colour values of its nearest non-black
neighborer, in the pixel.

II. Using Medial filter,

Simple apply a medial filter on the image, but do not count black pixels when applying on each window of image.


T wave analyis and extraction using discrete wavelets

First, this is not an easy task to automate. An algorithm that works great on clean ECGs can fail miserably on noisy ones. That said, I would suggest you look at the approach used to identify the ST segments in the Long Term ST Database. It is one of, if not the, largest ST segment database publicly available. They describe the approach they used, which includes automatic segmentation followed by expert review.

Also, take a look at the PhysioToolkit from PhysioNet. I'm pretty sure the have some automatic ECG annotation tools. I looked, but couldn't find them. Still, I think they are there somewhere.


How to Emulate Sample Frequency Offset for OFDM in Matlab?

You effectively want a interpolated waveform that is interpolated by 200,000 such that for each new sample with the 5 ppm offset you can select one additional offset to induce that time offset, for example

x[1], x[200,002], x[400,003], x[600,004]....

(Or equally one less if your sampling clock increased in frequency 5ppm).

One very simple approach to do this is to use linear interpolation to determine the value of each new sample given the adjacent two samples, or curve fitting to more adjacent samples using higher order polynomials which is what is done with detailed filter design approaches. If your waveform is sufficiently oversampled to start with, then a simple linear interpolation approach may be more than sufficient. (To assess this you could measure the error vector magnitude of the mid sample using linear interpolation compared to a resample by 2 and see if that error is less than your target SNR requirement).

With a 5ppm lower clock the new samples would be:

$$y[n]=x[n] + mod(n,200,000)(x[n+1]-x[n])/200000$$

For n as the index of the original N sample waveform 0 to N-1

For example, If the first 4 samples were

2, 5, 8, 11

First sample: 2

Second sample: $5+ 1(8-5)/200,000$

Third sample: $8+ 2(11-8)/200,000$


How do you measure the accuracy of a stereo system?

From a more pragmatic than philosophical perspective, possibilities include:

Comparing with ground truth data--these data could come from surveying's instruments, or be assumed to be "reliable" (e.g. you could assess the system's accuracy on the Eiffel tower and other well-known buildings for which measurements have been curated);

Verifying physical properties are met (e.g., floors are horizontal, walls are perpendicular, distances between floors or between windows are constant);

Leaving one ore more views out of the 3D estimation process and assessing how far the projections of the estimated 3D points are from the actual 2D positions in these views (ie., jackknifing).


What do these values in the freeverb diagram mean?

It's a Lowpass-Feedback Comb Filter, so let's dissect this one at a time:

Comb filter: You have delay line and add the output back to the original signal. That gives you a single echo
Feedback: You take some of the output of the delay line and put it back into the input of the delay line. That gives you a series of echo's that repeat regularly. You need to put a gain < 1 on this, otherwise it will repeat forever. That gain is called the feedback gain and controls how fast (or slow) the echoes decay.
You put a lowpass into the feedback path to model the typical absorption in a room that's higher at high frequencies and lower at low frequencies. You basically make the feedback gain frequency dependent.

$d$ is the pole location of the lowpass filter. That determines the cutoff frequency of your lowpass filter

$f$ is the feedback gain and determines how quickly the repeating echoes decay. That should be a function of the reverb time (at low frequencies) and the delay size


Filter Cascade: Additions and Multiplications per input sample

Their numbers work out right if you do the stated decimation before the stage, and you pay attention to the fact that decimations are cumulative.
So stage one operates at 1/3 of the input rate, stage two operates at 1/6, stage three operates at 1/12, and stage four operates at 1/48.  If you scale the adds and multiplies by those numbers, then their stated numbers work.


How to take acoustic measurements of a clasroom

It doesn't matter where your microphone is in the room when you take the measurement. What you need to do is find where the peaks are in the cross correlation (each peak corresponds to an echo at a certain lag/delay), figure out what the corresponding inverse filter needs to be to cancel out the sound.

If you're trying to do echo cancellation in MATLAB and have the signal processing toolbox, take a look at the example in xcorr by typing:

  doc xcorr

for a simple example.

Edit 1:To find the inverse filter required, you need to calculate the transfer function between the source and the mic. (edit: sentence in here was misleading) Once your filters are set up, as long as nothing is changing you should be able to cancel the echoes produced by the room.

Edit 2: To get the room impulse response you typically use a long LFM chirp (some people pop a balloon) and a microphone, then calculate the transfer function between the input signal and the recorded signal. Here's an example paper of what you can do: http://www.artalabs.hr/papers/im-aaaa2003.pdf

Edit 3:

  You do realise that reverberation (or echoes as you described) are dependent on location and even time?

you said you are making an "impulse response invariant" filter. So as long as neither your source or microphone are moving, the IIR filter created from your Room Impulse Response measurement will cancel out your echoes...just make sure the mic and speaker stay in the same spot.


A Hilbert Transformer is an acausal linear time-invariant (LTI) filter.  Being LTI it has all of the things an LTI filter has: the output is the convolution of the input against an impulse response.  And it has a frequency response with magnitude and phase.

But being acausal, that means that the impulse response begins to respond before the impulse hits the input (i.e. it can see future input).  That, of course, is impossible in a real-time system so the impulse response of the Hilbert Transformer is delayed to make it causal.  That delay causes a linear trend to the phase component.

The magnitude frequency response of the Hilbert Transformer is 1 (as linear gain) or 0 dB for virtually all frequencies (DC and Nyquist are the exceptions). But without the delay, a Hilbert Transformer shifts the phase for all positive frequencies by -90° (or $-\frac{\pi}2$ radians) and all negative frequencies by +90°.  That puts a phase discontinuity at 0 Hz (DC).  Likewise, you will see the opposite phase discontinuity at Nyquist.


Duality of the continuous-time Fourier transform - derivation and notation

You've answered your first question yourself: if you use the notation $X(j\omega)$ then you'll need to use $X(jt)$, because you simply replace $\omega$ by $t$.

The sign error in your derivation comes from the fact that you forgot to change the sign of the limits of the integral. After substituting $\omega=-u$ you should get

$$X(t) = -\int_{\infty}^{-\infty} x(-\omega) e^{j\omega t} \mathrm{d}\omega = \int_{-\infty}^{\infty} x(-\omega) e^{j\omega t} \mathrm{d}\omega$$


People interpret the term "resolution" more than one way.  One way is the ability to separate closely spaced frequency lines, and another is to estimate the the "true" frequency of a single tone.  Most DFT interpolators fall into the true value of a tone category.  The problems are not independent of each other but you should understand the problem you are trying to solve.

There is more than one way to interpolate intermediate frequencies from the DFT.

From Peter K's web page, the Fourier Coefficient section links to a survey prepared by Eric Jacobsen.

I've used

Macleod, Malcolm D. "Fast nearly ML estimation of the parameters of real or complex single tones or resolved multiple tones." IEEE Transactions on Signal processing 46.1 (1998): 141-148.

which is a complex interpolator with good result

The basic idea is that there is a largest bin that contains the frequency of interest and the bin leakage is asymmetric. One uses a "fit" to some formula as an interpolation.  THis implies that you need a nominal "high" SNR, and a few bins of separation. Macleod's technique assumes that you used a boxcar window which can be a problem where the noise is far from flat.

You might want to post plots for your problem for further advice.

The other categories  on Peter K's page are also worth looking at.


Butterworth filter's gain formula does not agree with R's `signal` package

The function butter() computes the coefficients of a discrete-time ("digital") Butterworth filter, whereas the gain formula you used is valid for a continuous-time ("analog") Butterworth filter. According to the R documentation you can use butter() to compute an analog filter using the plane argument.
The gain of a discrete-time Butterworth filter is obtained by the bilinear transform, which substitutes the "analog" frequency variable by $\tan(\omega/2)$, where $\omega$ is the frequency in radians, normalized by the sampling frequency. For a discrete-time unit gain lowpass filter of order $n$, the squared gain is given by
$$G^2(\omega)=\frac{1}{1+\left(\frac{\tan(\omega/2)}{\tan(\omega_c/2)}\right)^{2n}},\qquad |\omega|\le\pi\tag{1}$$


How are pulse trains generated in a radio transmitter?

Arbitrary Waveform Generators are a special type of test equipment that can generate, as the name implies, arbitrary waveforms as output signals that can be fed into transmitters or other circuits. They vary significantly in sophistication.


Relation between Laplace and Fourier transforms

The Laplace transform evaluated at $s=j\omega$ is equal to the Fourier transform if its region of convergence (ROC) contains the imaginary axis. This is also true for the bilateral (two-sided) Laplace transform, so the function need not be one-sided.

As for real and imaginary parts, since $s$ is a complex variable, both the Laplace and the Fourier transform generally have real and imaginary parts. Take as a simple example the function $x(t)=e^{-at}u(t)$, with $a>0$, where $u(t)$ is the unit step function. The Laplace transform is

$$X_L(s)=\frac{1}{s+a}\tag{1}$$

Since $a>0$, the ROC of $X_L(s)$ contains the imaginary axis, and the Fourier transform of $x(t)$ is simply obtained by evaluating $X_L(s)$ on the imaginary axis $s=j\omega$:

$$X_F(\omega)=X_L(j\omega)=\frac{1}{j\omega+a}\tag{2}$$

Since $s=\sigma+j\omega$ is generally complex, not only the Fourier transform but also the Laplace transform $(1)$ has a real and an imaginary part:

$$X_L(\sigma+j\omega)=\frac{1}{\sigma+j\omega+a}=\frac{\sigma+a}{(\sigma+a)^2+\omega^2}-j\frac{\omega}{(\sigma+a)^2+\omega^2}\tag{3}$$

Only when evaluated on the real axis $s=\sigma$ ($\omega=0$) does the imaginary part vanish.


Decomposition of $H(z)$ as maximum-phase, minimum-phase

So let's check that:
$$H(z) = H_{\rm min}(z) H_{\rm max}(z) H_{\rm uc}(z) \tag{1}$$
where:
$$
H_{\rm min}(z) = \frac{z-1/2}{z^3}\tag{2}
$$
and
$$H_{\rm max}(z) = z-2\tag{3}$$
and
$$H_{\rm uc}(z) = z-1\tag{4}$$

Substituting $(2)$, $(3)$, and $(4)$ into $(1)$:

\begin{align}
H_{\rm min}(z) H_{\rm max}(z) H_{\rm uc}(z) &= \frac{z-1/2}{z^3} (z-2) (z-1)\\
 &= \frac{z-1/2}{z^3} \left(z^2 - 3z + 2\right)\\
&= \frac{1}{z^3} \left(z^3 - 3z^2 + 2z - \frac{1}{2}\left(z^2 - 3z + 2\right)\right) \\
&= \frac{1}{z^3} \left( z^3 - \frac{7}{2} z^2 + 5 z - 1\right)  \\
&\not= H(z)
\end{align}

So something is wrong!

Try:
$$
H_{\rm min}(z) = 1-1/2z^{-1}\tag{2A}
$$
and
$$H_{\rm max}(z) = 1-2z^{-1}\tag{3A}$$
and
$$H_{\rm uc}(z) = 2(1- z^{-1})\tag{4A}$$

Now:
\begin{align}
H_{\rm min}(z) H_{\rm max}(z) H_{\rm uc}(z) &= \left(1-1/2z^{-1}\right) \left(1-2z^{-1}\right) \left(1-z^{-1}\right) 2\\
&= \left(1 - \frac{5}{2} z^{-1} + z^{-2}\right) \left(1-z^{-1}\right) 2\\
&= 2 - 5 z^{-1} + 2z^{-2} - \left(2z^{-1} - 5 z^{-2} + 2 z^{-3}\right)\\
&= 2 - 7 z^{-1} + 7 z^{-2} - 2z^{-3}\\
&= H(z)
\end{align}

I believe your interpretation of Oppenheim and Schafer (OS) is incorrect. And I can find no reference that says the poles have to be outside the unit circle for a maximum phase system. This one and the one I referenced in the comments both only mention the zero locations.

And this is what my copy of OS says:

which also only mentions zero locations.

I believe exercise 5.63 is in error.  The definition used in the body of the book is:

which again does not mention pole locations.

Also, having
$$
H_{\rm max}(z) = H_{\rm min}\left(z^{-1}\right) z^{-M_i}
$$
ensures that $H_{\rm max}$ is causal for FIR $H_{\rm min}$.

So you have two options: either the exercise 5.69 definition is wrong, or the relationship between $H_{\rm max}$ and $H_{\rm min}$ is wrong.


The standard meaning of white noise includes an insistence (whether implicit or explicit) that the mean is $0$. Thus, what you want to prove is trivially true: since
$$Y[n] = \sum_{k=-\infty}^\infty h[n-k]X[k] = h\star X \big\vert_{n},$$ the linearity of expectation (the notion that $E[aX]=aE[X]$ and that the expectation of a sum is the sum of the expectations) gives us that
$$\mu_Y[n] = E\left[\sum_{k=-\infty}^\infty h[n-k]X[k]\right]
= \sum_{k=-\infty}^\infty h[n-k]E\big[X[k]\big]$$
where all the $E\big[X[k]\big]$ have value $0$, and so both sides have value $0$.

More generally, if the $X$ process has constant mean $\mu_X$, then
$$\mu_Y[n] = \sum_{k=-\infty}^\infty h[n-k]E\big[X[k]\big]
= \sum_{k=-\infty}^\infty h[n-k]\mu_X = \mu_X\sum_{k=-\infty}^\infty h[n-k]
= \mu_X \sum_{i=-\infty}^\infty h[i]$$ is also a constant.


Image Processing/ Computer Vision: How is separability of filters implemented?

Szeliski isn't saying what you think they're saying.

The vertical part of the separable filter is:
$$
\mathbf{v_2} = 1/4\begin{bmatrix}1\\2\\1\end{bmatrix}
$$
not
$$
\mathbf{v_1} = 1/2\begin{bmatrix}-1\\0\\1\end{bmatrix}.
$$

You can see this by defining:
$$
\mathbf{h}= 1/2\begin{bmatrix}-1&0&1\end{bmatrix}
$$
and looking at

$$ \mathbf{K_1} = \mathbf{v_1} \mathbf{h} = \begin{bmatrix} 1/4 & 0 & -1/4\\
0 & 0 & 0\\
-1/4 & 0 & 1/4
  \end{bmatrix}$$

and

$$ \mathbf{K_2} = \mathbf{v_2} \mathbf{h}= \begin{bmatrix} -1/8 & 0 & 1/8\\
-1/4 & 0 & 1/4\\
-1/8 & 0 & 1/8
  \end{bmatrix}$$

You can see this effect in the image I included in your question: the vertical edges are highlighted, but the horizontal ones disappear. This is because the $\mathbf{v_2}$ filter is a low-pass filter, not enhancing edges.


Does down-mixing a signal before taking the fourier transform aid in signal detection?

If you are doing are doing this in both the analog and digital domain, then there may be significant differences due to trade-offs between the noise sources in the analog stage(s) in the filters and the mixer oscillator frequency, and the corresponding sampling jitter and quantization noise in converting to digital.

Even strictly in the digital domain, the bandpass and lowpass filter responses may have different effects on how potential stop-band noise might affect your detection accuracy, compared to how windowing (rectangular or otherwise) affects the continuous frequency response of each FFT bin.


The sequences $e^{j\omega_0n}$ are eigensequences of discrete-time LTI systems, i.e., the response to such a sequence is the same sequence scaled by a complex constant (the eigenvalue). This can be shown as follows:

$$\begin{align}y[n]&=(h\star x)[n]\\&=\sum_{k=-\infty}^{\infty}h[k]x[n-k]\\&=\sum_{k=-\infty}^{\infty}h[k]e^{j\omega_0(n-k)}\\&=e^{j\omega_0n}\sum_{k=-\infty}^{\infty}h[k]e^{-j\omega_0k}\\&=e^{j\omega_0n}H(\omega_0)\end{align}$$

where $h[n]$ is the system's impulse response, and $H(\omega)$ is its frequency response. A completely analogous proof can be given for continuous-time systems.


Classifying animated GIFs as animation or film

It's a studied problem. Here are some relevant papers:

Detecting cartoons: a case study in automatic video-genre classification
A color-action perceptual approach to the classification of animated movies
Detecting Cartoons: Automatic Video Genre Classification


How to determine the "variability" in a frequency range?

Daniel,

Upon re-reading your question, it appears that what I have learned to be known as the 'Gabor Bandwidth" might be useful to you in this case, for you trying to measure 'spectral variability'. (Dilip provided a good answer to my question on Spectral Moments here).

When I studied it further, the Gabor-Bandwidth seems to really just be a measure of how 'spread' the spectrum is from its mean. (Hence the manipulation of moments).

Take a look and see what you think.


Explicit a succession using z inverse transform

If you are allowed to define your system to be causal so that it is non-zero only for $n \ge 0$, then you can use the unit step function $u[n]$ so that your system becomes:

$$y[0] = 1$$
$$y[n+1] = 2y[n] + 3u[n]$$

Which does have a Z-transform representation.


Downsampling oversampled image (smallest size which retains image 'detail' within acceptable limits)

There are many options yet it is highly dependent on the properties of the image being resampled.  Images downsampled with bicubic interpolation are regarded as smoother and with fewer interpolation artifacts. You may also try nearest neighbor interpolation and bilinear interpolation. Make sure a low pass filter is applied before further processing with downsampling strategies, otherwise the aliasing artifacts may be induced.

Besides, you may also consider vectorization algorithm  in which a resolution independent vector representation is created before rendered to the image at a desired resolution. it is mainly used in up-sampling but the vectorization method is helpful in preserving the feature connectivity in down-sampling as well.

Another interesting blog on comparison between different image re-size approaches.

EDIT

Regarding the design of low pass filter served as a pre-processing method before downsampling, you can consider this way:

Suppose a low pass filter $g$ works on the image $x$, and the filtered image is $u = x*g$. After that, simply decimate $u$ by 2, you get $v$. Now you want to recover the original image. You can zero-upsample the $v$ into $w$ where $w_{2n} = v_{n}$, and $w_{2n-1} = 0$. Then another filter $h$ is applied, trying to obtain $\hat{x} = w * h$. After Z-transform, $U(z) = X(z)G(z)$, $\hat{X}(z) = W(z)H(z)$.

Since $W(z) = 1/2(U(z) + U(-z))$, we have $\hat{X}(z)=1/2(H(z)G(z))X(z) + 1/2(H(z)G(-z))X(-z)$

A good downsampling is equivalent to a good recovery to some extent. In order to best recover the original image ($\hat{X}(z)=X(z)$), you are trying to make $H(z)G(z)$ as close to $2$ as possible, while $H(z)G(-z)$ as close to $0$ as possible. The design of low pass filter $G(z)$ is based on this criteria.


Generally, this should not be possible because:

you project a 3D information onto a 1D signal, and the system is underdetermined in general,
you add a non-linearity (the norm), that further hinder restoration.

However, using the triangle inequality, you can bound $B$ from below:

$$ \|A+B\| \le \|A\|+\|B\|$$

hence

$$ B_{\inf} =\max(\|A+B\| -\|A\|,0)\le \|B\|$$

In the example you gave, if the signals $A$ and $B$ are weakly related, outcomes can become interesting:

signal = (sum(C.^2)).^(0.5);
signalInf = max((signal-(sum(A.^2)).^(0.5)),0);
signalG = (sum(B.^2)).^(0.5);

plot(t,[signal',signalInf',signalG'])
title('Observed signal');
legend('A','Binf','B')

From here, you can recover putative peaks from $ B_{\inf}$. However, it is unlikely to work every time.

Could additional filtering help on $S(t)$? Without knowledge about the relative spectra of $A$ or $B$, I don't know (yet).


fft and inverse fft error size and applying impulse response matlab

I don't know your overall goal but if you want to test the MSE associated with FFT / IFFT then you should perform the following

num_of_samples = 2819519 ;
fft_len         = 2^nextpow2(2819519) ;

X = fft(signal, fft_len);
y = real( ifft(X, fft_len) );
mse = sum( (signal - y(1:num_of_samples) ).^2 )/num_of_samples


I will be using MATLAB as part of my solution.  The basic algorithm will be this:

Read in the image and convert to black and white
Invert the intensities and fill in the holes.  There are speckles in the image where when you convert the image to black and white, the objects are not solid.  We want to ensure these are solid
Count how many unique objects there are using bwlabel.

You'll need the Image Processing Toolbox for this algorithm.  If you don't have this and are using MATLAB, let me know and I'll edit my post.

Without further ado:

% Read in image, convert to black and white - Link comes from your image posted here
im = imread('http://i.stack.imgur.com/lBGU1.png');
imBW = im2bw(im, 0.3); %// Specify manual threshold of 0.3

% Invert intensities and fill in holes
imBWFilled = imfill(~imBW, 'holes');

% Count how many unique objects there are
[L,num] = bwlabel(imBWFilled);

% Show final image and display number of objects counted in the title
imshow(imBWFilled);
title(['Total number of objects: ' num2str(num)]);

L contains a map where each pixel is an ID of which object that belongs to.  0 means that the pixel belongs to the background while any pixel that has values greater than or equal to 1 means that the pixel belongs to a particular object associated with that ID number.  num gives you the total number of objects seen in the image.  The output of this code thus gives:

If you don't have MATLAB and want to compute this using another language, let me know.  You can find the number of objects by considering an image as a connected graph.  If you don't have bwlabel, you can use any graph searching algorithm (breadth-first search, depth-first search, etc.) to help you compute this.  Start with any pixel that belongs to an object, and perform BFS / DFS to visit those pixels that belong to the object.  You then set all of these pixels to belong to an ID number.  When the queue / stack is empty, you then select another pixel that belongs to an object and repeat the algorithm.  You stop when you have visited all pixels that belong to objects.  The total number of IDs you have issued will essentially be how many unique objects you have in your image.

Edit - July 14th, 2014

Seeing your comments, you wish to be able to distinguish between what is a nut and what is a bolt.  As such, we can simply add on top of the current code.  The algorithm that I have developed is loosely based on circularity:

For each object, find the centre of mass.  This is simply taking all of the X and Y co-ordinates and averaging them.
For each object, find their boundaries / perimeters
For each centre of mass we have, calculate the distance between this point and all of the object boundary points.
Find the difference between the maximum and minimum of these distances - denote this the range
We will have N ranges for N objects.  Simply take a look at these ranges and threshold the array.  If you take a look at the bolts, they have a longer vertical direction than horizontal direction.  For the bolts, the vertical and horizontal direction is pretty much the same.  As such, the range should give us an idea of what is a nut and bolt because the difference between the maximum distance to the centre and the minimum distances should be very large in comparison to the nut.  As such, for each object we have, check to see if the range is below a certain value (nut) or above a certain value (bolt).

Without further ado, here is the code:

%%% For each object in the image, find the centre of mass
centres = zeros(num,2);
% Cycle through each unique object label and extract (X,Y) co-ordinates
% that belong to each object.  Compute centre of mass for each.
for n = 1 : num
    bmap = L == n;
    [rows,cols] = find(bmap == 1);
    centres(n,:) = [mean(cols) mean(rows)];
end

% Find boundaries of all objects
bwBound = bwperim(imBWFilled, 8);

% For each object, find the distances between the centre of mass with all
% of the pixels along the boundary for each object.  Find the range (max -
% min).
ranges = zeros(num,1);
for n = 1 : num
    bmap = L == n; % Obtain all pixels for an object
    boundPix = bwBound & bmap; % Logical AND with boundaries map to extract
                               % only those pixels around the perimeter
    [rows,cols] = find(boundPix == 1); % Find these locations
    % Compute the distances between the centre of mass with these points
    dists = sqrt((cols - centres(n,1)).^2 + (rows - centres(n,2)).^2);

    % Find the difference between the maximum and minimum distances
    ranges(n) = max(dists(:)) - min(dists(:));
end

This isn't complete yet.  I stopped here so you can see what the ranges for all of the objects look like:

ranges =

37.9615
63.9613
54.9266
 5.0716
 4.1578
 6.3114
 7.2356
41.6381
10.3123
34.0938
 5.0021
67.3290

As you can see, there are some ranges that are quite small ($< 10$) while there are some ranges that are large.  To be safe, let's choose a threshold of 15.  As such, for those ranges that are less than 15, this is classified as a nut while those that are larger are classified as a bolt.  When we're done, we simply count how many of these fall within each range, and those are the number of nuts and bolts we have.  Let's continue with the algorithm:

% Find those object IDs that have less than a range of 15.  These are the
% bolts
% The rest are nuts
indBolts = find(ranges < 15);
indNuts =  find(ranges >= 15);

% Total number of nuts and bolts
numBolts = numel(indBolts);
numNuts = numel(indNuts);

This will give you 6 and 6 respectively as you expected.  Now to finish up the algorithm, I'm going to give you an added bonus.  For each unique object we have, I'm going to colour their interiors to be a certain shade of gray.  The background will be black, while a nut is gray and a bolt is white.  The code to do this is:

finalMap = uint8(zeros(size(imBW)));

for n = 1 : numBolts
    finalMap(L == indBolts(n)) = 128;
end
for n = 1 : numNuts
    finalMap(L == indNuts(n)) = 255;
end

figure;
imshow(finalMap);
title(['Number of Nuts: ' num2str(numNuts) ', Number of Bolts: ' num2str(numBolts)]);

The final image we get is:


Synchronizing data sent over sound beeps with Goertzel Detection, C#

You can't simply skip ahead some amount of time, and won't be able to skip ahead when decoding in real time.

You need to simply decode as the values come in, and switch your output depending on what comes in.  Change states only on input changes.
So, like this:

Feed input signal into the bank of goertzel filters,  the detected tone is the filter output that goes over a defined limit.
When detected tone changes to 12kHz, set message start flag
When detected tone changes to 13khz, send output bit 1
When detected tone changes to 14kHz, send output bit 0
When detected tone changes to 10kHz, end of message.
Repeat from step 1 until end of message.

Working that way, you should stay synchronized and have fewer errors.  You should also be able to shorten the beep lengths, though you may need to work with the length of the Goertzel filters to make it work with shorter beeps.

You will find that you will have to make a compromise between the length of the beeps, the length of the goertzel filter, and the separation between the individual tones in order to get fast, reliable communication.

To get faster communication, you will need shorter filters which in turn require larger separation between the tones to get reliable detection.

What you are doing is called frequency shift keying (FSK,) and has been used to transmit data over audio paths for a long time.

Noise will always be a problem.  There are two ways to fight it, and you might have to use them both.
First, you make your goertzel filters as long as possible, but still short enough to react within your desired bit rate.
Second, most systems encode the data with checksums.  Checksum codes can be used that allow you to recover from bit errors - they don't just say "bad data," they can also say "bad data, bit X should be value Z."  There are of course limits.  The better the checksum and error correction, the more extra information you have to send, so the slower your effective data rate.

A third thing you could do is to lower the data rate and extend the length of the goertzel filters, but you are already running a very low rate.


Combining 2 sinusoids of equal amplitude with different frequencies into 1 new "wave"

Let us assume the following signal form $$x(t)=\sin(2\pi t f(t))$$ where $f(t)$ is some "instantaneuous frequency", which should increase linearly from $w_1$ to $w_2$ over the time.

The time should end, when the sine wave has done one period, i.e. its argument became $2\pi$. At this time, $f(t)=w_2$. From that we can calculate:

$$2\pi w_2 t_e=2\pi\\t_e=1/w_2,$$ where $t_e$ is the end time of the simulation. Hence, we have $$f(w)=w_1+\frac{w_2-w_1}{t_e}t=w_1+(w_2-w_1)w_2t.$$ We can plot this function:

w1 = 5
w2 = 10
T = 1/w2

Fs = 100000

t = np.arange(0, T, 1/Fs)

f = np.linspace(w1, w2, len(t))
plt.plot(t, np.sin(2*np.pi*f*t))
plt.title('Graph')
plt.grid(True)

Now, regarding the question, where the zero is. The sine goes to zero, when its argument becomes $\pi$. Hence, we have

$$2\pi (w_1+(w_2-w_1)w_2t)t=\pi,$$

which we can solve by (its a quadratic equation, so nothing special, but I use mathematica for convenience):

$$\text{Solve}\left[t (t \text{w2} (\text{w2}-\text{w1})+\text{w1})=\frac{1}{2},t\right]$$

$$\left\{\left\{t\to \frac{\sqrt{\text{w1}^2-2 \text{w1} \text{w2}+2 \text{w2}^2}-\text{w1}}{2 \left(\text{w2}^2-\text{w1} \text{w2}\right)}\right\},\left\{t\to \frac{\sqrt{\text{w1}^2-2 \text{w1} \text{w2}+2 \text{w2}^2}+\text{w1}}{2 \left(\text{w1} \text{w2}-\text{w2}^2\right)}\right\}\right\}$$

There are two solutions, because for negative t, we can also end up with a solution (both $f(t)$ and $t$ would be negative to yield 1/2 as the argument).

To answer the comment from @Jason R: This indeed can also be seen as a frequency-modulation, given by

$$x(t)=\sin\left(2\pi\int_0^t\phi(\tau)d\tau\right)$$

when setting

$$\phi(\tau)=w_1+2(w_2-w_1)w_2\tau.$$

In this case we have $tf(t)=\int_0^t\phi(\tau)d\tau$.


I would recommend doing some proper reading of books/tutorials etc on the Fourier transform and the Discrete Fourier Transform (DFT).  You will also want to look at filters and probably convolution for the bandpass filter.  But regardless, here's some info that should help with your current state:

Use matplotlib to plot the FT that you've calculated:

fs = 44100 # Or whatever the actual sample rate is (Hz)
f = linspace (fs,len(data), endpoint=False)
mpl.plot (f, abs (fftres))
title ('Magnitude spectrum of the signal')
xlabel ('Frequency (Hz)')

Another function you will want to look at is fftshift, which makes it simple to 'roll' your spectrum around so that DC is in the centre.

The frequencies that are sampled by the DFT are between $-f_s/2 < f < f_s/2$ for sampling frequency $f_s$.  And the FFT produces the same number of samples as the number of data points provided to it.  The FFT function returns the spectrum so that the DC (constant/average) value is first in the array.  Use the function fftshift to adjust the array if you want it ordered  $-f_s/2 < f < f_s/2$.

I will now repeat: Fourier transforms and filtering are a huge field, and while it will definitely be interesting and useful to play with the functions for a bit like you are doing, you really need to read up on what is going on if you want to actually do anything useful.

Aside: why are you starting at index 1 in your loop, rather than 0?


Alternative to Orthogonal Matching Pursuit (OMP) Algorithm

The main advantage of OMP is that the residual is orthogonal to the current solution.

Let's say you select all $k$ columns from $A$ (also called atoms) at once and let us also presume that $A$ is an overcomplete basis (this is more or less the standard in OMP literature).

Now, with your method, if the atom that correlates the most with your measurements $y$ is linear-dependent with $p < k$ other atoms in $A$ you will end-up with an $k-p$-sparse signal, because $p$ entries will be more or less redundant. The same argument can of course be extended to less correlated atoms. You might also be lucky and never see the phenomenon.

Let's take the same example but with OMP this time. During the first iteration you would select the atom that correlates the most with measurements $y$. After that you compute the coefficient in $x$ such that the new residual is orthogonal to the current measurements approximation. In other words you got the most of the information provided by the currently selected atom so during the next iteration you are very likely to pick an atom that contains fresh information (ask yourself what would happen with linear dependent atoms in this case).

Here is a list of atom selection look-ahead strategies based on OMP and OLS that you might find interesting to read: POMP, LAOLS and POLS.


Filtering performance on Poisson noise with quadratic data-fidelity

For large intensities / large "bins", i.e. "areas for which events are counted and accumulated", Poisson processes lead to nearly Gaussian distributed individual values -- basically, without trying to derive this, I think that's  the application of the CLT on a lot of realization of a point process.

EDIT: Shameless plug: I really really like the wikipedia page on "Relationships among probability distributions; notice the arrow between the Poisson and the normal distribution:


Matlab - How to create an autocorrelogram using a spike-train

When I run that line of code, however, I find that there's a peak towards the length of the timestamps vector.

So I'm assuming you're running:

r = xcorr(x);

and you're perplexed by the peak at r(length(x)) ?  That's because xcorr returns a vector that is the autocorrelation from a lag of -length(x) to a lag of +length(x) (unless told otherwise) so that the middle value is the zero lag.

However you seem interested in finding correlations between the pulse times.  That's a whole other kettle of fish... especially when there are pulse trains of different periods or there are missing data. To see the complexity that can happen, check out this paper.

What are you trying to achieve?

To see that there's a big difference between the autocorrelation function of the pulse arrival times and the actual pulses themselves, I've generated an example in R.  The resulting autocorrelation functions are plotted below.

--

Ok, so you have a periodic(ish) pulse train and you want to figure out what the pulse period is.

The autocorrelation of the pulse times will not help. This is because the times are just a sorted array, so you'll just get a straight line as the autocorrelation.
The autocorrelation of the pulses themselves should help.  However, you'll get a peak at 0 lag, and at lags of $nP$ where $P$ is your pulse period and $n$ is an integer.
The post you reference is looking at the histogram of the inter-pulse times (time differences between successive pulses).  Depending on how varied your pulse train is, this might be the best approach.

If I change my simulation to use a pulse that is periodic (with zero noise on the pulse times!), then I get the second lot of code below, and the results are graphed here:

For this particular run, the Tperiod = 9.350267.  As you can see, both the ACF of the pulses and the histogram method get a peak close to where we want it (though quantization effects mean we don't get exactly at 9.350267; better selection of "bins" might improve that).

R Code Below

#26547

Npulses <- 100
Tmax <- 1000

pulseTimes <- sort(runif(Npulses,0,Tmax))

# install.packages("Hmisc")
library(Hmisc)

xc_times = acf(pulseTimes, plot = FALSE)

times <- 0:.1:Tmax

pulses <- times*0

for (pulseNo in 1:length(pulseTimes))
{
  idx <- which.min(abs(pulseTimes[pulseNo] - times))
  pulses[idx] = 1
}

xc_pulses = acf(pulses, plot = FALSE)

plot(xc_times$lag,xc_times$acf, type="l", ylim=c(-0.1,1))
lines(xc_pulses$lag, xc_pulses$acf, col="red")
legend(15, 0.9, c("ACF of times", "ACF of pulses"), lwd=c(2.5,2.5),col=c("black","red"))

Second version of R code below

#26547

Npulses <- 100
Tmax <- 1000

Tperiod <- runif(1,9,10)
Tstart <- runif(1,0,1)

#pulseTimes <- sort(runif(Npulses,0,Tmax))
pulseTimes <- Tstart + (0:Npulses)*Tperiod

# install.packages("Hmisc")
library(Hmisc)

xc_times = acf(pulseTimes, plot = FALSE)

times <- 0:.1:Tmax

pulses <- times*0

for (pulseNo in 1:length(pulseTimes))
{
  idx <- which.min(abs(pulseTimes[pulseNo] - times))
  pulses[idx] = 1
}

xc_pulses <- acf(pulses, plot = FALSE)

histogram <-  hist(diff(pulseTimes), breaks=xc_pulses$lag, plot=FALSE)

plot(xc_times$lag,xc_times$acf, type="l", ylim=c(-0.1,1))
lines(xc_pulses$lag, xc_pulses$acf, col="red")
lines(histogram$mids, histogram$counts / max(histogram$counts), col="green")

legend(12, 1.0, c("ACF of times", "ACF of pulses", "Histogram of inter-pulse times"), lwd=c(2.5,2.5, 2.5),col=c("black","red", "green"))


MFCC - Significance of number of features

Number of filter banks

One of the last steps in the MFCC's calculation is measuring the energy in the filter banks. We do that because want to reduce the dimensionality of our input vector (amplitude spectrum), as well as capture its envelope. Those triangular filters are spaced over the Mel scale:

This means that we have very good resolution in low frequencies. Exactly opposite is true for higher frequencies. We do that because MFCC's are suited for speech-related tasks and most of the information is located in lower frequencies (i.e. formants).

So how many filter banks do we actually want? For example HTK by default is using 20 filter banks. You might increase that number, especially if you are dealing with signals that contain a lot of closely spaced frequencies and you want to resolve them between each other. It's totally up to you - in the end what really matters is the classification performance.

Number of coefficients

After taking the logarithm of energies in each filter bank, the last step (or second last if you are doing some liftering) is to calculate the Mel Frequency Cepstral Coefficients. We do that by fitting the cosines to calculated log energies using the DCT. This captures the periodicity in the reduced spectrum. Figure below should help you in understanding that process. You can imagine that for number of coefficients equal to number of filter banks this would correspond to capturing the alternating energy between each of the filter banks.

So how many MFCC's do we want to calculate? HTK by default uses 12 and in most applications it is more than enough. In general we don't want too many coefficients because:

It's all about reducing the dimensionality of our feature space.
One DCT properties is that it de-correlates and keeps most of the information in first few coefficients.

Again it all depends on your application and you should adjust this number based on the recognition performance. From my experience the gain from increasing the number of MFCC's was negligible compared to introducing the $\Delta$ and $\Delta\Delta$ coefficients.


Essentially, isn't this just looking for local maxima in the image? Have a look at the following code. It dilates the image, i.e. putting each pixel the value of the maximum of its neighborhood. Then, it checks which pixels actaully were the maximum and marks them with a circle:

import skimage
import skimage.io
import skimage.morphology
import numpy as np
import skimage.draw
import scipy.ndimage

I = skimage.io.imread("Capture.PNG", as_grey=True)

selem = np.ones((9,9))
D = skimage.morphology.dilation(I, selem)
plt.figure(figsize=(12,12))

plt.subplot(221); plt.imshow(I); plt.title("Original image")
plt.subplot(222); plt.imshow(D); plt.title("dilated image")

maxima = I==D
plt.subplot(223); plt.imshow(maxima); plt.title("Detected Maxima regions")

labels, nlabels = scipy.ndimage.label(maxima)

plt.subplot(224); plt.imshow(labels)
print (nlabels)

C = scipy.ndimage.measurements.center_of_mass(maxima, labels, 1+np.arange(nlabels-1))
I2 = I.copy()
for x, y in C:
    rr, cc = skimage.draw.circle_perimeter(int(x), int(y), 5, shape=I2.shape)
    I2[rr,cc] = 1
plt.imshow(I2); plt.title("Center of mass of maxima"); plt.colorbar();

what is missing is e.g. a threshold to ignore small maxima.


How to deal with a negative pole (unstable) in the pre-filter of a control system?

In general, in order to stabilise a system, however complex, if you have the transfer function $G(s)$, you introduce a feedback loop with a new function $F(s)$.

Write the closed loop transfer function for the new system with the added $F(s)$, and then find $F(s)$ in order for the new system to be stable. This is like the first exercise in any control book in order to offer an example of stabilising a system via negative Feedback.

Check the book from Ogata on control engineering for reference.


Discrete time non-linear time invariant system dynamics descriptions (state-space or input-output relationship)

Well, any input-output representation obviously admits a state-sapce form. for your equation in $y[k]$ you can easily construct one as follows. Create a "shift" system (an integrator chain) as
$$
\begin{aligned}
x_1[k+1] &= x_2[k],\\
x_2[k+1] &= x_3[k],\\
&\vdots\\
x_n[k+1] &= y[k]
\end{aligned}
$$
In this way indeed you have $x_n[k] = y[k-1]$, $x_{n-1}[k]=y[k-2]$, ..., $x_1[k] = y[k-n]$. You can do the same by defining another shift system to carry the information in $u[k]$, namely
$$
\begin{aligned}
z_1[k+1] &= z_2[k],\\
z_2[k+1] &= z_3[k],\\
&\vdots\\
z_m[k+1] &= u[k]
\end{aligned}
$$
Thus, with $\xi := col(x_1,...,n_n,z_1,...,z_m)$ the IO representation is equivalent to
$$
\begin{aligned}
\xi[k+1] &= A\xi[k],\\
y[k] &= h(\xi_1[k],\dots,\xi_n[k],u[k],\xi_1[k],\dots,\xi_m[k])
\end{aligned}
$$
The converse is instead not trivial. Pick the LTI system
$$
\begin{aligned}
x[k+1] &= \begin{bmatrix}2 & -1\\0 & 0.1\end{bmatrix}x[k] + \begin{bmatrix}1\\1\end{bmatrix}u[k],\qquad (1)\\
y[k] &= x_2[k]
\end{aligned}
$$
Then the output is the same as
$$
y[k] = 0.1 y[k-1] + u[k],\qquad (2)
$$
however, (2) and (1) are not equivalent.
In general, given a IO representation you can always find a state-space representation, however for general state-space forms a IO representation might be not enough to describe the whole state information. You should find a diffeomorphism among the two representations to assess equivalence


Square wave FFT results show frequency that seems too low?

The frequency resolution of your FFT is determined by its length and the sampling rate.

In your second plot you show Fs = 10e6 and 3200 samples. Assuming you don't zero-pad these samples when you pass them to the FFT this gives you a frequency resolution of 3125Hz:

$$ 3215 = \frac{10 \times 10^6}{3200} $$

So at lower frequencies you have frequency bins as 0Hz, 3125Hz, 6250Hz and 9375Hz. The largest peak in the FFT will be in the bin closest to the signal fundamental frequency.

If you want better resolution you either need to capture more samples, or append a load of zeros to the end of your time series (zero-padding).

I would recommend reading about windowing for FFT analysis also, especially if you go down the zero-padding route.

As an aside: the reduce_sample_array function seems to decimate the time series by simply discarding samples. You really should be doing proper down-sampling here, to avoid aliasing.


Impulse response of a 3x3 PSF - how to find analytical expression for fourier transform of a 3x3 matrix?

You can use the 2D DFT formula

$$H(\omega_1,\omega_2) = \sum_{n_1} \sum_{n_2} h[n_1,n_2] e^{-j(\omega_1 n_1 + \omega_2 n_2)} $$

and simply the trigonometric algebra to get a closed form analytic expression for the 2D-DTFT. However, as @LaurentDuval has already mentioned, your 3x3 kernel is separable and one set of 1D filters is this

$$f[n_1]  = [\frac{1}{4} , \frac{1}{2}, \frac{1}{4}]^T $$
$$g[n_2]  = [\frac{1}{2} , 1 , \frac{1}{2}] $$

Then from the Fourier transform we know that

$$ h[n_1,n_2] = f[n_1]g[n_2] \implies H(\omega_1,\omega_2) = F(\omega_1) G(\omega_2) $$

Analytic expressions for 1D-DTFT's can be obtained from
$$H(\omega) = \sum_{n}h[n] e^{-j \omega n} $$

Assuming 3x3 kernel has zero phase, then we see that :
$$F(\omega_1) = \frac{1}{4} e^{j \omega_1 } + \frac{1}{2} e^{j 0 } + \frac{1}{4} e^{-j \omega_1 } = 0.5 \cos(\omega_1) + 0.5 $$

similarly
$$G(\omega_2) = \frac{1}{2} e^{j \omega_2 } + 1 e^{j 0 } + \frac{1}{2} e^{-j \omega_2 } = \cos(\omega_2) + 1 $$

Hence:
$$H(\omega_1,\omega_2) = 0.5 \left( \cos(\omega_1) + 1 \right)(\cos(\omega_2) + 1) $$


Alamouti Scheme. MISO. STBC for real constellation

Your analysis is correct. I would write it as follows: if you define the transmitted matrix as
\begin{equation*}
X = \begin{bmatrix} x_1 & -x_2 \\ x_2 & x_1
\end{bmatrix}
\end{equation*}
with one column transmitted at each symbol interval (that is, first you transmit $[x_1, x_2]$ and then $[-x_2, x_1]$), then the system can be written as $$Y=hX$$ with $h = [h_1, h_2]$, or equivalently it can be written as $$Y = Hx$$ with
\begin{equation*}
H = \begin{bmatrix} h_1 & h_2 \\ -h_2 & h_1
\end{bmatrix}
\end{equation*}
and $x=[x_1, x_2]^T$.

Of course, you can also transpose all vectors and matrices and arrive at equivalent results.


Time invariance and linearity of recursive system

The system is both linear and time-invariant.

Theory says linearity means that for an input $$ax_1[n] + bx_2[n]$$ your output would be $$ay_1[n] + by_2[n]$$ where $y_1[n], y_2[n]$ are the outputs of $x_1[n], x_2[n]$, respectively.

In your case, an input $x_1[n]$ gives output $$y_1[n] = y_1[n-1] - y_1[n-2] - x_1[n]$$ same for an input $x_2[n]$, thus $$y_2[n] = y_2[n-1] - y_2[n-2] - x_2[n]$$

Multiplying the first equation by $a$ and the second by $b$ we get $$ay_1[n] = ay_1[n-1] - ay_1[n-2] - ax_1[n]$$ and $$by_2[n] = by_2[n-1] - by_2[n-2] - bx_2[n]$$

Adding together, we get
$$\begin{align}
ay_1[n] + by_2[n] &= ay_1[n-1] - ay_1[n-2] - ax_1[n] + by_2[n-1] - by_2[n-2] - bx_2[n] \\
ay_1[n] + by_2[n] &= ay_1[n-1] + by_2[n-1] - ay_1[n-2] - by_2[n-2] - ax_1[n] - bx_2[n])
\end{align}$$
which is the given system with $x[n] = ax_1[n] + bx_2[n]$ and $y[n] = ay_1[n] + by_2[n]$. So it's linear.

For time-invariance, we first delay the input by $k$, thus we get $$y[n] = y[n-1] - y[n-2] - x[n-k]$$

Can you show that for a delayed output (with same delay $k$), the same result is obtained?


Using SURF algorithm to match objects on MATLAB

People generally use Euclidean distances, just like you do and it works as it is, pretty well. However, if you require further robustness, I propose the following:

Compute the matching both ways: image1->image2 and image2->image and then retrieve only the agreeing matches. This introduces some consensus.
Assuming you have planar patches, you run a homography using RANSAC. The good thing here is that RANSAC is robust to outliers, given a predifined tolerance. This way prune away a lot of matches, and you will be left with a meaningful set of inliers.

Using these two together should really give you robust and repeatable matches.


How to implement bandpass filter on complex valued signal?

You can do different things: For example,

use the (frequency-symmetric) real-tapped bandpass that firwin gives you, and after applying that, apply a complex high pass (a Hilbert filter, essentially) to kill all negative frequencies. That sadly leaves you with the original problem (finding a complex-tapped filter using scipy)
Do the same as above, but already convolve the real-valued band pass filter taps with the Hilbert filter taps, so that you only need to apply one filter
Do the "usual" lowpass-to-bandpass transform trick: You design a real-valued low pass filter (e.g. using firwin) with impulse response (==taps) $h_{LP}[n]$ that has the same passband and transition widths as your desired complex band pass. Afterwards, you multiply that $h_{LP}$ with the shift you need to move the 0-centric low pass up to your $f_\text{center}$-centric band pass: $$h_{BP}[n] = h_{LP} \cdot e^{j2\pi \frac{f_\text{center}}{f_\text{sample}}n}\text.$$ The idea is that you shift the low pass in frequency by convolving with a dirac impulse in frequency domain, which is equivalent to multiplying with the complex sinusoid of said frequency in time domain. These answers (1,2) might be of interest to you.


Can we start a wave $\sin(2x-4)$ from position 4 in the horizontal axis?

You can shift the origin to $x=a$ by replacing the argument $x$ by $x-a$:

$$\sin(2(x-a))\tag{1}$$

Since for $\sin(2x-4)=\sin(2(x-2))$ you have $a=2$, the origin is shifted to $x=2$. If you want to shift it to $x=4$, you get from Eq. $(1)$

$$\sin(2(x-4))=\sin(2x-8)\tag{2}$$


Intertial navigation on android phone with Kalman filter

The quick answer is that its not really possible to get the position information from a phone IMU with gyroscope, compass and accelerometer. Even with a Kalman filter its too unreliable to get any linear acceleration signal that can be integrated to give position. I have looked at a lot of Kalman filter papers for IMUs and at best, the linear acceleration is treated as a nuisance term. The three sensor IMU with Kalman filter does however give pretty solid 3D rotation estimates.


We must produce s + 3 images in the stack of blurred images for each
  octave, so that ﬁnal extrema detection covers a complete octave.

For $s=3$ this means you will have $s + 3 = 6$ blurred images (the Gaussian images shown in the paper in Figure 1 on the left). Having $6$ Gaussian images will result in $5$ DoG images (shown in Figure 1 on the right). This will allow you to do the extrema detection on $s=3$ scales (using the method shown in Figure 2).


Detecting a frequency pattern in the frequency domain

Try

xcorr(trainingfft,datafft,'option')

to specified a normalization for the cross-correlation. See the file page for more details of the option.


finding energy content of signal from "energy spectral density" function

The so called auto-correlation function that's used in the example (which was found in a previous example) is not correct.

First note that it's the deterministic auto-correlation computation, compared to its theoretical counterpart which requires that the signals be random for their ACF to be defined at all.

Anyway, so the following integral signifies one definition of deterministic ACF of a given deterministic signal $x(t)$ as :

$$R_{xx}(\tau)  = \int_{-\infty}^{\infty} x(t)x(t-\tau) dt$$

applying this definition to the signal $x(t) = e^{-a t} u(t)$ with $a >0$ we have:

$$
\begin{align}
R_{xx}(\tau)  &= \int_{-\infty}^{\infty} x(t)x(t-\tau) dt \\ \\
&= \int_{-\infty}^{\infty} e^{-a t} u(t) e^{-a (t-\tau)} u(t-\tau) dt \\ \\
&= e^{a \tau} \int_{-\infty}^{\infty} e^{-2 a t} u(t) u(t-\tau) dt \\ \\
&= e^{a \tau} \int_{\max\{0,\tau\}}^{\infty} e^{-2 a t} dt \\ \\
&= \begin{cases}{ e^{a\tau} \int_{\tau}^{\infty} e^{-2at} dt ~~~,~~~\tau \geq 0 \\ \\ \\
e^{a\tau} \int_{0}^{\infty} e^{-2at} dt ~~~,~~~\tau < 0 }\end{cases}\\ \\
&= \begin{cases}{ e^{a\tau} \frac{e^{-2a\tau}}{2a} ~~~,~~~\tau \geq 0 \\ \\ \\
e^{a\tau} \frac{1}{2a} dt ~~~,~~~\tau < 0 }\end{cases}\\ \\
&= \begin{cases}{ \frac{e^{-a\tau}}{2a} ~~~,~~~\tau \geq 0 \\ \\ \\
\frac{e^{a\tau}}{2a} dt ~~~,~~~\tau < 0 }\end{cases}\\ \\
R_{xx}(\tau) &=   \frac{e^{-a |\tau|}}{2a}  \\ \\
\end{align}
$$


Orthogonality provides an interesting backbone to the structure of the filter-banks (FB). First, from an analysis FB, the synthesis FB is very direct, so it can ease implementations. Second, the orthogonality often allows faster implementations, as there is "little redundancy" in computation. Third, orthogonality ensures that matrices are well-conditioned, this may reduce risks of error propagation (esp. with quantization). Fourth, energy is preserved, so operations performed in the dual domain can be measured (energetically) in the primal domain. Fifth: it's a great deal of simplification when you need to get tractable proofs for their design (energy compaction for instance) or performance (as for instance Gaussian noise remains Gaussian after an orthogonal transformation). Note that these five points are not completely orthogonal.

To provide a little more details:

Orthogonality, and more generally paraunitarity (cf. Paraunitary matrices by Barry and Ted Hurley), i.e. a square matrix $U(z)$ satisfying $U(z)U^∗(z^{−1}) = 1$, are practical, as the inverse is somehow a conjugate-transpose
For orthogonal bases in an $N$-dimensional space, the first orthogonal vector has about  $N-1$ degrees of freedom, the second one $N-2$, etc. So orthogonal matrices have a structure (see the Orthogonal group), which can be used to factorize them, hopefully with less operations (Fast multiplication of orthogonal matrices). For instance, $f 3 \times 3$ orthogonal matrices can be represented by  three Givens rotations, but other designs are possible
Determinants are $\pm 1$ (and of course the same of their inverses). When a matrix contains close vectors, projections split evenly of them, which has effects on 1) simple inversion 2) iterations (like in computerized tomography) where powers of matrices (like FFT)  over millions of cells/pixels/voxels may diverge.
When energy is preserved, a structured signal (in the appropriate basis or frame) may come out of unstructured noise. This is a rationale behind "sparsifying transformations", and how thresholding values can denoise. Somehow behind Stein estimates (Stein Unbiased Risk Estimation), with orthogonaly invariant estimators, sample eigenvalues are modified (for the best) yet sample eigenvectors untouched.
All block orthogonal transforms are quite easy to deal with.

a) Overlap and windows are more involved. An early attempt is the LOT (Lapped Orthogonal Transform), and their generalizations (GenLOT - Generalized Lapped Orthogonal Transforms , GULLOT - Generalized Unequal Length Lapped Orthogonal Transforms , etc.). They can be seen as DCT (or Hadamard) matrices with additional degree-of-freedom. By decomposing them with rotations, it is possible to optimize their design by imposing constrains, like "Coding gain", "DC leakage" (Linear phase paraunitary filter bank with filters of different lengths and its application in image compression, 1999), that can be related to their performance in data compression. Trivia: read Properties of the Eigenvectors of Persymmetric Matrices with Applications to Communication Theory to see how correlation structures impose properties of eigenvectors
b) Some noise properties can be estimated in a data sparsifying domain. A classical example in the wavelet domain (and filter-banks generally) is the MAD (median absolute deviation) estimator of the Gaussian noise in a combination of subbands.

However, orthogonality has limits, imposes additional constraints that sometimes spoil other desirable properties. Plus, the human sensory systems care little about  strict orthogonality, and you can gain degrees of freedom (even in computational efficiency) by discarding strict orthogonality. The filterbanks used in the JPEG2000 image compression standard are not orthogonal, but biorthogonal. [UPDATE] In some cases, you may want integer or dyadic rational operations (for memory or speedup), and this is rarely compatible with orthogonality. In video coding, a lot of decorrelating operations are not strictly orthogonal.

On most practical cases though, a lot of systems remain at least "close to orthogonal" in some sense.


Memory of a filter used on real time digital signal processing

If memory is your primary concern then an IIR filter is the way to go.  It uses very little memory and is pretty computationally efficient.

FIR filters' memory length is one less than the number of filter taps.  This usually isn't that much, but it can add up if there are a number of filters and/or they are long.


LoG filter creating additional maxima in scale space

Without exactly seeing your code it is hard to know for certain, but I suspect this is because you are not tracking your maxima well through scale space.

The position of the maxima will move in scale space. At the smallest scales you get the 4 small blobs around each circle. These are due to digitization and possibly noise/irregularities in the boundary edge.

These maxima will expand as the scale increases and merge together to form the "correct" blobs around the circles.

However, the other side of the peaks will also merge together resulting in the 12 medium size blobs and finally these will merge to form the four largest blobs.

Note at each stage here the number of maxima decreases 36->21 (9+12) -> 13 (9+4)

The problem is as you only look for local maxima in scale-space you lose  a lot of the relationship between the maxima.

One option may be to look for all local maxima at each scale and track these through scales. However, I'm not sure how straightforward this would be or how you determine which are the significant peaks using this approach.

A possibly better option may be similar to what you suggested by thresholding the responses. This article does basically the same thing as you. It applies a threshold on the mean blob volume over all scales. Which it claims it better than thresholding the peak response.


It is possible that two decoders will have different results and be compliant with standard. The standard "ISO/IEC 11172-4 (part 4) compliance testing" specify two classes of decoders:
"accurate" and with "limited accuracy" *).

The accuracy is measured by comparing decoded stream with reference stream (reference input and output is specified in standard). If the RMS of difference is below $2^{-15}/\sqrt{12}$ the decoder is classified as accurate and  if below $2^{-11}/\sqrt{12}$ it is classified as "limited accuracy" one.

So if you are going to use different decoders or even the same decoder on different processor architectures you have to be prepared for this kind of discrepancy in output. Do not know your application or intend but maybe you can mimic the standard and work with RMS of error too?

*) I do not have text of the standard in English in front of me so actual name used in standard may be different.


My signal is a time-signal (displacement vs time), but the filter takes the input in terms of cut-off frequency. Does this mean I have to convert the signal into frequency domain, or does the filter take care of it? I am bit confused on how this actually works.

You design or specify the filter in the frequency domain, but the actual filtering can be directly accomplished in time.  There are FFT based filtering algorithms that are typically more efficient, but there is no requirement to transform your signal to the frequency domain.

On a related note, I don't know the frequency of each motion that I am trying to filter out i.e. the intuition is that the frequency for different motions be different, but how do I determine by cut-off values? Is it usually done by trial and error, or is there a more scientific way to accomplish this?

Anything that doesn't require trial and error is suspect, scientific methodology wise. I suggest that you try to be more specific about what filtering out different motions actually means.  Some kinds of phenomena are a superposition of features which are often separable, of which some can be "filtered out" and others don't lend themselves to that.

You might want to use a tool known as the Short Time Fourier Transform (STFT) and see if you can eye ball any patterns associated with desired/undesired motion.

Finally, try Google Scholar.


Determining Stability of a continuous time system using Laplace Transform

ROC is a strip of plane which does not contain any poles. The region is determined by the location of poles.

The first signal, $x(t) = e^{-t} u(t)$ , has the Laplace transform $X(s) = \frac{1}{s+1}$ and its ROC is $ \mathcal{Re}\{s\} > -1 $ , since $x(t)$ was defined to be causal.

The second signal, $y(t)$, is obtained by advancing $x(t)$  by $1$ unit in time. Hence $y(t) = x(t+1)$. By Laplace transform property it's seen that $Y(s) = \frac{e^{s}}{s+1}$ , and it has the same ROC with $X(s)$, since they are both right sided sequences. However, $y(t)$ is not causal despite being right sided, due to the advance of $1$ unit to the left.


Differences between soft knee and hard knee in Dynamic Range Compression (DRC)?

Regarding advantages - there is only one I can think of... The soft knee makes the compressor less noticeable. It sounds more natural, especially when you are dealing with higher compression ratios (like 4:1, etc.).

On the other hand, if you want to participate in the "Loudness War", then you should rather use hard knee, as it allows you to squeeze out even more loudness.

Regarding implementation. Assuming the:

$T$ - threshold

$W$ - soft knee width (soft knee extends from $T-\frac{W}{2}$ up to $T+\frac{W}{2}$)

$R$ - compression ratio

The soft knee compressor is governed by the equation:

$$
y_G = \begin{cases}
x_G & x_G-T<\frac{-W}{2} \\
x_G+\frac{\left(\frac{1}{R}-1\right)\left( x_G - T+ \frac{W}{2}\right)^2}{2W} & |x_G-T| \le \frac{W}{2} \\
T+\frac{x_G-T}{R} & x_G - T > \frac{W}{2}
\end{cases}
$$

Here is how the compression curve will look like:

If you want to achieve the hard knee, then just drop the middle equation:

$$
y_G = \begin{cases}
x_G & x_G < T \\
T+\frac{x_G-T}{R} & x_G > T
\end{cases}
$$

Obviously then you will get the following compression characteristic:

Whenever you have to do something about Dynamic Range Compression, I suggest you to take a look at the article below. Those guys did a great job of putting together all the info. Nonetheless I highly recommend reading the articles in references.

  Giannoulis D. et al. - Digital Dynamic Range Compressor Design - A
  Tutorial and
  Analysis


Transfer function of an Exponential system in Z domain

This system doesn't have a transfer function, so there's no way to find it. Only linear time-invariant (LTI) systems can be described by a transfer function. In discrete-time, the $\mathcal{Z}$-transforms of input and output sequences of an LTI system are related by

$$Y(z)=H(z)X(z)\tag{1}$$

where $H(z)$ is the system's transfer function, which is the $\mathcal{Z}$-transform of its impulse response. Equation $(1)$ is another way to state that the output sequence is given by the convolution of the input sequence with the system's impulse response. If the system is not LTI there is no function $H(z)$ which relates $X(z)$ and $Y(z)$ as in Eq. $(1)$.

Your system is definitely non-linear so whoever asked you to find its transfer function is either clueless or wanted to fool you.


Positive and negative frequencies in DFT due to frequency folding, or due to negatively indexed frequencies?

I think that considering the DFT from a linear algebraic point of view has some value, so I will attempt to introduce the foundations.

We will assume that our signal is a vector of $N$ complex entries.

$\mathbb{C}^N$ is the vector space of vectors with $N$ complex entries. Let $\mathbf{u}_0,\mathbf{u}_1,\ldots,\mathbf{u}_{N-1}$ be vectors in $\mathbb{C}^{N}$ defined by
\begin{equation}
\begin{split}
\mathbf{u}_k
&=~
\frac{1}{\sqrt{N}}\left(\begin{array}{c}
\exp(2\pi \mathsf{j}\times 0\times(k/N))\\
\exp(2\pi \mathsf{j}\times 1\times(k/N))\\
\exp(2\pi \mathsf{j}\times 2\times(k/N))\\
\exp(2\pi \mathsf{j}\times 3\times(k/N))\\
\vdots\\
\exp(2\pi \mathsf{j}\times (N-1)\times(k/N))
\end{array}\right)\\
&=~
\frac{1}{\sqrt{N}}\left(\begin{array}{c}
1\\
e^{2\pi\mathsf{j}k/N}\\
(e^{2\pi\mathsf{j}k/N})^2\\
\vdots\\
(e^{2\pi\mathsf{j}k/N})^{N-1}
\end{array}\right),
\end{split}
\end{equation}
for $k = 0,1,2,\ldots,N-1$, where $\mathsf{j} = \sqrt{-1}$.

Every entry of $\mathbf{u}_0$ is $1/\sqrt{N}$, so $\mathbf{u}_0$ might be considered as a sampled DC signal.
The entries of $\mathbf{u}_1$ are samples of a complex exponential with fequency $\frac{1}{N}$,
The entries of $\mathbf{u}_2$ are samples of a complex exponential with fequency $\frac{2}{N}$,
and so on, up through frequency $\frac{N-1}{N}$.

$\mathbf{u}_0,\mathbf{u}_1,\ldots,\mathbf{u}_{N-1}$ form an orthonormal basis for $\mathbb{C}^{N}$, which means that each $\mathbf{u}_k$ has norm 1, they are all orthogonal to one another, and each vector in $\mathbb{C}^{N}$ can be represented unambiguously as a linear combination of them. An important result of this is that, if $\mathbf{x}\in\mathbb{C}^{N}$, then there is exactly one list of complex numbers $c_0,c_1,\ldots,c_N$ such that
\begin{equation}
\mathbb{x} = c_0\mathbf{u}_0 + c_1\mathbf{u}_1 + \cdots + c_{N-1}\mathbf{u}_{N-1}.
\end{equation}

The coefficients mentioned above are the entries of the DFT of $\mathbf{x}$:
\begin{equation}
\mathbf{x} = X[0]\mathbf{u}_0 + X[1]\mathbf{u}_1 + \cdots + X[N-1]\mathbf{u}_{N-1}.
\end{equation}
We might interpret $X[0]$ as the strength of the DC component of $\mathbf{x}$, $X[1]$ as the strength of the component of $\mathbf{x}$ with frequency $\frac{1}{N}$, and so on. Since $\mathbf{X} = \mathsf{DFT}\mathbf{x}$ has complex entries, there is some phase information attached to each "strength".

So far, we have considered only components of non-negative frequencies. What if we would rather view $\mathbf{x}$ as a combination of negative and positive frequency components? Consider a component of frequency $-\frac{k}{N}$ for $0< k \leq \frac{N}{2}$:
\begin{equation}
\mathbf{u}_{-k}
~=~
\frac{1}{\sqrt{N}}\left(\begin{array}{c}
\exp(2\pi \mathsf{j}\times 0\times(-k/N))\\
\exp(2\pi \mathsf{j}\times 1\times(-k/N))\\
\exp(2\pi \mathsf{j}\times 2\times(-k/N))\\
\exp(2\pi \mathsf{j}\times 3\times(-k/N))\\
\vdots\\
\exp(2\pi \mathsf{j}\times (N-1)\times(-k/N))
\end{array}\right).
\end{equation}
The $\ell^{\textrm{th}}$ entry of this vector is
\begin{equation}
\begin{split}
u_{-k}[\ell] &=~ \frac{1}{\sqrt{N}}\exp\left(2\pi\mathsf{j}\times\ell\times\frac{-k}{N}\right)\\
&=~ \frac{1}{\sqrt{N}}\exp\left(2\pi\mathsf{j}\times\ell\times\frac{-k}{N}\right)\times\underbrace{\exp\left(2\pi\mathsf{j}\times\ell\times\frac{N}{N}\right)}_{1}\\
&=~ \frac{1}{\sqrt{N}}\exp\left(2\pi\mathsf{j}\times\ell\times\frac{N-k}{N}\right)\\
&=~ u_{N-k}[\ell].
\end{split}
\end{equation}
In other words, the negative-frequency component $\mathbf{u}_{-k}$ is exactly the same as the positive-frequency component $\mathbf{u}_{N-k}$.

Suppose that $N = 2M$ for some positive integer $M$. Then
\begin{equation}
\begin{split}
\mathbf{x} &=~ X[0]\mathbf{u}_0 + \cdots + X[N/2-1]\mathbf{u}_{N/2-1} + X[N/2]\mathbf{u}_{N/2} + \cdots + X[N-1]\mathbf{u}_{N-1}\\
&=~ X[0]\mathbf{u}_0 + \cdots + X[M-1]\mathbf{u}_{M-1} + X[M]\mathbf{u}_{M} + \cdots + X[N-1]\mathbf{u}_{N-1}\\
&=~ \underbrace{X[0]\mathbf{u}_0 + \cdots + X[M-1]\mathbf{u}_{M-1}}_{\textrm{non-negative-frequency components}}
+ \underbrace{X[M]\mathbf{u}_{N-M} + \cdots + X[N-1]\mathbf{u}_{N-1}}_{\textrm{negative-frequency components}}
\end{split}
\end{equation}
For a full decomposition, one can choose the frequency-sets
\begin{equation}
-\frac{N/2}{N},-\frac{N/2-1}{N},\ldots,-\frac{1}{N},0,\frac{1}{N},\ldots,\frac{N/2-1}{N}
\end{equation}
or
\begin{equation}
0,\frac{1}{N},\ldots,\frac{N-1}{N},
\end{equation}
each of which consists of $N$ distinct frequencies. In truth, one can choose other frequency-sets of $N$ frequencies, too, but these are the ones to which we have attached some intuition over the decades.

MATLAB's fft gives the DFT with all non-negative frequencies. To convert the output of fft to the vector of coefficients for negative, zero, and positive frequencies, one applies the fftshift function.

All of this and much more is explained from a linear algebraic point of view in Linear Algebra, Signal Processing, and Wavelets - A Unified Approach (PDF) by Øyvind Ryan of the University of Oslo.


How to compute Laplace Transform in Python?

Update : Given the approach started in the OP's Github code I have this suggestion :
Observe that the unilateral Laplace Transform given as:
$$X(s) = \int_0^\infty x(t)e^{-st}dt$$
Is just the Fourier Transform of a causal function with a weighting exponential:
$$X(s) = \int_0^\infty x(t)e^{-(\sigma+j\omega)t}dt$$
$$X(s) = \int_0^\infty [e^{-\sigma t}x(t)]e^{-j\omega t}dt$$
Which is the Fourier Transform of $e^{-\sigma t}x(t)$.
So to proceed, learn how to produce surface plots in python, and then index through $\sigma$ within the Region of Convergence (see below) and compute the FFT of $e^{-\sigma t}x(t)$  to create your complex surface values given $\sigma$ and $\omega$.

The challenge with what you are trying to do is that the Laplace Transform is a function of the complex variable "s", so for each possible value of "s" (which is simply the set of all complex numbers) the Laplace Transform would have a complex result with a magnitude and phase.
So when you try to plot this, you could plot the magnitude and phase as two separate plots but these would be 3D surface plots showing either the result of the magnitude or the phase over a 2D surface, similar otherwise to how we plot the magnitude and phase of the Fourier Transform as a 2D function over the line representing all frequencies.
For the graphical representation of the Laplace Transform, we typically just show the locations where that function goes to infinity (poles) or is zero (zeros). In fact every other sample on the surface is uniquely defined by the pole and zero locations, so that is all we need to show to define it.
This is best shown with an example. Since I already have the graphic for this particular case, consider the time domain function of a decaying sinusoid given by the formula below and the plot below that where we see in the dashed red line the envelope for the decaying function $2e^{-0.2t}$. $u(t)$ is the step function which is $0$ for time $t<0$ and $1$ for time $t \ge 0$.
$$ x(t) = u(t)2e^{-0.2t}sin(0.5t)$$

To get the Laplace Transform (easily), we decompose the function above into exponential form and then use the fundamental transform for an exponential given as :
$$\mathscr{L}\{u(t) e^{-\alpha t}\} = \frac{1}{s+\alpha}$$
This is the unilateral Laplace Transform (defined for $t = 0$ to $\infty$), and this relationship goes a long way since we can describe the response of any causal linear system using such exponential forms.
So the equation above, assuming $t>0$, and using Euler's identity becomes :
$$ x(t) = u(t)2e^{-0.2t}sin(0.5t)$$
$$ = 2e^{-0.2t}\frac{(e^{+j0.5t}-e^{-j0.5t})}{2j}$$
$$ = je^{-0.2t}(e^{+j0.5t}-e^{-j0.5t})$$
$$ = je^{-(0.2-j0.5)t} -je^{-(0.2+j0.5)t}$$
Which we can then easily take the Laplace Transform to get a function of the complex variable $s$ :
$$ X(s) = \frac{j}{s+0.2+j0.5}-\frac{j}{s+0.2-j0.5} = \frac{1}{(s+0.2+j0.5)(s+0.2-j0.5)}$$
The graph of the magnitude of this is shown below, which is a surface plot for all values of $s$, where $s$ is a complex variable with real and imaginary components traditionally described as $s = \sigma +j\omega$.  So we simply choose a particular complex value $s$, and plot the magnitude of the result of $X(s)$.

Caveat : This is indeed the magnitude of the Laplace Transform if the Laplace Transform exists. For causal functions, the Laplace Transform does not exist for all values of s whose real part is to the left of the leftmost pole. Why? Because for any values of $s$ in that region (left of the leftmost pole), the Laplace Transform (given by the integral) will not converge (grows to infinity). While for all values of s at and to the right of the leftmost pole, the transform will converge to be the magnitude on the surface plotted above. So the proper graphic would only be valid to the right of the right most pole. This can be confusing as $X(s)$ certainly converges as was done to make this plot, but it is in the transformation itself that the result cannot be obtained.
This is very clear if we consider the envelope in our example function $x(t)$ which was given by $2e^{-0.2t}$ for all $t>0$, and the Laplace Transform $X(s)$ for $s = -1$ :
$$X(s) = \int_0^\infty x(t)e^{-st}dt$$
$$X(s= -1) = \int_0^\infty 2e^{-0.2t}e^{t}dt = \int_0^\infty 2e^{0.8t}dt$$
which will not converge since the function $e^{0.8t}$ continuously grows larger for larger $t$.
Hopefully after reading this you no longer feel the need to plot the Laplace Transform, and in practical application a plot of it is never used beyond showing the pole and zero locations. However what is very useful is knowing that the Fourier Transform is the Laplace Transform when $s = j\omega$. You can see this if you compare the two equations, and the small breakout in upper right-hand corner of the plot above is also showing this, which is the Frequency Response specifically. The one precaution is that the Fourier Transform is often given as a bilateral function (t extending from $-\infty$ to $\infty$) so to be truly equivalent unless the function is declared to be causal, we must be using the bilateral Laplace Transform for the two to be exactly identical (which is also seldom used). This explains why the Fourier Transform of a sine wave would appear as two impulses at $\pm \omega_c$ while as you might deduce from the plot above, the Laplace Transform of the sine wave would appear more as tent poles on the $j\omega$ axis. You would see this as well with the Fourier Transform if you did the transform of $u(t)cos(\omega t)$ specifically.


What interpolation methods can I use to get the tightest fit for these curves?

If you have a time series of images, I recommend the CONDENSATION algorithm for the interpolation of curves (see also the wikipedia page). There are a number of implementation available, notably in opencv.


Error introduced by smoothing filter at the end of a processing block

A moving average filter is a FIR filter of equal weights 1/N and a FIR filter is a degeneralization of a IIR filter.
In a non realtime implementation, it is functionally equivalent if you consider the filter to be N samples delay, N samples lookahead or N1 samples delay and N2 = N-N1 samples look ahead. Pick whatever time philosophy makes sense in your head and minimize tedious indexing.
Matlabs filter() and filtic() functions provide the means for carrying the  result from block to block for IIR filtering. Sadly only single channel support.
If you use something like:
conv2(ones(5,1), 1, X, ‘full’);
You will get an output vector that is larger than the input vector. Those surplus samples are exactly what you need to add to the initial samples of the next block.


How to design a lowpass filter using Hamming window and sinc function

If you can’t make the lengths of your Sinc kernel and window function the same, you can zero-pad the shorter of the two before you multiply them together.


How to removed ripples cause by zero-padding in fft?

Windowing the time domain function prior to zero padding will reduce the ripples caused by zero padding (due to the otherwise sharp truncation of the time domain function). Without additional windowing you are equivalently multiplying your longer time domain function with a rectangular window, which Fourier Transform is a Sinc function. To reproduce this exactly in frequency the frequency must extend to infinity, which is doesn't in the discrete system so aliasing of these tails in frequency occurs, resulting in ripple. Choose a window whose Fourier Transform has tails that are much lower at the aliasing frequencies and the ripple will be significantly reduced! (I personally like the Kaiser window as you can adjust the trade between width of the main lobe (frequency resolution) and depth of the tails (dynamic range) by adjusting the Beta factor in the window construction.


the FFT is a fast method of doing the DFT.    yes, multiplication in one domain corresponds to convolution in the reciprocal domain.  and since the DFT is a totally circular operation, any consequential convolution done in one domain by the DFT by multiplication in the other domain is circular convolution.  you don't have a choice about it.  to make a circular convolution look like linear convolution, you must zero-pad sufficiently.

how many zeros to pad?  need to consider how fast convolution (overlap-add or overlap-save) works.  if $N$ is the DFT length, the non-zero lengths of the two sequences being convolved, must not add to more than $N+1$ otherwise there will be overlapping and aliasing of non-zero components.


How to find the perfect regularization number?

Approaches can be any of the following:

Model the noise that you expect, for ex: gaussian (least squares in Minimum variance unbiased estimator for a linear signal model in presence of gaussian noise). Based on this model try and estimate the noise variance, the regularization term should be close to noise varainace.
Deploy machine learning techniques based on signal parameters such as dynamic range, frqeuency, FFT noise floor etc. to estimate noise
Follow a recursive approach with adjusting the regularization to minimize a suitable objective, ex: gradient descent etc.

In general your regularization is as good as how close it is to the actual noise or estimation of noise. This goes back to the idea of whether a frequentist approach is better or Bayesian approach. The answers is if the modelling of the parameter is accurate to it's actual behaviour then the Bayesian approach is better but if you model a parameter with a distribution but it is actually not following that closely then Bayesian would lead to suboptimal result and better go for a frequentist approach.