Video compression - why not just sending the changed parts?

You could do that, in theory, if you knew "what are the parts" (or objects) frame-to-frame in a video. But:

objects move, change scale, rotate, have color changes, they even can be occluded.
the segmentation is a complicated task in general.

MPEG 4, at least, allows to manipulate objects instead of pixels, see MPEG-4 Natural Video Coding - An overview. If objects and backgrounds are already decoupled, like filmed with a green screen:

then the object contents can be processed separately. Otherwise, what is more important is that, visually, the representation is accurate enough. Thus, a block from some object in a frame  can be cast to another block of a different object in a different frame, as long as the result is satisfactory. Thus, in a lot of cases, the block matching is sufficient, and more importantly, relatively fast to compute.


What is the "bilinear interpolation kernel" in personlab paper?

When I read «bilinear interpolation kernel» in image processing, I assume that they mean bilinear image imterpolation:
y(n+t) = (t-1)x(n) + tx(n+1)

If you decide upon a fixed uniform upsampling factor and standard dsp «zero-fill then lowpass filter» upsampling, the convolution kernel turns out to be samples of a triangular waveform.

-k


is there a way to apply a median filter diagonally?

The median filter can be applied on many pixel-set shapes, and indeed this was used to emulate of fasten full 2D median filters by a combination of leaner medians.

In Matlab, your request just amounts to gather diagonal elements in a vector. The simplest could be:

extract each sliding square patch of the image
extract the diagonal as a matrix with diag.m
compute the median

If you want the antidiagonal, flip the square patch with flipud, fliprl.
If you want other diagonals than the $\pm 45°$, you have to use discrete lines, for instance with the Bresenham algorithm.

If you want to remain high level, I'd suggest you to use gblk: a pedestrian data block extractor at MatlabCentral to parse an image with overlapping blocks.


what Graphical DSP programming tools are there for ASIO4ALL?

ASIO4ALL is a hardware independent low latency ASIO driver for WDM audio devices. As being a (wrapper) driver it does not support programming other but change of settings (device/channels selection, bit depth, sample rate, I/O buffers (latency)) through control software.
Regarding DSP programming tools, could you mean Creative/SoundBlaster specific kxproject driver/software?
Edit 1:
Creative/SoundBlaster DSP chips are listed here.
There are also software based methods as well as like EqualizerAPO  (system-wide equalizer for Windows) which lets you route channels, EQ filtering, convolution etc..


Scenes generally exhibit a fade-to-black transition. You could capture those frames with simple image processing tools, which are found in many libraries, such as OpenCV.

If you want to rely on the change of the mis-an-scene to be robust against sudden changes, then of cours algorithms for detecting temporal differences is more appropriate. For that, you might want to consider:

http://vis.uky.edu/~cheung/courses/ee639_fall04/readings/spie99.pdf
http://www.mitsubishi-electric-itce.fr/uk-rce/pubdocs/VIL07-D071.pdf

If you want a simple algorithm do as described here:
https://stackoverflow.com/questions/4801053/video-scene-detection-implementation

For harder cases, you might even try GMM-like background models to test if the transition was a temporarily quick shot or a complete change.


Is it true that each and every signal has multiple frequency components?

The only signal, that really has just one frequency component, is an infinite sine signal. Limiting the signal duration in any way is bound to produce other frequency components, as time limiting can be thought of as multiplying with a rectangle window, which translates to convolution with an si-function in the frequency domain, thus introducing new frequency components, however infinitesimal.

As for your example: the 60Hz AC power "signal" is messed up with all kinds of stuff coming from devices connected to the grid, that "backfire" into it. Also, the 60Hz are not really constant, they vary ever so slightly around 60Hz.


Adding phase to a signal that has multiple frequencies

Caveat: many behavioral models of amplifiers exist, none of which is perfect, and I don't know what is best for your particular amplifier.  But the most common model is basically your option 2.  I have never heard of anyone using option 1.

Let $x(t)$ be the input to the amplifier in complex baseband form.  A very common amplifier model is to assume that the gain is a memoryless function $g(|x(t)|)$ of the input amplitude, and the phase shift is also a memoryless function $\phi(|x(t)|)$ of the input amplitude.  The complex baseband output is therefore

$y(t) = g(|x(t)|)\exp(j\phi(|x(t)|))x(t) $

This model assumes the signal is relatively narrowband.

There is no need to decompose the signal into its Fourier coefficients to perform this transformation, and indeed this only complicates the matter.


It is not really a convolution, more a weighting. Your expression is $0$ homogeneous in $x$, as you divide $x$s values by $x$s.  Specifically:

it is not defined if $\sum_{-K,K} x[i+k] = 0$
if $x[k]$ is locally a constant, the output does not depend on it (because it will be equal to $\sum h[k]$)
if $h[k]$ is an averaging filter ($h[k] = \frac{1}{2K+1}$), the output will be $\frac{1}{2K+1}$

So globally I would not consider it as a  classical filter for $x$. But it could be seen as an adaptive weighting for $h$, based on a unit-normalized window $\chi$ with $\chi[n] = \frac{x[n]}{\sum_{|k|\le K} x[n]}$. As such, it is easily implemented in Matlab or Python. So far, I don't have a use for it. What was your purpose?


The more independent data you have, the more constrained are the possible solution sets that could produce that data, usually.  If any higher frequency content in the possible solution sets is constrained to not be completely arbitrary (which data derived from sub-pixel shifted sampling might so constrain), then the solution sets could possibly becomes sharper due to the steeper edges that can be reproduced using said higher frequency content.


Artificial Neural Network Example Problem

The tutorial states that each neuron uses the sigmoid activation function:

That funciton is:
$$
\sigma(x) = \frac{1}{1+ e^{-x}}
$$
so just substitute $x = 0.755$ into that to get 0.68026719669.


It is helpful to realize that the auto-correlation of $X[t]$ does not depend on $t$, just on the time difference $\tau$. So we can write it as $r_X(\tau)$. Your result, which is correct, can then be written as

$$E(Y[t+\tau]Y[t])=\frac14\left(r_X(\tau)+r_X(\tau+1)+r_X(\tau-1)+r_X(\tau)\right)\tag{1}$$

Using the expression for $r_X(\tau)$, which is correctly given in your question, $(1)$ can be simplified to

$$E(Y[t+\tau]Y[t])=\begin{cases}\frac{p}{2}(1+p),&\tau=0\\
\frac{p}{4}(1+3p),&\tau\in\{-1,1\}\\p^2,& \text{otherwise}\end{cases}$$

which is of course also independent of $t$.


What would be the effect of repeating samples in time domain on frequency domain?

It's quite benefical to view the mathematical expression of the process as well.

In discrete-time interpolation, a mere $N$ times repetition of a single sample (pixel) into the new enlarged block is an interpolation with zero-order hold kernel. This is mostly applied for simple fast but crude image zooming.

In 1D, the kernel is $$h_0[n] = \begin{cases}{ 1 ~~~,~~~ 0 \leq n \leq N-1 \\ 0 ~~~,~~~ \text{ otherwise } } \end{cases} \tag{1}$$

The interpolated signal can be expressed as :

$$y[n] = \sum_{k} x[k] h_0[n-kN] \tag{2}$$ this can be expressed as a convolution like

$$y[n] = \left(\sum_{k} x[k] \delta[n- Nk] \right) \star h_0[n] \tag{3}$$

where the signal in parenthesis is the expanded version of $x[n]$:

$$x[n] \rightarrow \boxed{ \uparrow N} \rightarrow x_e[n] = \sum_{k} x[k] \delta[n- Nk] \tag{4}$$

Applying DT Fourier transform on Eq(3) yields:

$$ Y(e^{j\omega}) = X_e(e^{j\omega}) \cdot H_0(e^{j\omega})  \tag{5} $$

where $X_e(e^{j\omega}) = X(e^{j\omega N}) $ due to the expansion operation in Eq(4) and $H_0(e^{j\omega})$ is the DTFT of the rectangular pulse $h_0[n]$ which is:

$$H_0(w) = e^{-j\omega (N-1)/2} \frac{ \sin( \omega N /2 ) } { \sin( \omega / 2 ) } $$

then the result becomes:

$$ Y(e^{j\omega}) = X(e^{j\omega N}) \cdot e^{-j\omega (N-1)/2} \frac{ \sin( \omega N /2 ) } {\sin( \omega / 2 ) }  \tag{6} $$

Whose magnitude is :
$$ \boxed{ |Y(e^{j\omega})| = |X(e^{j\omega N})| \cdot \left| \frac{ \sin( \omega N /2 ) } {\sin( \omega / 2 ) } \right| }  \tag{7} $$

Note that the $X(e^{j \omega N})$ is composed of repeated images of $X(e^{j \omega})$. And the effect of the interpolation kernel $h_0[n]$ is to remove the images as much as possible with that crude frequency response $H_0(e^{j \omega}) $.

The following Octave/Matlab code demonstrates the operation:

clear all; close all; clc;

N = 8;    % interpolation
L = 32;   % original signal length
x = 5*( sin(0.3*pi*([0:L-1]-L/2))./(pi*([0:L-1]-L/2))).^2;
x(isnan(x))=5*0.3*0.3;

h0 = ones(1,N);

xe = zeros(N*L,1);  % expanded signal
xe(1:N:N*L) =  x;

yi = conv(xe, h0);  % the interpolated signal: y[n] = xe[n] * h0[n]

figure,subplot(2,1,1);
stem(x);
title('original signal');
subplot(2,1,2);
stem(yi);
title('zero-order hold interpolated output');

figure,subplot(2,1,1);
plot(linspace(0,2,1024),abs(fft(x,1024)),'b');
title('Spectrum magnitude of x[n]')

subplot(2,1,2);
plot(linspace(0,2,1024),abs(fft(xe,1024)),'b');
hold on
plot(linspace(0,2,1024),abs(fft(h0,1024)),'r');
plot(linspace(0,2,1024),abs(fft(yi,1024)),'g','linewidth',2);
legend('|X_e(w)|','|H_0(w)|','|Y_i(w)|');

with the following outputs:

and the spectrum:

as you can see, the nulls of the periodic sinc (the Dirichlet kernel) are at the centers of the repeated spectra (the images) of $X(w)$. However this is not enough to suppress them all unlike an ideal lowpass filter (or one of its adequate approximations) which would completely eliminate the image spectra and yield only the one inside the main lobe; hence yielding a perfect interpolation compared to this crude zero order hold kernel giving.


Disagreement between Analytic Convolution and FFT Convolution

After playing with the notebook for a while, I figured out the following:

The result from scipy is shifted by one pixel in each axis (probably due to its definition of the "same" mode). Or, alternatively, both the analytic result and pyfftw results are shifted the opposite way.

After working around that, the result is exactly (within error) the same as the pyfftw result. I did this as follows (though this drops some data on the edges):

B_padded = np.pad(B, ((0, 1), (0, 1)), 'edge')[1:, 1:]
C1 = scipy.signal.convolve(A, B_padded, 'same', 'direct') * delta**2

The analytic solution (which peaks at 4.0) is not quite right for your input data (31x31 pixels at 1.0, edges at 0.5, and corners at 0.25). The peak value (where the image and the window align perfectly) should be:

box_bulk = 1.0**2 * (31 * 31)
box_edges = 0.5**2 * 31 * 4
box_corners = 0.25**2 * 4
peak = (box_bulk + box_edges + box_corners) / 1024 * 4  # ~= 3.876

which is the answer produced by both scipy and pyfftw.

Gist (also includes plotting tweaks).

(disclaimer: friend of asker)


Since you have frequency bucket information, you essentially have a frequency histogram.

If you treat histograms as n dimensional vectors, you can use dot product to get a similarity value that goes from 0 to 1 and could tune that so that sounds that have a value greater than say 0.6 are considered a match.

For more information check this out:  http://blog.demofox.org/2015/02/13/writing-a-basic-search-engine-aka-calculating-similarity-of-histograms-with-dot-product/

Alternately you could get a little more advanced and use support vector machines, which require training data of both negatives and positives to find a way to differentiate hits and misses.


Extended Kalman Filter for Linear Systems with Non Linear Measurements

When you start comparing sub optimal approaches, it usually comes down to trial and error.

A KF essentially operates on means and variances. An EKF does the same but the nonlinear measurement and or states necessitate some additional considerations.

The first thing one needs to determine is if the transformation results in a unimodal probability distribution where there is a well defined mean and variance.  If not, a particle filter will be a better solution. If the unimodal assumption is realistic, then the bias associated with how you map or update a mean is the next problem. In some cases a linearization is good enough and sometimes the transformation needs to be used and that includes things like Runge Kutta ODE solvers, which most text books fail to mention. In your case, a measurement nonlinearity can either use a transformation or a linearization based on a Taylor expansion to propagate a mean.

The final problem is how to propagate the variance.  Most EKF will use a Taylor Expansion most typically, just the first derivative, but there are second and higher order filters as well.  These enter as complexity considerations as well. The text edited by Gelb has a good number of covariance update heuristics in the latter chapters, including those for discontinuous transformation.

Gelb, Arthur, ed. Applied optimal estimation. MIT press, 1974.

The UKF is an attractive approach to covariant update.

Step size is also important. Your measurements may not be at regular time steps.

In summary, sub optimal is less than optimal but which suboptimal is better is going to depend on the application.  Essentially you need to try some alternatives and judge which is better.  There are a lot of alternatives and which is best isn’t easily predictable.


Binarization of image for varying background and foreground colors for OCR

This problem is well studied. It is termed scene text recognition. Of course, the state-of-the-art methods involve deep learning. In particular, I like the works, which use a bi-directional LSTM together with features extracted from Convolutional Neural Networks. This allows for end-to-end training and frees the segmentation or binarization requirement. In particular, I like the work:

  Reading Scene Text in Deep Convolutional Sequences -
  Pan He, Weilin Huang, Yu Qiao1, Chen Change Loy2, Xiaoou Tang2
  AAAI 2016


Compensating the Group Delay of an Analog Filter using DSP

It's often easier to design FIR filters for compensating group delay. At the same time they could also compensate the magnitude if necessary. The easiest method is to use a complex least squares method, which boils down to solving a system of linear equations for the filter coefficients.

The difficult part is to choose an appropriate bulk delay in the definition of the desired phase response. A good starting point is to choose the average delay to be close to half the filter length.

Take a look at this answer for an example of FIR group delay equalization.

You can use the Matlab/Octave function lslevin.m for designing FIR filters with prescribed magnitude and phase responses.


Comparing PSD two signals of different length

I am not interested in the common frequencies I am interested in the energy content.

You might be fooling yourself. Parseval's theorem clearly states that the energy in time is the same as the energy in frequency domain. Thus, if you're just after energy, your question's title makes little sense.

So instead of

samples -> DFT -> |·|² -> Σ

you can just

samples- > |·|² -> Σ

and get exactly the same result.

That clears up the question: If you're dividing the sum of magnitude square samples by the number of samples, you get an "Energy per time", which is power. That's a fair comparison for two signals. But then, you're comparing powers, not energies, and that might contradict your problem statement.

general advice: Taking a few guesses here, namely that you are a student relatively fresh into DSP, and you will need to hand in some kind of report or go through some kind of exam or lab demo on this stuff.

Just write down exactly what you're looking for, in mathematical notation. That will do two things:

it's much easier to spot such "terminology" mishaps like Energy vs Power, because you're suddenly sitting in front of a paper that says "looking for energy $P$" and you'll start to wonder why it is called $P$.
It's much easier to apply the things you've learnt in theoretical classes when your problem is already in the same "language". For example, I'm pretty optimistic that you've heard of energy and power signals, and even Parseval might have come up. Especially in such normalization questions like you have it, just seeing and discussing the formula helps a lot – because, later on, you will have the chance to define things yourself, and often, things are really easy to name and understand when  you see them as formula and go "oh, this is a sum of magnitude squares – let's call this energy", and vice versa, someone tells you about energy in images, and you'll instantly realize that an image is just a boring 2D signal, and thus, energy will just be "sum over all rows over all colums of color intensity squares".


How to downscale an image with interpolation by a factor larger than interpolation kernel size?

It is arguably inefficient. You could achieve the same results with a single processing step using a single kernel. You can do this by convolving all of the kernels together, then applying your resize using this larger kernel.

A similar process to this is used in most digital painting packages. The resizing usually offers some fixed kernel sizes (e.g. bilinear / bicubic - these may be referred to as 'smoother' / 'sharper') as well as some form of intelligent resize (sometimes called 'smart' or 'best quality') where the resampling kernel size is a factor of the image size.

If you need more information you should add some details about what precisely you are trying to achieve.


According to the properties of Z-transform, the ROC of sum of two signals is the intersection of their individual ROCs.

If you look at the signal
$$ x[n] = (\frac{1}{3})^n u[n] + (2)^n u[-n-1] $$

you see that it consists of two signals; one right sided $(\frac{1}{3})^n u[n]$ with a pole at $z = 1/3$ and whose Z-transform is $X(z) = 1/(1 - (1/3)z^{-1})$ and with region of convergence ROC1: $ |z| > 1/3$ and the other is a left sided one $(2)^n u[-n-1]$ with a pole at $z=2$ whose Z-transform is $-1/(1-2z^{-1})$ and with ROC2 as $ |z| < 2$. Therefore, the intersection of their ROCs is
$$ROC_x = ROC1 \cap ROC2 = \{ 1/3 < |z| \} \cap \{|z| < 2 \} = \{ 1/3 < |z| < 2 \} $$

similarly for the signal $y[n]$.


Working out integration of basic signal functions

So you've already figured out the $u(t-T)$ trick:  the unit step (delayed) just changes the limits of integration.

The other trick is that $\delta(t-T)$ just "samples" the value of the rest of the integral's argument. So, for example, $\int x(t) \delta(t-T) dt = x(T)$ where $T$ is somewhere within the integral's limits.


Confusion about how to calculate signal to noise ratio(SNR)?

The SNR is the ratio of the powers. Note that for zero-mean signals, the variance equals the power, so the SNR is the ratio of the variances of signal and noise, but only if both are zero-mean.


Why is a wavelet transform implemented as a filter bank?

First of all the basic idea of wavelet transforms lies in multi-resolution analysis. What this means is that the signal is looked at from different scales.

It is probably easier to understand this with images (which are 2D signals). The idea of multi-resolution is like zooming in and out of a reference signal and compared with windows of the image according to the amount of zoom.

Similarly for 1D you are trying to see how the mother wavelet at different scales (think of it as compressing and expanding in the 1D) compares with the signal at different delays (at different points on the discrete 1D).

Now when you look at a given set of discrete points you want to transform based on a mother wavelet it easy to see them as filters (based on mother wavelets) that are applied at different scales.

For a better understanding of wavelets and DWT I suggest you read Polikar's tutorial : Wavelet Tutorial


Filtering noise from FFT, where noise is known to be shifted delta function

High-pass filter the input with a zero-phase filter to eliminate the signal of interest, leaving only the high-passed noise. If the signal is periodic, you can do this in the frequency domain and convert back to time domain by IFFT. Find the time of the peak and deduce the magnitude of the unfiltered delta function from the magnitude of the filtered delta function (the peak value or root mean square). Knowing the delay and magnitude, subtract the noise from the original input signal. Do FFT.


Why do you get instantaneous phase by calculating correlation with a carrier signal?

This is the analytic signal. See wikipedia for better details.

Basically you are taking the carrier correlation, $f(x)$, and its Hilbert transform, $(h \ast f)(x)$ and combining them into a single complex valued signal, the analytic signal

$$f_A(x) = f(x) - i (h \ast f)(x)$$

which can be written as the complex exponential

$$f_A(x) = A(x) e^{i \phi(x)}$$

Why this is useful is that it performs a split of identity into seperate amplitude $A$ and phase $\phi$ components. Actually the Hilbert transform requires the entire signal to compute, but because you know the function of your signal already you only need the one sample.To localise the phase, you would use a bandpass filter.


Frequency estimates from three-dimensional data

I have [...] three-dimensional positional data [...], unevenly sampled at approx. 50 Hz, with timestamps.

As a first step, it would be good to interpolate them so that they appear to be sampled at regular $\frac{1}{50Hz}$ intervals.

  I am trying to find the dominant frequencies in the signal/tremor.

In this case, even a simple Discrete Fourier Transform (DFT) over a long enough segment (say for instance 10 times the longest period you are trying to discover) would suffice.

  However, the dominant frequency may not be constant.

In which case, you would have to consider something like a spectrogram or more generally, obtaining the frequency response over time (i.e. sliding windows).

The only "problem" that the Magnitude representation has is that $d = \sqrt{\sin(\theta_1)^2 + \sin(\theta_2)^2}$ is constant for the right choice of $\theta_1, \theta_2$, even if the individual components oscillate.

The danger of this happening in exactly that way depends on how the data are acquired but the point remains that by obtaining the magnitude you are introducing some interference in the measurements.

So, use with caution.

Now, the way that a multidimensional DFT works is by first of all assuming that the dimensions the signal is measured over are orthogonal and then (as a byproduct of that) repeatedly applying the transform to the "remaining" dimension.

In the one dimensional case, you have the way a quantity evolves in time, as a time series. The application of the DFT here is straightforward and it decomposes the quantity over time into a sum of sinusoids over time.

In the two dimensional case, you apply the DFT over the "rows" of a two dimensional matrix holding your "signal" and then once more over the columns of the already transformed rows from the previous step. The result of this process is a two dimensional spectrum where the equivalent frequency bin (from the one dimensional case) is now a frequency ring. And this is because in the two dimensional case, it is not enough to ask "which frequency... (?)", you also have to specify "...along a particular direction". For more information about this, please see here and here.

In the three dimensional case, you apply DFT to the "rows", you then apply DFT to the "columns" of the transformed "rows" and then you apply DFT once again along the "depth" rows (or, really, the remaining dimension) of the previously transformed data. This returns a spatial representation of your data where the "frequency ring" (which used to be the "frequency bin") is now a "frequency shell", that is, a hollow sphere. It is not enough to ask "which frequency..." now, you have to specify the direction on the surface of a sphere.

You may be wondering "so what?" by now and that is putting it mildly.

The point here is that if you do a three dimensional DFT you will also be able to infer the most dominant direction along which the tremor is happening too. And that might be "correlateable" with other parameters of the health condition. That is, different brain circuits deteriorating, leading to tremors along specific directions.

To do that, do a three dimensional DFT, shift the spectrum, so that the low frequencies are towards the centre of the described "cube", take the magnitude of the complex result and find the maximum value (that is not at DC). The maximum value will be at some point $m,n,k$. The "angle" between a vector denoting the "forward" direction and the $m,n,k$ vector would give you the direction that this movement is happening.

The tremor might be up-down, diagonal, back to front, circular, etc. A three dimensional DFT will characterise this periodic movement fully.

Hope this helps.


In CP-OFDM, the transmitted signal is the output of the IDFT appended with a cyclic prefix. No pulse-shaping is used.

In so-called "new waveforms" proposed as candidate for 5G, pulse shaping is heavily used.

In FBMC, each subcarrier (i.e each complex exponential) is convoluted with a pulse with a certain shape (usually a bell shape with beginning and end at zero amplitude). As a consequence of convolution, the amplitude of each filtered subcarrier at the beginning and at the end of the multicarrier symbol goes to zero.
Consequently, in case of multi-path, this zero-valued beginnings and ends avoid ISI, since the part of the previous symbol overlapping most with the current one is its end, which as we said is zero-valued (or almost) and the current symbol has little energy at its beginning (zero or almost).

Hope this helps.


take any of the pinking filters you can find and swap the zeros with the poles.

e.g. here's my pinking filter with 3 real poles and 3 real zeros (all inside the unit circle):

$$ H(z) = A \frac{(z-q_1)(z-q_2)(z-q_3)}{(z-p_1)(z-p_2)(z-p_3)} $$

$p_1$ = 0.99572754,     $q_1$ = 0.98443604

$p_2$ = 0.94790649,     $q_2$ = 0.83392334

$p_3$ = 0.53567505,     $q_3$ = 0.07568359

just swap the poles with the zeros.


Extracting amplitudes of components from FFT

Once you use parabolic or other interpolation for a frequency peak location estimate, you can use an offset Sinc kernel, centered at that estimated frequency, to interpolate the magnitude of that peak location.  Sum each of the complex components of all the adjacent component values of the FFT result that are above the local noise floor (weighted by the offset Sinc kernel).  Then take the magnitude of that complex result.

If you use any window function (other than rectangular of FFT width) use the transform of that window function instead of a Sinc kernel.


In digital image processing, an image resizing changes only the array size (dimensions) of the image. It does not care about dpi. It just upsamples or downsamples the image array.

DPI (Dot Per Inch) is considered only when the actual physical size of the displayed (or printed) image is of concern.

When you display a scaled image on the same physical monitor, then the apparent dpi of the image content will be changing too. This is especially important on text processing with font rendering.

Some programs try to keep image physical display size the same, while scaling them or when displayed on different monitors, by incorporating a dpi parameter to be used as well.

[edit for answer below]
In those examples from Gonzalez's book, images are considered with their physical sizes when displayed (or printed). Assuming square pixels, DPI of a physical display device is computed by
$$ \text{DPI} = \frac{ \text{ number of horizontal pixels on the monitor}}{ \text{ horizontal length (in inches) of the monitor}} $$

for example for a typical 15.6 inches, FHD (1920 x 1080 ), 16:9 aspect ratio laptop, the DPI will be
$$ \text{DPI} = \frac{ 1920 \times \sqrt{16^2 + 9^2}}{15.6 \times 16}  = 141.212 $$

That's the DPI of the monitor.

Physical size of a displayed image with a array size of $M \times N$ is given by
$$ \text{H} = \frac{ M }{ \text{DPI}} $$ horizontally and

$$ \text{W} = \frac{ N }{ \text{DPI}} $$ vertically in inches...

For amonitor with a given DPI, an image with smaller array size $M,N$ will display smaller on the screen. If you upsample (aka interpolate, or scale) the image array from $M,n$ to, say, $3M, 3N$, then the displayed image size will also be tripled.

And finally; interpolation of an image (or any digital data) means computing the sigal values that's physically in between its existing samples. In that respect you are decreasing the pixel distance when you interpolate. However, interpolation will not give you new data; it will only make exsting data smoother. That's exacty why you cannot get authentic FHD video simply by interpolating from an SD original source...


Which is the total energy of the product of two discrete energy signals

You don't give enough information for an expression of the combined energy in your expression.
$$\begin{eqnarray}
& &\sum_{i=0}^N \left| x_s(i)\cdot\left( 1+K\,w(i) \right) \right|^2 \\
&=&\sum_{i=0}^N \left( x_s(i)\cdot\left( 1+K\,w(i) \right) \right)^\star \times\left( x_s(i) \cdot\left( 1+K\,w(i) \right) \right)\\
&=&\sum_{i=0}^N \left( x_s(i)^\star\,x_s(i) + x_s(i)^\star\,K^\star\,w(i)^\star x_s(i)\,K\,w(i) + x_s(i)^\star x_s(i)\,K\,w(i)+x_s(i)\,x_s(i)^\star\,K^\star w(i)^\star\right)\\
&=&\sum_{i=0}^N \left( |x_s(i)|^2 \times \left(1+2\,\Re(K\,w(i))+|K|^2\,|w(i)|^2 \right) \right)\\
&=&E_s+\sum_{i=0}^N|x_s(i)|^2\left(2\,\Re(K\,w(i))+|K|^2\,|w(i)|^2 \right)
\end{eqnarray}$$

This means you need to know the correlation of the square magnitude of $x_s$ with the real part of $w$ and the square magnitude of $w$. The correlation of just $x_s$ and $w$ does not play any role.


How does the Meyr-Oerder Timing Algorithm work?

For PAM signals, the squaring creates a spectral component at $\dfrac{1}{T}$.  Consider the following trigonometric identity to understand why this is the case:

$$\cos^2\left(2\pi \dfrac{1}{2T}\right) = \dfrac{1}{2}+\dfrac{1}{2}\cos\left(2\pi \dfrac{1}{T}\right)$$

The phase of that spectral component at frequency $\dfrac{1}{T}$ is a measure of the evenness (or oddness) of the symbol relative to the reference time $t=0$ of the fourier transformed squared samples.  A phase of zero corresponds to the symbol being centered at that reference time.

The Fourier transform is used to find the phase of the $\dfrac{1}{T}$ spectral component, which is an unfiltered, scaled estimate of the symbol timing offset.


Show Equivalence Between Multiplication in Time Domain to Convolution in Frequency Domain

Let's assume you have 2 signals: vX and vY.
So:

clear();

numSamplesX = length(vX);
numSamplesY = length(vY);

numSamplesConv = numSamplesX + numSamplesY - 1;

vTimeDomainConv = conv(vX, vY);

vFrequencyDomainConv = ifft(fft(vX, numSamplesConv) .* fft(vY, numSamplesConv), 'symmetric');

max(abs(vTimeDomainConv - vFrequencyDomainConv)) %<! Should be < 1e-12

This is the Convolution Theorem for Discrete Signals to show convolution in time domain is equivalent to element wise multiplication in frequency domain.

If you want to show element wise multiplication in time domain can be done using the convolution in frequency domain you need to either interpolate the time domain signal to length of linear convolution or just use circular convolution in frequency domain:

clear();

numSamples = 6;

vX = randn(numSamples, 1);
vY = randn(numSamples, 1);

vXY = vX .* vY; %<! Time Domain Mulaiplication
vDxy = fft(vXY); %<! The DFT of the Time Domain Multiplication

vDx = fft(vX);
vDy = fft(vY);

% Applying circular convolution
vXYC = cconv(vDx, vDy, numSamples) / numSamples; %<! Normalization

% Comparing the result in Frequency Domain
max(abs(vXYC - vDxy)) %<! Should be < 1e-12

% Time Domain
vXYFromFrequency = ifft(vXYC, 'symmetric');

% Comparing the result in Time Domain
max(abs(vXY - vXYFromFrequency)) %<! Should be < 1e-12


Histogram matching (specification) in Python

If I understood right, you are stuck in matching a given histogram into a desired one and creating a new image from this matched histogram obtained by your filtering method.

I would first suggest you to get rid of all the unnecessary stuff (including the python code) above and isolate your problem as "histogram matching" in mathematical terms.

Now I will describe how you can approximately match a given histogram to a desired one in two setps: by first converting it into that of a uniform (equalized) one, and in then converting this uniform one to the desired one. Finally we shall combine these two steps to get the answer.

The method is based on converting a random variable into another by means of a transform $G(\cdot)$. You will replace the random variable $X$ with image intensity function $I(n,m)$.

Consider a uniform random variable $X$ whose CDF is $F_X(x) = U(x) = P(X < x) = x$.  We want to generate a new rv $Y$ from $X$, by means of a transform $Y=G(X)$, such that its CDF is the desired one. We shall therefore find this transform $G(x)$.

Under suitable conditions: $$F_Y(y)= P(Y< y) = P(G(X)<y) = P(X< G^{-1}(y)) = F_X(G^{-1}(y))=G^{-1}(y) $$

from which we deduce that the transform that converts uniform $X$ to arbitrary $Y$ with a desired CDF is $G(x) = F_Y^{-1}(x)$ from which we obtain: $$ \boxed{Y = F_Y^{-1}(X)}$$

Therefore the transform that converts a uniform RV $X$ to an arbitrary (desired) RV $Y$ is $$ \boxed{G(x) = F_Y^{-1}(x)}$$

The transform that converts an arbitrary RV $Y$ to a uniform RV $X$ is $$ \boxed{ G^{-1}(x)=F_Y(x) } $$

Finaly, the transform that converts arbitrary RV $Y$ to another arbitrary RV  $Z=T(Y)$ is $$ \boxed{ T(y)=F_Z^{-1}(F_Y(y)) } $$

We are almost done, we will finish by converting above algorithm to Image manipulation and also clarify how you will overcome the problem of not having those CDF's as formulas.

Replace RV $Y$ and $Z$ with i/o images $Iy(n,m)$, $Iz(n,m)$, also replace CDFs of $Y$, $Z$ with that of "computed from histograms". I denote it as $H_y(i)$ and $H_z(i)$ where $H_y(i)$ denotes the cumulative histogram obtained from image at hand, and $H_z(i)$ computed from the desired (to be matched) histogram. Under suitable conditions and correct scalings to ensure intensity levels to remain within valid limits, therefore we can make the definiton as $I_z = T(I_y)$ or more explicitly $I_z(n,m) = H_z^{-1}(H_y(I_y(n,m)))$

Your only problem is about the inverse cumulative histogram of $Z$ (or $I_z$) and unfortunately I don't know a method which may guarantee perfect inversion. As it is a CDF from a discrete histogram there might be empty bins and hence some indices of inverse $H_z(i)$ might have no mappings at all. It is up to you how to best overcome this problem.

However, as I understood from your code, you add pixels to empty/weak bins so that $H_z(i)$ becomes invertable...

Below I will put some matlab code "as is" : There are many methods to match histograms and I don't think below is the most accurate or the efficient one.

Iy = imread('cameraman.tif');  % Intesity between [0-255]
hy = imhist(Iy);         % compute image histogram hy(i)
figure,stem(hy);
title('original image histogram');

S = size(Iy);
Imax = 255;             % 8 bit
K = Imax / (S(1)*S(2)); % Scale factor
Hy = zeros(1,256);      % compute CDF Hy(i) of input image
Hy(1)=hy(1);
for i=2:256
  Hy(i) = Hy(i-1) + hy(i);
end
Hy = K*Hy;
figure,stem(Hy);
title('CDF of original input image');

% Get the desired image and histograms:
I_des = imread('tire.tif');
figure,imshow(I_des);
S = size(I_des);

h_des = imhist(I_des);
figure,stem(h_des);
title('desired image histogram');

h_des(h_des < (S(1)*S(2)/(Imax*8))) = h_des(h_des < (S(1)*S(2)/(Imax*8))) + 32; % HERE  paddings
Hz = zeros(1,256);  % This is desired CDF
Hz(1)= h_des(1);
for i=2:256
  Hz(i)=Hz(i-1)+h_des(i);
end
K_des = Imax / (Hz(256));
Hz = K_des*Hz;
figure,stem(Hz);
title('Desired CDF Hz');

% MATCH is performed
Imatch = (Iy);     % Processed image from matching
for i=0:255       % for each intensity level
 ind = (Iy==i);   % find those index with an intensity value of "i"

 j = Hy(i+1);

 looping=1;  %  implement you own method
 while looping
    k = find( (Hz > (j-1)) & (Hz < (j+1)));
    if length(k)>0   % for practical matters...
      looping =0;
    else
      j = j + 1;    % when there is no match, update a little
    end
 end

 Imatch(ind) = k(1); % Adjust the original intensity now...
end

figure, imshow(Imatch);
title('Image Matched to desired histogram');

h_match = imhist(Imatch);
figure, stem(h_match);
title(' Matched Histogram, not a good match ???');


Is it possible to reduce the complexity of radix-2 FFT if the input vector contains identical elements?

if

$$\text{FFT}\{[a, b, c, d]\} = [p, q, r, s]$$

then

$$\text{FFT}\{[a, b, c, d, a, b, c, d]\} = [p, 0, q, 0, r, 0, s, 0]$$

or
$$=[2p, 0, 2q, 0, 2r, 0, 2s, 0]\text,$$

depending on normalization.

Edit:
Using this property, we can make FFT nearly 2x faster on this specialized data. But the complexity is still $O(n \log n)$. I mean, big O notation doesn't count constant factor.


What are the keywords to find an algorithm that detects embossed text from a 3d model?

One approach would be to use 3D reconstruction and then use deformable ICP (iterative closest point) to recover the deformation map of the embossed text. The deformation map would transform the CAD model to the deformed scene. The inliers you compute during correspondence estimation would also allow you to segment these characters. (This might be similar to computing optical flow)

For reconstructing your scene, I would recommend photometric stereo techniques. They are pretty accurate in reconstructing small details. For large objects, they would not work. Some intro texts are here and here.

Also, in such cases, you could use an appropriate lighting to emphasize the text. A low angle ring light would be an example.


You can find a nice tutorial for time-frequency analysis in Numerical python by Johansson, chapter 17. link to github repository.

You can also check the scipy.signal.spectrogram.

import numpy as np
from scipy import signal
from scipy.fft import fftshift
import matplotlib.pyplot as plt

# Generate a test signal, a 2 Vrms sine wave whose frequency
# is slowly modulated around 3kHz, corrupted by white noise
# of exponentially decreasing magnitude sampled at 10 kHz.

fs = 1e4
N = 1e5
amp = 2 * np.sqrt(2)
noise_power = 0.01 * fs / 2
time = np.arange(N) / float(fs)
mod = 500 * np.cos(2 * np.pi * 0.25 * time)
carrier = amp * np.sin(2 * np.pi * 3e3 * time + mod)
noise = np.random.normal(scale=np.sqrt(noise_power), size=time.shape)
noise *= np.exp(-time / 5)
x = carrier + noise

fig, ax = plt.subplots(2, figsize=(8, 7))

f, t, Sxx = signal.spectrogram(x, fs)
ax[0].pcolormesh(t, f, Sxx)
ax[0].set_ylabel('Frequency [Hz]')
ax[0].set_xlabel('Time [sec]')

# Note, if using output that is not one sided, then use the following:
f, t, Sxx = signal.spectrogram(x, fs, return_onesided=False)
ax[1].pcolormesh(t, fftshift(f), fftshift(Sxx, axes=0))
ax[1].set_ylabel('Frequency [Hz]')
ax[1].set_xlabel('Time [sec]')
plt.savefig("fig.png")
plt.show()


The problem of polyphonic music transcription has received a huge amount of attention in the research community in recent years and I would go as far as to say that for single instrument polyphonies (piano, guitar etc...) the results are very good. Here are a few papers/authors who have looked at this problem deeply. Derry Fitzgerald has done a lot in the area, a lot of his NMF work on source separation produces accurate transcriptions. Anssi Klapuri has looked even more specifically at the problem of chord detection within multi instrument polyphonies. And finally, although not published, Mikel Gainza has developed very accurate chord transcription algorithms for commercial music soon to be released in a guitar based software product Riffstation.
The publications in the linkes here should give you a good idea of how the Polyphonic music transcription landscape currently lies.


Multidimensional FFT with mixed directions

if I understood it right, given your 3D volumetric data $x[m,n,p]$, first you want to take a 2D forward DFT (via 2D FFT) in the first two variables $m,n$ (i.e., for each m-n plane along 3rd dimension) and then take a 1D inverse DFT of vectors along the third dimension.

Note that in the first operation you replace each data plane with frequency plane and create a new volumetric data, then you will be taking inverse DFT of this resulting volume, the MATLAB/OCTAVE code below demonstrates a simple example:

A = zeros(2,2,3);      % create a 3D data array
A(:,:,1) = [1,2;3,4];  % first plane of 3D data
A(:,:,2) = [1,2;3,4];  % second plane of 3D data
A(:,:,3) = [1,2;3,4];  % third plane of 3D data

B = fft2(A)        % it'll take 2D-DFT independently for each plane of A

C = ifft(B,[],3)  % we indicate that ifft() be taken along 3rd dimension...

As you can see, since after the first step we have same frequency data for each plane of B, the third step produced impulses for each vector of C along 3rd dimension, due to inverse DFT.

I believe this is easy enough? But still in two steps (you can merge them into one line but still two calls)...


First note that:
$$ \cos(2\pi 50 t) \longleftrightarrow 0.5 \delta(f+50) + 0.5\delta(f-50) $$

$$\sin(2\pi 150 t) \longleftrightarrow 0.5 j \delta(f+150) -j 0.5\delta(f-150)$$

Hence the baseband spectrum is:
$$ X(f) = 0.5 j \delta(f+150) + 0.5 \delta(f+50) + 0.5 \delta(f-50) - 0.5 j\delta(f-150) $$

Then a shift to right by 100 Hz yields (without sampling weight):
$$ X(f-100) = 0.5\left[j \delta(f+50) +\delta(f-50) + \delta(f-150) - j\delta(f-250) \right]$$

Shift left by 100 Hz:
$$ X(f+100) = 0.5\left[ j\delta(f+250) +\delta(f+150) + \delta(f+50) - j\delta(f-50) \right]$$

Shift right by 200 Hz:
$$ X(f-200) = 0.5\left[ j \delta(f-50) +\delta(f-150) + \delta(f-250) - j\delta(f-350) \right]$$

Shift left by 200 Hz:
$$ X(f+200) = 0.5\left[ j\delta(f+350) +\delta(f+250) + \delta(f+150) - j\delta(f+50) \right]$$

The sum of all the shifts for $n=0,\pm 1,\pm2$ :
$$X_s(f) = 0.5\left[
j\delta(f+350) +\delta(f+250) + \delta(f+150) - j\delta(f+50) +
j\delta(f+250) +\delta(f+150) + \delta(f+50) - j\delta(f-50) +
j \delta(f+150) +\delta(f+50) + \delta(f-50) - j\delta(f-150)  +
j \delta(f+50) +\delta(f-50) + \delta(f-150) - j\delta(f-250) +
j \delta(f-50) +\delta(f-150) + \delta(f-250) - j\delta(f-350) \right]
$$

After the 100 Hz lowpass filter only those impulses inside $-100<f<100$ band remain:
$$LPF\{ X_s(f) \} = 0.5\left[ - j\delta(f+50) + \delta(f+50) - j\delta(f-50) +
\delta(f+50) + \delta(f-50) + j \delta(f+50) +\delta(f-50) + j \delta(f-50)
\right] $$

Which can be further simplied into:
$$LPF\{ X_s(f) \} =  \delta(f+50) + \delta(f-50)  $$

Hence we arrive (without sampling weight)
$$ \boxed{ z(t) = 2\cos(2 \pi 50 t)} $$

Therefore we can conclude that the output of the LPF will be proportional to $\cos(2\pi 50 t)$.


Non-linearly additive spectral subtraction?

It's difficult to suggest an algorithm without knowing more about the nature of these signals, but if it really is as simple as the above situation, then a spectral subtraction approach which uses the reference 'Single Substance' signal to activate a binary mask might do the job for you.


Interpolated FIR filter (from Oppenheim and Schafer's Discrete-Time Signal Processing, 3rd ed)

The reason why "interpolated FIR filter" is an appropriate name is because the filter $h_2[n]$ has $M_1-1$ zeros between its (non-zero) filter taps, and the convolution with the filter $h_1[n]$ interpolates those $M_1-1$ values.

Eq. $(104)$ generally only represents an FIR filter if both $h_1[n]$ and $h_2[n]$ are FIR filters.


The frequencies at FFT output are relative to the input sampling rate fs and the FFT size N. The $k^{th}$ index corresponds to a frequency $k*fs/ N$. Pay attention that values above $N/2$ correspond to negative frequencies


Estimate gap between two adjacent images

There are multiple papers dealing with this issue, using similar methods (extrapolation beyond the image edges to create a virtual overlap).

Alignment and Mosaicing of Non-Overlapping Images by Yair Poleg & Shmuel Peleg

Building a Mosaic from Non-Overlapping Images by Benjamin Choi

A Novel Technique for Non-overlapping Image
Mosaicing based on Pyramid Method by Achala Pandey & Umesh C. Pati


Parametric EQ vs multiple bandpass filters

a graphic EQ need not mean cascaded shelves.  it could mean cascaded peak/cut bell filters.  maybe with shelves as bookends, maybe not.

the bell filters are symmetrical in log freq (or the pitch scale) and are naturally spaced equally in log frequency.  you can draw any reasonable shape with a sum dB of these cascaded peak/cut filters.  but nothing too wild.

you could also do this more analytically with a sum of FIRs, all time-aligned.  make a series of overlapping complementary window functions (like a Hann) but operating on log frequency instead of linear frequency.  but that frequency response can be interpolated and sampled in linear frequency and an FIR can be derived from that magnitude response with a delay of some sufficient constant value.


Inverse DFT in terms of original signals for product of two DFTs

Again in this problem you can make use of the duality property to derive the output result without much involvement in the DFT sum formula. Before beginning, however, it's necesary to remember the time-frequency relations of two basic operations, namely, the expansion and compression blocks.

The EXPANSION block is the following:

$$ x_N[n] \longrightarrow \boxed{ \uparrow 2 } \longrightarrow y_{2N}[n] $$

shown operating on a time-domain signal $x_N[n]$ of length $N$ and producing the output $y_{2N}[n]$of length $2N$ given by the relation :

$$y_{2N}[n] = \begin{cases}{ x_N[n/2] ~~,~~~n=2m ~,~~~~~~ m=0,1,..,N-1 \\ ~~~~~0~~~~~~~~,~~~ n = 2m+1 ~, m=0,1,..,N-1} \end{cases} \tag{1}$$

The frequency domain relation between $N$-point DFT $X_n[k]$ of $x_N[n]$ and $2N$ point DFT $Y_{2N}[k]$ of $y_{2N}[n]$ is:

$$ Y_{2N}[k] = X_N[ (k)_N] ~~~,~~~ k = 0,1,2,...,2N-1 \tag{2}$$

Verbally stated, $Y_{2N}[k]$ is a cascade (mere repetition) of $X_N[k]$ and can be shown with the following matrix symbolism, each X,Y are (say) row vectors :

$$ Y = [ ~~~-~ X ~ -~~ | ~~~-~  X ~~~-~ ] \tag{3}$$

The COMPRESSION block is the following:

$$ x_N[n] \longrightarrow \boxed{ \downarrow 2 } \longrightarrow y_{N/2}[n] $$

shown operating on a time-domain signal $x_N[n]$ of length $N$ and producing the output $y_{N/2}[n]$of length $N/2$ given by the relation :

$$y_{2N}[n] = x_N[2n] ~~~, ~~~ n = 0,1,..,N/2-1 \tag{4}$$

The frequency domain relation between $N$-point DFT $X_N[k]$ of $x_N[n]$ and $N/2$ point DFT $Y_{N/2}[k]$ of $y_{N/2}[n]$ is:

$$ Y_{N/2}[k] = 0.5( X_N[ k ] + X_N[k + N/2]) ~~~,~~~ k = 0,1,2,...,N/2-1 \tag{5}$$

Verbally stated, $Y_{N/2}[k]$ is obtained from $X_N[k]$ by folding: adding second half of $X_N[k]$ onto its first half.

$$ Y = [ ~-~ {}_{N/2}X ~-~] + [~-~ X_{N/2} ~-~ ] \tag{6}$$

Then the DUALITY principle states that, if frequency domain sequences are processed by those expansion and compression blocks, then their effects in the time-domain seuqences will be the same as above that happened on the frequency-domain sequences of time-domain expanion and compression. Note that, when different sample rates (or number of DFT points) are used in the input and output of the processing blocks, which is the case for the expander and compressor examples, then the duality relation is scaled by the ratio of sampling rates which would cancel the 0.5 at the Eq(5) in our case.

By using these concepts, we can describe the time-frquency domain relation of the problem you have provided.

Now observe that your output signal $$Z[k]=Y(e^{j\omega})|_{\omega = \frac{2\pi}{N} 2k} = [Y(0) , Y(2), ... , Y(30) | Y(1) , Y(3) , ..., Y(29) ]$$ can be given by the following block diagram (with $N=31$):

$$ Y_N[k] \rightarrow \boxed{ Y_N -|- Y_N } \rightarrow W_{2N}[k] \rightarrow \boxed{ \downarrow 2} \rightarrow Z_N[k] $$

Notice that the first blcok is given in (3) and related to the expander operation, where as the second block is clearly a compressor. Lets denote the time domain signals by $y_N[n]$, $w_{2N}[n]$ and $z_N[n]$, and proceed from left to right as follows:

First, given the relation beween the DFT sequences $Y_N[k]$ and $W_{2N}[k]$ then their time-domain originators is the expander operation and hence we have from Eq(1):

$$w_{2N}[n] = \begin{cases}{ y_N[n/2] ~~,~~~n=2m ~,~~~~~~ m=0,1,..,N-1 \\ ~~~~~0~~~~~~~~,~~~ n = 2m+1 ~, m=0,1,..,N-1} \end{cases} \tag{7}$$

And secondly, given the compressor block applied on the frequency signal $W_{2N}[k]$, we conclude its effect, by using duality, on the time domain signal $w_{2N}[n]$ as:

$$ z_{N}[n] = w_{2N}[n] + w_{2N}[n+N/2] ~~~,~~~n=0,1,...,N-1 \tag{8}$$ namely the folding operation on $w_2[n]$.

Combining Eq(8) and Eq(9) we conclude the relation between N-point sequences $z[n]$ and $y[n]$ as:

$$\boxed{z[n] = \begin{cases}{ ~~~y_N[n/2] ~~~,~~~n=2m ~,~~~ m=0,1,..,N-1 \\ y_N[n/2 + N/2],~~~ n = 2m+1 } \end{cases} }$$

Eventhough the description seems really long and somewhat abstract, the operation is really neat and quick as the following Matlab/OCtave code demonstrates  it:

N = 31;
y = 1:N;   % just some signal

Y = fft(y,N);    % obtain YN[k]
W = [Y, Y];      % W = [ Y | Y] , w[n] = y[n/2]: is expanded verison of y[n]
Z = W(1:2:end);  % Z[k] = W[2k]; compress W[k] to get Z[k]

z = ifft( Z, N); % obtain z[n] from Z[k]

% Part-II: Use the time-domain relation as described in the answer to obtain z[n]
we(1:2:2*N) = y ;           % we[n] = y[n/2]: expanded version of y[n]
we = [we , 0];              % Append (2-1) zeros to the end of we if N is odd.
z2 = we(1:N) + we(N+1:2*N); % folding we[n] to obtain we[n]

% Display:
[z' , z2' ]


Skel to infinity morphological operations on 3D volume

As a start, you could look at the distance transform. The distance transform will assign each voxel the minimal distance to the volume boundary. Afterwards you could do a non-maximum suppression. The results should look similar to a skeleton.

I bet there are more sophisticated algorithms available for 3d skeleton extraction, but this is not my field of research, so I cannot help you with that.


I'm not sure I understand the delta-equals notation, but I'll assume it means definition.

It's important to understand that you're doing a variable substitution here.  In the strictest sense, the $t$ in $u(t-1)$ is not the same $t$ as in $u(t)$.  It may be easier for you to understand if you write your definitions as functions of $z$ instead of $t$:

$$ u(z) \triangleq
\begin{cases}
0, & \text{if }z<0 \\
1, & \text{if }z\ge0
\end{cases} $$

$$ \begin{align}
r(z) &  \triangleq
\begin{cases}
0, & \text{if }z<0 \\
z, & \text{if }z\ge0
\end{cases} \\
& = z \ u(z) \\
\end{align}$$

So, $z=t$ in $u(t)$ and $z=t-1$ in $r(t-1)$, etc.  You're changing a function of $z$ into a function of $t$.

The two functions you're dealing with have a breakpoint in them where the function changes behavior.  In both cases it happens when $z=0$.  So, looking at $u(z)$ for example, $u(z)=0$ for all values of $z$ less than zero and $u(z)=1$ for $z=0$ and all values greater than zero-- just saying what the function says, but in words.

When I use $u(t)$, I am substituting $z=t$ so $u(t)=0$ for all times before $t=0$, and $u(t)=1$ from time $t=0$ onward.

Now think about what's happening in the function argument with all the $(t-1)$, $(t-2)$ values...  If you were to look at $u(t-4)$, then, it has the effect of shifting the entire function to the right in time.  In this case:
$$ u(t-4) =
\begin{cases}
0, & \text{if }t-4<0 \\
1, & \text{if }t-4\ge0
\end{cases} $$
or, rearranging the conditions a bit:
$$ u(t-4) =
\begin{cases}
0, & \text{if }t<4 \\
1, & \text{if }t\ge4
\end{cases} $$
so you have a step function that is zero for all times before $t=4$, and one for times $t=4$ onward.

With the $r(z)$ function, with $z=t-2$ you substitute for $z$ in all places, so you have:
$$ \begin{align}
r(t-2) &  =
\begin{cases}
0, & \text{if }t<2 \\
t-2, & \text{if }t\ge 2
\end{cases} \\
& = (t-2) \cdot u(t-2) \\
\end{align}$$
which is zero for all times before $t=2$ and then a line starting at $(t=2,x(t)=0)$ with a slope of positive 1 and a y-intercept of -2.

The coefficients have the effect of changing the slope and corresponding y-intercept but the t-intercept remains the same.  $2\cdot r(t-2)$ has a slope of 2 and a y-intercept of -4 but still crosses the $t$ axis at $t=2$.

The x(t) function is a linear combination of such functions, so you must add up the values at each point.  The easiest way to see this on graph paper would be to draw each term separately one above the other and then add them vertically.

With the two components you you've highlighted, $−r(t−1)+2\cdot r(t−2)$:  There is a slope of zero before $t=1$ because $u(t)$ is zero slope.  The first term starts a ramp of slope of -1 starting at $t=1$ and the second term adds to that a ramp of slope +2 starting at $t=2$.  Where those two terms overlap (everywhere from $t=2$ and to the right) the combined slope of those terms is +1.


Why are 2nd order edge detectors such as Laplacian more susceptible to noise?

Generally 2nd derivative is more sensitive to noise than 1st derivative. The 2nd derivative is usually accompanied by zero crossing detection, so it works better when grey level transitions are smooth. When the noise level is high, the gray level fluctuates and may result in many zero crossings.

Besides, there is less chance of multiple edge responses with 2nd derivative edge detectors.


Echo Cancellation over Voice Telephone Network

A usual arrangement of echo cancelers in a full-duplex system is shown in the following figure (from Digital Communication by Lee and Messerschmitt):

The echo cancelers try to remove any residual transmitted signal in the receiver of the same end. Since the input to the echo canceler is the transmitted signal, it can only remove signal components in the received signal that are correlated with the transmitted signal. It cannot remove any other signal components, such as the signal that is transmitted from the far end (which is the desired received signal). Consequently, the error signal can never become zero, and the desired signal (transmitted from the other end) will not be influenced by the adaptive filter.


This sum appears quite often in DSP.

\begin{align}
\sum_{n=0}^{N-1} \exp(-j\alpha n) &\stackrel{(a)}{=} \frac{1- \exp(-j\alpha N)}{1 - \exp(-j\alpha )}\\
&= \frac{e^{-j\alpha N/2}(e^{+j\alpha N/2}- e^{-j\alpha N/2})}{e^{-j\alpha /2}(e^{+j\alpha /2}- e^{-j\alpha /2})}\\
&\stackrel{(b)}{=} e^{-j\frac{\alpha}{2}(N-1)} \frac{\sin(N\alpha/2)}{\sin(\alpha/2)}
\end{align}
where the last step $(b)$ comes from identity $e^{+j\theta}-e^{-j\theta} = 2j\sin(\theta)$ and be careful at $(a)$ as Matt's comment.

Substitute $\alpha=k\omega_0=k\frac{2\pi}{N}$, $\sin(N\alpha/2)=\sin(k\pi)=0$ if $k\in\mathbb{Z}$ hence so is the sum.

Update for $\alpha=0$

It is trivial that $\exp(-j\alpha n) \stackrel{\alpha=0}{=} 1$ then the sum equals $N$.

It can be thought as the case of the final result after $(b)$  when $\alpha \to 0$ then
$$e^{-j\frac{\alpha}{2}(N-1)} \frac{\sin(N\alpha/2)}{\sin(\alpha/2)} \stackrel{\alpha \to 0}{\to} N$$
where the limit is evaluated by L'Hôpital's rule.

Further interesting discussion can be found in Laurent Duval's answer.


Getting the bitstream of an unknown RF signal

First I found the sine wave extrema (black markers at the bottom) and stretched and superimposed a sine wave graph (red) on top of the signal graph with the extrema in the same locations. The signal matches (0) the sine wave over some of the non-zero segments, and over the other non-zero segments it matches the sine wave with a sign flip (1). This gives a binary sequence.

A sign flip is equivalent to a 180 degree phase shift, so this is a form of binary phase shift keying (BPSK, 2PSK, also phase reversal keying, PRK) with the useful symbol interval duration (of each sinusoid piece) seemingly equal to the guard interval duration (of each zero-valued segment that separates symbols). It is unusual that here the symbol rate is higher than the sine wave frequency, with a ratio of about $20/7\approx2.86$.


Could modulation/coding schemes be generated from an AI?

As per my knowledge, one of the recent papers on this topic can be found here, where the authors used machine learning algorithms to generate the optimal constellation shaping. Surprisingly, the results were in line with an old work: F. R. Kschischang and S. Pasupathy, "Optimal nonuniform signaling for Gaussian channels," in IEEE Transactions on Information Theory, vol. 39, no. 3, pp. 913-929, May 1993.

Moreover, autoencoders, neural networks, reinforcement learning and deep learning are being applied for end-to-end learning (especially in the case of wireless communications), where the mathematical model of the wireless communication channel is not very accurate. One such article can be found here.

Also, for the case of polar codes, deep-learning-based decoders are very popular nowadays at least in academia.


concurrent multiple sound sources seperation for TDoA estimation

Shuo,

Here is a couple of insights for your system:

1.- The different sound sources must play different signals (beacons), that may be resolved at the receiver device. A simple example would be to transmit a set of tones (one by speaker) with different frequencies. A slightly more complex solution would be to transmit other kind of signals with low cross-correlation (e.g. Gold, Kasami or Zadoff-Chu sequences).

A different approach at all would be to transmit actual modulated messages that somehow don't interfere each other at the receiver. This is usually even more complex.

2.- Regarding the method for distinguishing the different signals at the receiver, the typical solution is to use a set of matched filters (one filter for each beacon signal). Then, determine the TDoA by measuring the differences between the time indexes in which you find peaks at the matched filters' outputs.

On the other hand, if you really need to use modulation/encoding, a simple set of matched filter won't be enough, you will probably need to implement actual demodulation schemes, which may include synchronization algorithms and other practical restrictions.

According to my personal experience, and considering your introduction lacking much detail, you may successfully use a set of (frequency separated) tones at the transmitter, and a set of matched filters at the receiver. At least as a first approach, and if the performance is not good enough for you, you may use more advanced techniques. I honestly believe that in most cases no further cumbersome techniques are required.

I hope this helps!


You can compute the instantaneous frequency of your data in Matlab as follows:

z = hilbert(x);
N = length(z);
f = angle(z(2:N).*conj(z(1:N-1)))/2/pi;

To get units of Hertz, multiply f by the sample rate.

There are other, similar estimators that do a similar thing. This matlab code can help.

Eric Jacobsen also has a couple of pages here and here.


How to set covariance for Bivariate Logistic Distribution

The covariance of any multivariate distribution is always a matrix. Maybe it is of help to you to know what this matrix contains:

-In general, each entry i,j of the matrix gives you the correlation between the RVs i and j.

-So a diagonal element i,i will simply be the variance of the RV i.

-If RV i and j are independent, then the non-diagonal entry i,j is zero (they are uncorrelated). The reverse of this statement is not generally true, but it is for Gaussians iirc.

I don't know what you are trying to do, but perhaps the trace of the matrix or its determinant is of use to you?


STFT window width vs. detectable Hz range

Regarding the first paragraph of your question: The number of samples used in a discrete Fourier transform (DFT) depends on your desired DFT bin spacing (the frequency difference between adjacent DFT bins) in the frequency domain. If that desired bin spacing is Fo Hz, and the time-domain signal sample rate is Fs samples/second, then the number time samples, N, is:

N = Fs/Fo samples.             (1)

For example, at Fs = 44100, if you select N = 4410 samples then the DFT bin spacing will be:

Fo = Fs/N = 44100/4410 = 10 Hz.        (2)

However, for the traditional radix-2 fast Fourier transform (FFT) N must be an integer power of two. So using the above Eq. (2) when N = 2^12 = 4096 then Fo = 44100/4096 = 10.766 Hz.

Regarding the second paragraph of your question: I don't fully understand that paragraph but if Fs = 44100, someone performing 1024-point FFTs (N = 1024) produces spectra where the FFT bin spacing is:

Fo = Fs/N = 44100/1024 = 43.06 Hz.        (3)

Does that Fo = 43.06 Hz FFT bin spacing seem correct based on other information you see from the NMF chapter?


Orthogonal Basis for a 2D Signals (Compressive Sensing)

If you have a function $ f \left[ m, n \right] \in \mathbb{R}^{M \times N} $ then the DFT of those functions is an orthogonal basis of functions in $ \mathbb{C}^{M \times N} $.

So if we have:

$$ F \left[ k , l \right] = \frac{1}{\sqrt{M N}} \sum_{m = 0}^{M - 1} \sum_{n = 0}^{N - 1} f \left[ m, n \right] {e}^{-j 2 \pi \left( \frac{k}{M} m + \frac{l}{N} n \right)} $$

So the set of $ \left\{ g \left[ m, n \right] = \frac{1}{\sqrt{M N}} {e}^{-j 2 \pi \left( \frac{k}{M} m + \frac{l}{N} n \right)} \right\}_{k, l = 0}^{M - 1, N - 1} $ is the orthogonal basis (Actually Orthonormal).


I'm not a mathematician but if the multiple composition $ f_1(f_2(f_3 \dots)))$
is equivalent to a sequence of single compositions $g_2(g_1)$, and another $g_3(g_2)$ and so forth.

The generalization looks reasonable for  $x, f_{1}, f_{2} \dots f_N \in \cal{R}$ and $f_i$ differentiable with no discontinuities.


What transformations have been applied to this image to yield four derivatives?

This to me looks exactly like the two-dimensional Haar wavelet transform. All four quadrants are a part of the transform, even though different operations have been performed on them. Upper left quadrant is a result of low-pass filtering with a 2x2 two-dimensional rect function (Haar wavelet) and decimation be a factor of 2, upper right is the output of a horizontal high-pass filter, lower left is the output of a vertical high-pass filter, and lower right is the output of a diagonal high-pass filter.


I think it's instructive to try to understand how this error function is derived and why it makes sense. First, the desired impulse response $h_d[n]$ is given implicitly by the input-output relation of the desired system:
$$y[n]=\sum_{k=1}^{5}x[n-k]=\sum_{k=-\infty}^{\infty}h_d[k]x[n-k]\tag{1}$$
The given error function is the error function minimized by Prony's method for the design of IIR filters. We try to approximate a given transfer function $H_d(z)$ by an IIR filter $H(z)=B(z)/A(z)$:
$$B(z)\stackrel{!}{=}H_d(z)A(z)\tag{2}$$
For the given example we get (with $B(z)=1$)
$$1\stackrel{!}{=}\left(h_d[0]+h_d[1]z^{-1}+\ldots+h_d[N-1]z^{-(N-1)}\right)\left(1+a_1z^{-1}+z_2z^{-2}\right)\tag{3}$$
For equation $(3)$ to be satisfied, we need all coefficients associated with negative powers of $z$ to vanish:
$$\begin{align}h_d[1]+h_d[0]a_1&\stackrel{!}{=}0\\
h_d[2]+h_d[1]a_1+h_d[0]a_2&\stackrel{!}{=}0\\
h_d[3]+h_d[2]a_1+h_d[1]a_2&\stackrel{!}{=}0\\\vdots\end{align}\tag{4}$$
In practice, we can solve $(4)$ in an approximate way by minimizing the sum of the squares of the left-hand side, leading to the given error function
$$E=\sum_{n=-\infty}^{\infty}\big(h_d[n]+h_d[n-1]a_1+h_d[n-2]a_2\big)^2\tag{5}$$
Finally, the optimal coefficients $a_1$ and $a_2$ are obtained by taking the derivative of $(5)$ w.r.t. $a_1$ and $a_2$ and equating them to zero. This results in two linear equations with two unknowns.


Following up on my deleted answer... If you take a filled ellipse and project all the points onto the $x$ axis, more points will be projected near the origin than on the extremes, in a circular-shaped distribution. Not a gaussian distribution, and not the uniform distribution I mentioned in the 1-D analogy in my deleted answer. The resulting distribution actually has pdf $p(x) = \sqrt{(1 - (\frac{x}{r})^2)}$, and from there you can compute that the standard deviation is $\frac{r}{2}$.

Thus, if data is uniformly distributed in the interior of an ellipse of radii $a, b$ (whose axes are the $x$ and $y$ axes), the standard deviation of the $x$ coordinate is $\frac{a}{2}$ and of the $y$ coordinate is $\frac{b}{2}$. So the correction factor you need to use is simply 2.

Here is a worked example in python recovering the center (translation matrix), rotation matrix, and radii of an ellipse from points randomly sampled from its interior:

import numpy

# Generate some points distributed uniformely along an ellipse
x = 2 * (numpy.random.rand(10000, 1) - 0.5)
y = 2 * (numpy.random.rand(10000, 1) - 0.5)
d = (x / 0.5) ** 2 + (y / 0.25) ** 2
inside = numpy.where(d < 1)[0]
x = x[inside]
y = y[inside]
data = numpy.hstack((x, y)).T

# Now rotate by 0.5 rad and translate it to (4, -8)
angle = 0.5
rotation = numpy.array([
    [numpy.cos(0.4), -numpy.sin(0.4)],
    [numpy.sin(0.4), numpy.cos(0.4)]])

data = numpy.dot(rotation, data)
data[0, :] += 4
data[1, :] -= 8

# Step 1: center the data to get the translation vector.
print 'Translation', data.mean(axis=1)
data -= numpy.reshape(data.mean(axis=1), (2, 1))

# Step 2: perform PCA to find rotation matrix.
scatter = numpy.dot(data, data.T)
eigenvalues, transform = numpy.linalg.eig(scatter)
print 'Rotation matrix', transform

# Step 3: Rotate back the data and compute radii.
# You can also get the radii from smaller to bigger
# with 2 / numpy.sqrt(eigenvalues)
rotated = numpy.dot(numpy.linalg.inv(transform), data)
print 'Radii', 2 * rotated.std(axis=1)


Bilinear transform with pre-warping for systems other than classical filters

The bilinear transform ("Tustin's method") is indeed mainly used for transforming frequency selective filters with magnitude responses that are optimal with respect to some criterium, such as Butterworth, Chebyshev, or Cauer filters.

For more general systems, the bilinear transform is usually not the best choice because of the frequency warping, which mainly affects frequency domain behavior close to Nyquist. For frequency-selective filters, cut-off frequencies can be pre-warped such that the cut-off frequencies of the resulting discrete-time system are correct, but for more general systems it is usually not clear which frequencies should be matched and which shouldn't. Also, for filters with approximately piece-wise constant magnitude responses, frequency warping is not a big problem as long as the cut-off frequencies are mapped correctly. For general systems one would usually prefer the error between the continuous-time system and the discrete-time system to be distributed more evenly in the frequency domain, which is not possible with the bilinear transform.

There are several methods for converting continuous-time systems to discrete-time systems. All of them have advantages and disadvantages, and it depends on the application which one is to be preferred. Take a look at Matlab's c2d function, which implements several methods. Also take a look at this question and its answers.

For approximating a measured response, you would directly design the corresponding system in the discrete-time domain. For general frequency responses, the frequency sampling method is very simple and effective. If an IIR approximation is desired, the Matlab function invfreqz could be used. It tries to fit a rational transfer function to given frequency response data.


What are features in signal processing using a more tangible or less abstract description

A feature is a number that describes one aspect of a signal.
Signals can be very complex, and the simplest analysis tools (like a time plot, a spectrum, or an energy measurement) don't tell you everything; in fact, for specific types of analyses, they almost don't tell you anything useful.
So, features are designed to describe very specific aspects of a signal.
To make an analogy, let's say you're an expert in ice cream and want to compare two brands. Just looking at a picture or the list of ingredients is not enough. You can taste them, but you can't accurately describe the taste to your colleagues (just like you can listen to sound but you can't describe exactly what it sounds like).
So, you come up with "features" about the ice cream: the density, the statistics of the nitrogen bubbles inside it, the sugar/fat ratio, the crystalinity, whatever. Some of the features will be useless, but some will turn out to be really good at objectively conveying something interesting about ice cream.


Absolute convergence of periodic sinc interpolation

I'm top-editing this since it answers the question directly.
The sinc series is fundamentally a $C/x$, so you can extract as many absolutely convergent series out of it as you want, but what is left over is still only conditionally convergent.  Also, you can rescale $x$ and it is still a $C/x$ series.
Saying you have a summation to or from infinity is an informality.  Formally, you have a finite sum to some value, and take the limit as that value goes to infinity.
Therefore, your first and second series should have been done like this:
$$ \lim_{L \to \infty} \sum_{k=-L}^{L} f(k) =  \lim_{L \to \infty} \left[  \sum_{k=0}^{L} f(-k) + f(k+1) \right]  $$
$$ =  \lim_{L \to \infty} \left[ \sum_{k=0}^{L} \left( f(-k) + f(2k+1) + f(2k+2) \right) + \sum_{k=0}^{L+1} f(-k-L-1) \right] $$
Likewise, your third should have added this:
$$ \sum_{k=0}^{L+1} \left( f(-k-L-1) + f(-k-2L-3) \right) $$
Sometimes it takes a while to get around to where you should have been in the first place.  I'm deleting the rest.  Whoever is curious can find it in the edit history.

Proceeding informally....
First rearrange it:
$$
\begin{aligned}
y_m &= \sum_{n=0}^{N-1} x[n] \sum_{k=-\infty}^{\infty}  \operatorname{sinc} \left( \frac{Nm}{M} - n - Nk \right) \\
&= \sum_{n=0}^{N-1} x[n] W_m[n]
\end{aligned}
$$
One way to look at that is a resampled value is a linear combination (weighted average) of the sample points.
Another way is that you now have $N$ separate infinite series, all of the form:
$$
\begin{aligned}
W_m[n] &=  \sum_{k=-\infty}^{\infty}  \operatorname{sinc} \left( \frac{Nm}{M} - n - Nk \right) \\
 &=  \sum_{k=-\infty}^{\infty}  \frac{ \sin \left( ( Nm/M - n - Nk ) \pi\right) }{ (Nm/M - n - Nk) \pi }  \\
\end{aligned}
$$
Even $N$ Case:
$$ W_m[n] =  \sin \left( ( Nm/M - n ) \pi\right) \sum_{k=-\infty}^{\infty}  \frac{ 1 }{ (Nm/M - n - Nk) \pi }  $$
Odd $N$ Case:
$$ W_m[n] =  \sin \left( ( Nm/M - n ) \pi\right) \sum_{k=-\infty}^{\infty}  \frac{ (-1)^k }{ (Nm/M - n - Nk) \pi }  $$
Clearly, both are cases of $C/x$ series and not absolutely convergent.  If $Nm/M$ is an integer all the terms are zero except perhaps the zeroth terms.
As for the second comment, if I remember correctly (and I've already proven I didn't remember well), doing it formally does away with all the rearrangement tricks.  And yes, If I remember correctly, absolutely convergent series are immune to rearrangement tricks.
This too:
A series converges if and only if the sequence of partial sums converges.
A sequence converges if and only if for any given $\epsilon$ there exist a $\delta$ so for every $k > \delta$ the absolute value of the difference of the limit and the sequence value is less than $\epsilon$.
Stamp it on your forehead for formal occasions.
Disclaimer:  Been a long time ...

As clearly as I think I can say it:
The only conditions for which the series in Olli's Eq (1) will converge absolutely is when all the terms heading towards infinity are zero, since then their absolute values are zero.  This happens when all the $x_n$ are zero (the trivial solution) or $Nm/M$ is an integer.  Both the even and odd cases under any different conditions can be rearranged to be summations of alternating monotonically decreasing sequences, therefore they converge conditionally since they diverge absolutely.

Epilogue:
There is no need to do the infinite summation at all.  Direct closed form expressions exist for the odd and then even case based on the interpolation functions found when considering an inverse DFT as a continuous function.  The derivation of the functions can be found in the epilogue of my answer here:
How to get Fourier coefficients to draw any shape using DFT?
The derivation is based on the definitions of the DFT, the inverse DFT, and a finite geometric summation.
Resampling the continuous function at $M$ evenly spaced (in the cycle domain) points can be done by a simple variable substitution.
$$ t = \frac{m}{M} 2\pi $$
The direct sample set to sample set equations are then as follows.
Odd case:
$$ y_m = \sum_{n=0}^{N-1} x[n]
\left[ \frac{   \sin \left( N \left( \frac{m}{M} - \frac{n}{N} \right) \pi \right) }
            { N \sin \left(   \left( \frac{m}{M} - \frac{n}{N} \right) \pi \right) }
\right] $$
Even case, evenly split Nyquist bin:
$$ y_m = \sum_{n=0}^{N-1} x[n]
\left[ \frac{   \sin \left( N \left( \frac{m}{M} - \frac{n}{N} \right) \pi \right) }
            { N \sin \left(   \left( \frac{m}{M} - \frac{n}{N} \right) \pi \right) }
\right]         \cos \left(   \left( \frac{m}{M} - \frac{n}{N} \right) \pi \right) $$
These are mathematically equivalent to taking the DFT of size $N$, zero padding it at the Nyquist frequency to size $M$ (splitting the Nyquist bin in the even case), then taking the inverse DFT to recover a $M$ point upsampled sequence.  All the upsampled points lie on the underlying continuous interpolation function no matter what the point count.
For the $N=2$, $M=4$ case:
$$
\begin{aligned}
y_0 &= x_0 ( 1 ) + x_1 ( 0 ) = x_0 \\
y_1 &= x_0 \left( \frac{ \sin( \pi / 2 ) }{ 2 \sin( \pi / 4 ) } \cos( \pi / 4 ) \right) + x_1 \left( \frac{ \sin( -\pi / 2 ) }{ 2 \sin( -\pi / 4 ) } \cos( -\pi / 4 ) \right) \\
&= \frac{1}{2} ( x_0 + x_1 ) \\
y_2 &= x_0 ( 0 ) + x_1 ( 1 ) = x_1 \\
y_3 &= x_0 \left( \frac{ \sin( 3 \pi / 2 ) }{ 2 \sin( 3 \pi / 4 ) } \cos( 3 \pi / 4 ) \right) + x_1 \left( \frac{ \sin( \pi / 2 ) }{ 2 \sin( \pi / 4 ) } \cos( \pi / 4 ) \right) \\
&= \frac{1}{2} ( x_0 + x_1 )
\end{aligned}
$$
Which should be the results you are expecting.
An infinite number of sinc functions can now take the day off.

Suppose that instead of doing halfsies on the Nyquist bin we apportioned them as $(1/2+g)$ and $(1/2-g)$, this would alter the continuous interpolation function as follows.
$$
\begin{aligned}
D(t_n) &= \left( \frac{1}{2} + g \right) e^{i(N/2) t_n } + \left( \frac{1}{2} - g \right) e^{i(-N/2) t_n } + \sum_{l=0}^{N-2} e^{i ( l - N/2 + 1 ) t_n } \\
&= \cos \left( \frac{N}{2} t_n \right) + i 2 g\sin \left( \frac{N}{2} t_n \right)  + \frac{  \sin( t_n  N /2 ) } { \sin( t_n / 2 ) } \cos( t_n / 2 ) - \cos( t_n  N /2 )  \\
&= \frac{  \sin( N t_n/2 ) }{ \sin( t_n / 2 ) } \cos( t_n / 2 ) + i 2g\sin \left( \frac{N}{2} t_n \right)
\end{aligned}
$$
The extra term introduced is purely imaginary.  That can be folded in, but I prefer to leave it separate when put back into the function definition.
$$
\begin{aligned}
z(t) &= \sum_{n=0}^{N-1} x[n] \left[ \frac{ \sin( N (t - \frac{n}{N}2\pi) / 2 ) } { N \tan( (t - \frac{n}{N}2\pi) / 2 ) } + i \frac{2g}{N}\sin \left( N (t - \frac{n}{N}2\pi) / 2  \right)  \right]
\end{aligned}
$$
If is obvious that any non-zero value of $g$ will add "energy" to the signal, thus the $g=0$ solution, corresponding to halfsies on the Nyquist bin, is the most natural solution, or lowest energy, out of a whole family of solutions of periodic bandlimited at $N/2$ functions.
The more significant convincer for me is that it also introduces imaginary values into what is other wise fully real set of weighting values.
Whether R B-J's series converges uniquely to this "natural" solution, or the "natural solution" is a unique solution (it is not) are two totally separate issues.

Olli, I hope this makes you smile.
Start with the discrete resampling formula for the odd $N$ case.
$$ y_m = \sum_{n=0}^{N-1} x[n]
\left[ \frac{   \sin \left( N \left( \frac{m}{M} - \frac{n}{N} \right) \pi \right) }
            { N \sin \left(   \left( \frac{m}{M} - \frac{n}{N} \right) \pi \right) }
\right] $$
Since the sequence of $N$ points is periodic $( x[n] = x[n+N] )$ and all the points are covered, we can shift the summation range to be zero centered.
$$ L = (N-1) / 2 $$
Also, the $m$th point can be located on the $n$ scale.
$$ w = m \frac{N}{M} = \frac{m}{M} N $$
Since the $M$ resampled points are evenly spaced along the cycle, they too can be arbitrarily shifted to be zero centered, though strictly that is not necessary.
Since "$t$" has already been used above, the scale of the domain of the continuous interpolation function, both will get new names.  "$z(t)$" and "$Y(\omega)$" describe the same function.  Plug all the defined values in.
$$
\begin{aligned}
y_m = Y(w) &= \sum_{n=-L}^{L} x[n]
\left[ \frac{   \sin \left(             \left( w - n \right) \pi \right) }
            { N \sin \left( \frac{1}{N} \left( w - n \right) \pi \right) }
\right] \\
&= \sum_{n=-L}^{L} x[n]
\left[
\frac{\frac{\sin \left( \left( w - n \right) \pi \right)}{ \left( \omega - n \right) \pi }}
{\frac{\sin \left( \frac{1}{N} \left( w - n \right) \pi \right)}{\frac{1}{N} \left( w - n \right) \pi }}
\right] \\
&= \sum_{n=-L}^{L} x[n]
\left[
\frac{\operatorname{sinc} \left( w - n \right) }
     {\operatorname{sinc} \left( \frac{1}{N} \left( w - n \right) \right)}
\right] \\
\end{aligned}
$$
Now it's time to take the big step, that is, big stroll out to infinity.  The cycle of $N$ points grows until one cycle spans negative to positive infinity.  As it gets bigger, the circular nature gets more remote.
$$
\begin{aligned}
\lim_{N \to \infty} y_m &= \lim_{N \to \infty}  Y(w) \\
&= \lim_{N \to \infty} \sum_{n=-L}^{L} x[n]
\left[
\frac{\operatorname{sinc} \left( w - n \right) }
     {\operatorname{sinc} \left( \frac{1}{N} \left( w - n \right) \right)}
\right] \\
&=  \sum_{n=-\infty}^{\infty} x[n]
\left[
\frac{\operatorname{sinc} \left( w - n \right) }
     {1}
\right] \\
&=  \sum_{n=-\infty}^{\infty} x[n] \operatorname{sinc} \left( w - n \right) \\
&=  \sum_{n=-\infty}^{\infty} x[n] \operatorname{sinc} \left( \frac{Nm}{M} - n \right)
\end{aligned}
$$
Now look at that.  The Whittaker–Shannon interpolation formula has been derived from scratch and we are right at your starting point.
The even case can be done similarly and ends up with the same formula.

Definition of DFT of $N$ samples
Inverse DFT used as Fourier Series Coefficients for interpolation function
Dirichlet Kernel form of interpolation function
Interpolation function used for $M$ samples
Even and Odd Discrete Weighted Average Resampling Formulas
N goes to infinity
Whittaker–Shannon emerges
Whittaker–Shannon applied to a repeating sequence of $N$
Convergence questioned

I hope realizing using step 7 to achieve what step 2 has already answered will put a smile on R B-J as well.  Your proof lies there.
For $ N = 2 $
$$
\begin{aligned}
y_m &= \sum_{n=0}^{1} x[n]
\left[ \frac{   \sin \left( 2 \left( \frac{m}{M} - \frac{n}{2} \right) \pi \right) }
            { 2 \sin \left(   \left( \frac{m}{M} - \frac{n}{2} \right) \pi \right) }
\right]         \cos \left(   \left( \frac{m}{M} - \frac{n}{2} \right) \pi \right) \\
&= \sum_{n=0}^{1} x[n] \cos^2 \left( \left( \frac{m}{M} - \frac{n}{2} \right) \pi \right) \\
&= x_0 \cos^2 \left( \frac{m}{M} \pi \right) + x_1 \sin^2 \left( \frac{m}{M} \pi \right)
\end{aligned}
$$
For $ x_0 = 1 $ and $ x_1 = -1 $
$$
\begin{aligned}
y_m &=  \cos^2 \left( \frac{m}{M} \pi \right) - \sin^2 \left( \frac{m}{M} \pi \right) \\
&=  \cos \left( \frac{m}{M} 2 \pi \right)
\end{aligned}
$$
I'm going to have to be done with this for a while.  Neat stuff.

Olli, thanks for the bounty points.
This little exercise has deepened my understanding of W-S considerably.  I hope that is true for you and Robert (and others) too.
It is still a precarious foundation though.  I wanted to convince myself that it would work for a sinusoid of any frequency.  To wit:
$$ x[n] = M \cos( \alpha n + \phi ) $$
$$
\begin{aligned}
x(t) &= \sum_{n=-\infty}^{\infty} x[n] \operatorname{sinc}(t-n) \\
&= \sum_{n=-\infty}^{\infty} M \cos( \alpha n + \phi ) \operatorname{sinc}(t-n) \\
&= \sum_{n=-\infty}^{\infty} M \cos( \alpha t + \phi - \alpha( t - n ) ) \operatorname{sinc}(t-n) \\
&= \sum_{n=-\infty}^{\infty} M \left[ \cos( \alpha t + \phi ) \cos( \alpha( t - n ) ) + \sin( \alpha t + \phi ) \sin( \alpha( t - n ) ) \right] \operatorname{sinc}(t-n) \\
&= M \cos( \alpha t + \phi ) \sum_{n=-\infty}^{\infty} \cos( \alpha( t - n ) ) \operatorname{sinc}(t-n) \\
& \qquad \qquad + M \sin( \alpha t + \phi ) \sum_{n=-\infty}^{\infty} \sin( \alpha( t - n ) ) \operatorname{sinc}(t-n) \\
&= M \cos( \alpha t + \phi ) \cos( \alpha( t - t ) ) + M \sin( \alpha t + \phi ) \sin( \alpha( t - t ) ) \\
&= M \cos( \alpha t + \phi ) \cdot 1  + M \sin( \alpha t + \phi )  \cdot 0  \\
&= M \cos( \alpha t + \phi )
\end{aligned}
$$
I seem to have accomplished my goal.  However, there is nothing in this proof that prohibits $\alpha \ge \pi$, though that is a condition for the validity of the theorem.  So, knowing that, you are okay.  If you didn't know that, the formula itself does not reveal it.  To me, that's troubling.

Reply to R B-J:
First off, no where is it stipulated that $x[n]$ must be real.  Even for a real valued function, you don't have to split the Nyquist bin halfsies to get a real interpolation function.  Just pick $g$ to be a multiple of $i$ above.
Suppose you have the function:
$$ z(\tau) = \sum_{k=-L}^{L} c_k e^{ik\tau} $$
Its band limit is $L$ or less.  Every k term, except 0, can be paired up with it's conjugate bin and the sum can be decomposed into a cosine and sine term.
let $ A = \frac{c_k + c_{-k}}{2} $ and $ B = \frac{c_k - c_{-k}}{2} $
$$
\begin{aligned}
c_k e^{ik\tau} + c_{-k} e^{-ik\tau} &= (A+B) e^{ik\tau} + (A-B) e^{-ik\tau} \\
&= 2A \cos(\tau) - i 2B \sin(\tau)
\end{aligned}
$$
For a regular bin, we can only say $X[k] = c_k$ if $k+N>L$, otherwise I have more than one k in the bin and cannot separate them.  At the Nyquist bin $X[k] = c_k + c_{-k}$
Think in terms of degrees of freedom.  For a complex signal, $c_k + c_{-k}$ has four and the Nyquist bin two.  Therefore there are two free.  Just enough to put a complex parameter on the Sine function at the Nyquist frequency.  With a real signal, $c_k + c_{-k}$ has two degrees of freedom and the Nyquist bin value restricts one of those leaving one left over.  Just enough for a real valued parameter times the Sine function to remain a real valued signal.
I showed earlier the translation between not doing halfsies and the consequence on the interpolation function.  Nothing prohibits that and it doesn't increase the bandwidth of the solution one iota.

R B-J asks:
//  //"But we do know A will be zero in the halfsies and W-S reconstructions."// how do you know that?  //
The halfsies is easy.  Without loss of generality, consider the $N=2$ case.
$$ x[n] = [1,-1] $$
$$ \frac{1}{N} X[k] = [0,1] $$
Halfsies on the Nyquist of 1.  Doing an unfurled inverse DFT with split Nyquist:
$$ x[n] = \frac{1}{2} e^{i\pi n} + \frac{1}{2} e^{-i\pi n} = \cos(\pi n) $$
Now allow $n$ to be real, call it $t$ to indicate the change.  This defines an interpolation function (still called $x$).
$$ x(t) = \cos(\pi t) $$
For every other even N, the unnormalized DFT will be (0,0,0,....,N), so the result remains the same.
For the W-S summation, look at the section where "omega" temporarily lived, the "Sinc is the limit of the Dirichlet Kernel" section.  The left side $y_m=Y(w)$ is known to be $ \cos( \pi w ) $.  I even did the specific $N=2$ case after the dependency list.  Just set "M=2" which makes $w = m$.  The limit reached at the end of the second chunk gives your summation.  Just reverse the order of the equation and you get:
$$ \sum_{n=-\infty}^{\infty} (-1)^n \operatorname{sinc}( w - n ) = \cos( \pi w ) $$
The fact that your summation is the limit of something is why proving it differently has been hard.
I think your time reversal argument is good, too.  The sampled points are time reversible on the discrete $n$ scale, but that does not mean the source x(t) is, but it does mean Y(w) is.
P.S.  From now on, when a fresh context can be established, I'm going to use $\tau$ for a $ 0 \to 2\pi $ cycle scale, $t$ to be on the sampling scale ($=n$).


I think you have in general the correct idea. But you have to remember that the encode order (and therefore the order in which the pictures will appear in the bitstream and be decoded), depends on how the references are.

if you have: I(1) B(2) B(3) P(4) B(5) B(6) P(7). etc...

You could have a structure where in the group [B(2) B(3) P(4)], the first B(2) has a reference to the past I frame, but a future reference to the following B(3) frame. In that case, the encode order start the same as having I(1) P(4), but the next cannot be B(2), because B(2) needs to wait for B(3) to be encoded. And this could change from repetition to repetition.

So in general the encode order of frames depends on how the referencing structure is. Since the standards do not specify this and allow freedom of choice when encoding, the GOP structure is in general totally arbitrary, so every specific case has to be analyzed.


Inverse Fourier transform - porting from MATLAB to Python

For a vector Y of length $N$, where $N$ is even, conjugate symmetry requires the following property to hold:
$$
Y[k] = Y[N-k]^*
$$
for $k=1,...,N-1$ and $Y_0$ is real. (See this.) Here $^*$ denotes the complex conjugate.

The array [1,2,3,3,2,1] this conjugate symmetry property doesn't hold. What did Matlab do when you used 'symmetric' with an input vector that wasn't conjugate symmetric? It pretends [1,2,3,3,2,1] is nearly conjugate symmetric$^\dagger$ and just replaces it with the closest conjugate symmetric vector, which is, [1,2,3,3,3,2].

If you try this in Python:

Y = np.array([1,2,3,3,3,2])
y = np.fft.ifft(Y)

your answer will match Matlab's answer with the 'symmetric' option for [1,2,3,3,2,1]. Note there there will be some rounding issues and you can verify that the imaginary part of y in Python is practically zero. Here's my output:

In [8]: np.fft.ifft(np.array([1,2,3,3,3,2]))
Out[8]:
array([  2.33333333e+00 +0.00000000e+00j,
        -5.00000000e-01 -1.48029737e-16j,
        -1.66666667e-01 +0.00000000e+00j,
         1.11022302e-16 +2.90348171e-16j,
        -1.66666667e-01 +0.00000000e+00j,  -5.00000000e-01 -1.42318434e-16j])

In [9]: np.real(np.fft.ifft(np.array([1,2,3,3,3,2])))
Out[9]:
array([  2.33333333e+00,  -5.00000000e-01,  -1.66666667e-01,
         1.11022302e-16,  -1.66666667e-01,  -5.00000000e-01])

In [10]: np.imag(np.fft.ifft(np.array([1,2,3,3,3,2])))
Out[10]:
array([  0.00000000e+00,  -1.48029737e-16,   0.00000000e+00,
         2.90348171e-16,   0.00000000e+00,  -1.42318434e-16])

And now to your main question: How do you port this to Python?

Step 1: Take the input vector Y=[1,2,3,3,2,1]

Step 2: Force it to be conjugate symmetric by looping through and changing the right half of the vector to be equal to the complex conjugate of the left half. Call this modified vector Y_new = [1,2,3,3,3,2].

Step 3: Compute np.fft.ifft(Y_new) and throw away its imaginary part. (You can do an assert and ensure the imaginary part has a very tiny norm.)

Alternatively, when $N$ is even, Steps 2 and 3 can be replaced by a quicker np.fft.irfft(Y[0:N/2+1]) as endolith suggests in the comment below.

See also: https://blogs.mathworks.com/steve/2010/07/16/complex-surprises-from-fft/ for an example when $N$ is odd.

$^\dagger$ AFAIK, mathematically, there is no such thing is "nearly conjugate symmetric." A vector is either conjugate symmetric or not! I wish Matlab's documentation was more explicit about what it is doing under to hood.


Apparently it's normalized so that the gain at Nyqusit is 0 dB. Potentially they just normalize to the maximum being 0 dB. Gain at Nyquist is given by

$$H(\omega = \pi) = \frac{b_0-b_1+b_2}{a_0-a_1+a_2}$$

whereas gain at DC is

$$H(\omega = 0) = \frac{b_0+b_1+b_2}{a_0+a_1+a_2}$$

Both values are real.


Warning: $|e^{j\omega}|$ is equal to $1$ if and only if $\omega$ is a real number. More generally, for $z$ complex,  $|e^{z}|=e^{\Re{z}}$.

This really depends in your prior knowledge, because functions like exponentials, sines or cosines, can be developed in different ways. Here, I'll assume you know about sines  and cosines, as you refer to angular frequency. Then, an historical way to define the exponential is via Euler formula, for instance:

$$e^{j\omega} = \cos \omega+j\sin \omega\,.$$

Then, through Argand's complex plane interpretation, you can interpret a complex number as a point with  2D coordinates $(x,y) \leftrightarrow  x+jy$: the real part refers to the $X$-axis, the imaginary one to the $Y$-axis. Thus, obviously, $e^{j\omega}$ has the coordinates of a point on the $(0,0)$-center circle with radius $1$, hence the magnitude  $|e^{j\omega}|$ is equal to the radius, or  $(\cos \omega)^2+(\sin \omega)^2=1^2$.

In more modern versions, $e^z$ is defined by a series (see Proof of complex conjugate symmetry property of DFT), and $\cos$ and $\sin$ are derived for it, so they don't come first.


Match (OCR?) a digital image of a product label to similar

First, you will need an image processing algorithm to detect local features of the digital images you have. 1 and 2 provide good reviews of such algorithms.:

After extracting the local features, a matching process between the extracted features and your set of labels is used to link the labels with acquired images. The second references goes over few matching techniques and performance metrics.


bibliography about implementation of IIR filters

Okay, according to the wikipedia article you cited (with an appropriate substitution in notation):

$$ \Big|H(j2\pi f)\Big|^2= 10^{\frac{2.00}{10}} \times \frac{12194^4\cdot f^8}{\big(f^2+20.6^2 \big)^2(f^2+107.7^2)(f^2+737.9^2)\big(f^2+12194^2\big)^2}$$

and the dB amplitude curve is this:

$$A(f)=10\log_{10}\left(\Big|H(j2\pi f)\Big|^2\right) $$

I hope you can see how I modified the language and that it is equivalent.  The purpose is to express this transfer function as an $H(s)$ and get where (and how many) the poles as zeros are.

Now you gotta factor some more shit out.  Oh, this is laborious and painful.  Sorry, unless someone wants to finish this answer, I will have to return to it.


Adding white noise to complex signal (complex envelope)

If you want a complex white noise sample with variance $N_0$, then you should generate two independent noise samples each with variance $\frac{N_0}{2}$. The two samples that you generate make up the real and imaginary parts of the desired complex noise sample. The resulting complex value will have the desired variance $N_0$. You should then sum this with your complex baseband signal to get the effect that you want.


What is the relation of the transition band's width and the filter order for the FIR windowing method

There are only heuristic formulas for estimating the filter order. For a Kaiser window (which is probably the most frequently used window for filter design) the required filter order can be estimated from [1]

$$M=\frac{A-8}{2.285\,\Delta\omega}\tag{1}$$

where $A=-20\log\delta$ ($\delta$ is the maximum deviation from the desired response), and $\Delta\omega$ is the (smallest) transition bandwidth. This formula is of course only valid for the approximation of ideal frequency selective filters (low pass, high pass, etc.). This formula is implemented in Matlab's kaiserord.m function.

[1] Oppenheim, Schafer, Buck, Discrete-Time Signal Processing, 2nd ed., p. 476


How to find peak frequency in an audio frequency spectrum?

The way of choosing the threshold is based off of what is called constant false alarm rate (CFAR) energy detection. The idea is that as the user you choose a desired probability of false alarm $\text{P}_{fa}$ to calculate what the threshold should be. If you have signal plus Gaussian noise:
$$ y[n] = x[n]+w[n]$$
where $w[n] \sim N(0, \sigma_w^2)$ and $\mathbb{E}\big(\big|x[n]\big|^2\big)=\sigma_x^2$. The CFAR detector takes $N$ samples of $y[n]$ as input, computes the energy $\frac{1}{N}\sum_{n=1}^N |y[n]|^2$, and compares it to a threshold to determine if there was signal present in that chunk. You can just as easily make the extension for this to work on the FFT, as you will have:
\begin{align}
\text{FFT}(y[n])&= \text{FFT}(x[n])+\text{FFT}(w[n]) \\
Y[k] &= X[k] + W[k]
\end{align}
The choice of threshold value has been worked out in many papers on energy detection, but for the case described here it is:
$$ \text{Threshold} = \bigg( \frac{Q^{-1}(\text{P}_{fa})}{\sqrt{N}} +1\bigg) \sigma_w^2$$
This intuitively should make sense. The threshold should depend on the noise power (if $\sigma_w^2$ increases, we need to increase the threshold to keep a constant $\text{P}_{fa}$), the number of samples (remember law of large numbers), and $\text{P}_{fa}$ (you already said that you observed how changing the threshold changes the false alarms).
There are tradeoffs to be had when choosing the parameters. For example, you could choose a very large $N$, but this would give you very little frequency resolution. Saying that there is a signal present between 17 kHz and 18 kHz is not very helpful for you (in fact there are multiple signals there), so you need to choose the $N$ small enough to get the resolution you want but not too small that you need a super high SNR for this to work properly. Perhaps choosing a moderate sized $N$, followed by a peak finder might be robust.
Aside: your signal has structure (periodic signals), use that to your advantage! Looking at the FFT isn't really giving you anything new, just looking at your signal in a different way. Take a look at the FFT of the autocorrelation function, you should have even cleaner peaks there.


What is the difference between monocular and stereo image?

Have you ever wondered why people and many animals have two eyes? By using two eyes we can have a perception of depth. The principle is very simple, while you are looking the screen keep your head fix and close one of your eyes, try to watch all the things behind your screen, now open your closed eye and close the other one and try to watch the same things, can you see the difference?

Each eye send information to our brain. This information is almost the same, if you could have two pictures of your both eyes and you try to overlap this pictures you will see that there are big displacements in the nearest objects to your eyes while the displacements in the objects that are far away will be very small. These observed displacements are translated to the depth of the objects.

Now other experiment that could help you to understand this concept is this: Try to keep your view fix close your eyes for a while, now open only one eye but try to don't move your view, try to touch the things around you, but remember keep your view fix, is it simple?. Now with one of your eyes closed, move your head and try to do the same, is it easier?. By moving your head, you are giving information to your brain of how far are the objects because in the nearest objects the observable displacements will be bigger. Of course you should consider that the eyes have also other kind of cells that help us to perceive depth. But in a simple stereo pair of cameras the information that we could recuperate is only the distance of the same point in each image, the disparity map, and there are a lot of ways to recuperate it, optical flow is very common. This kind of information is normally used to know the 3D space and the segmentation is applied comparing the the distance between the cameras and the objects.

While by using only a camera you can compare only the intensities changes and try to detect the main changes by using a lot of metrics. Is very common to use this kind of systems in tracking of objects by keeping the camera fix and trying to omit the things that aren't moving in the scene.


doesn't OFDM inherently deals with ISI? why cyclic prefix is needed?

OFDM (sub)carriers are only exactly orthogonal over a certain IFFT block length (and integer multiples).  Multipath echos will mix a subcarrier signal from the end of one symbol block into the beginning of the next, rendering the symbol in that carrier to be ambiguous over the initial one IFFT block length.  A Cyclic Prefix extends the length of each symbol beyond just that of one IFFT length so that the IFFT can be (hopefully) done over a portion of the signal past where all the orthogonal subcarriers have been messed up with multipath echoes.  Skipping past the messed up portion allows the (constant or otherwise) channel response to be mostly ignored by the receiver, while maintaining the orthogonality of all the subcarriers from each other.