What function to apply a Metamorphopsia (grid warp) translation?

What you are looking for is called a "Geometric Transformation".

When applied to image processing, it comes down to a re-mapping of pixel positions. In optical systems, they can be used to represent certain distortions caused by lenses. In fact, ray transfer matrices end up being applied as geometric transformations to images, since they represent exactly how is a pixel displaced from its original position as its light ray travels within an optical system.

There are two things that need a little bit more consideration with geometric transformations:

How to express the mapping of pixels; and
How to apply the interpolation.

Expressing the mapping of the pixels might be relatively "easy" for known distortions (as mentioned earlier) but what is depicted in the body of the question here seems to be arbitrary (?). In any case, you are going to have to somehow express the function that takes a pixel from location $x,y$ to location $x',y'$ in a way that makes sense to the problem at hand. The way that image processing software (such as Illustrator that you mention) does it, is by setting up a grid that divides the image into a set of patches. These lines can be piecewise linear or bezier curves. This "control" (in terms of a user interface element) is not exactly straightforward to build.

Interpolation is absolutely necessary here because the geometric transformation might "ask" for pixel values between known pixels. There are obvious ways out of this problem (e.g. linear interpolation) but it is very likely that you will need a different method that makes sense within your application.

In Python, geometric transformations can be handled by Scipy's geometric_transformation(). Its use is straightforward provided that you know the input and output pixel locations.

Hope this helps.

EDIT:

  How would I go about remapping these pixels considering I don't want singular pixels remapped but the entire areas 'pulled' in a direction.

This is what the geometric transformation does, as long as the input/output point pairs are provided.

  Will interpolation work with colour in my image?

Yes.

  Do you have any suggestions where I could learn how to do this practically? Arbitrary can come after I've sorted this first part out!

Yes, have a look at the links i provided. From the mathematical point of view, we are talking about a series of linear transforms that map a "rectangle" to a distorted "rectangle". See the documentation of the geometric_transformation().


Compute real signal from a discrete analytic signal

@user85858302: If designing a complex-valued bandpass filter is the solution to your problem, you can design such a filter without using the Hilbert transform. Start by using your favorite method for designing a real-valued N-tap FIR lowpass filter whose one-sided bandwidth is B Hz. That will give you N real-valued coefficients. Next, multiply those coefficients by a complex exponential whose frequency is Fc Hz. (Fc must be sufficiently less than half your data’s Fs sample rate.) That will produce N complex-valued coefficients that can be used in a tapped-delay line FIR filter block diagram. The complex-valued bandpass filter will have frequency response whose center frequency is Fc Hz and a bandwidth of 2B Hz.  Model this whole process in software to demonstrate what I'm sayin' here.

For those readers interested in both real-valued and complex-valued FIR filters, I suggest you have look at the following blog. You may well learn something new about FIR filters. http://www.dsprelated.com/showarticle/808.php


Extracting ridges in automatically in image

This should be doable with ssqueezepy's extract_ridges; try varying penalty and bw (see their docstrings). As last resort, feeding a cropped image that excludes region without ridges may work better, as the algorithm assumes the ridge spans the entire frame. You can automate this by finding indices at which column energies fall below a set threshold, e.g. np.where(np.sum(np.abs(stft_transform)**2, axis=0) < 0.1).
'Cropped' example:


Mathematical proof of relationship between derivative and convolution

There's no mathematical proof per se; the fact is that finite differences are one method of approximating derivatives. Remember the definition of a derivative:

$$
f'(x) = \lim_{\Delta x \to 0} \frac{f(x+\Delta x) - f(x)}{\Delta x}
$$

The above equation states that to calculate the derivative of a function $f$ at position $x$, one takes the difference between the function evaluated at $x$ and the function evaluated at some small distance $\Delta x$ away. The ratio of the difference in the function value between the two points and the distance traveled along the $x$-axis is defined to be the derivative, as the distance between the two points becomes infinitesimally small.

A finite-difference derivative method, then, is just a way of approximating the above relationship using some finite spacing instead of the infinitely-small spacing implied by the limit. This type of approach makes sense for real-world applications where you might not analytically know the function $f(x)$  (and therefore can't evaluate its behavior across arbitrarily-small intervals.

With that said, there are a number of ways to approximate the derivative of a signal in this way; one example is the causal first-order difference:

$$
x'(t) \approx \frac{x(t) - x(t-T)}{T}
$$

This is simple to implement, but performs best if $T$ is small relative to the inverse of the signal's bandwidth. A slightly more complicated method can provide better results in some cases (the second-order difference):

$$
x'(t) \approx \frac{x(t) - 2x(t-T) +x(t-2T)}{T^2}
$$

I've used the "backward second-order difference" form above to show how they could be applied using causal filters, but note that they will have a resulting delay commensurate with the order of the approximation.


Are Fast Wavelet Transform coefficients constant if the input signal frequency coefficients are constant?

[Specific answer] If you have a real sine or cosine, the discretization and finite support may entail issues in the possibility of having a constant: edge samples or periods that cannot be divided by powers-of-two can affect the  wavelet coefficients.
But more important: let us suppose the discrete sine period is well chosen. Then a low-pass or high pass filter on a sine will yield something close to a sine (because complex exponentials are invariant under linear systems). So even after downsamling, the coefficients in a wavelet subband are likely to behave as sinusoidal values.
If you expect "constant values", you are likely looking at a complexified and absolute-valued representation of the coefficients.
With a real signal and discrete wavelet, one will likely see ripples on several subbbands, instead of a constant.
[Generic answer] The purpose of most transforms (whether fast or not) is to help in handling or interpreting the content of some data, or classes of data. Invertible transformations keep all information, sometimes with redundancy.
Discrete wavelet transforms (DWT) are meant to capture information in a non-redundant way from not-so-stationary signals. They often  are really bad as capturing information from a purely periodic sine signal.
As most wavelet filter bank filters are not perfect, and the signal is of finite length, filtering and representing coefficients will suffer artifacts.  So in a DWT, a sine would spread on many subbands, and suffer from aliasing, which is well known. However, despite the above caveat, the magic of wavelets can however recover the original sine, from the spread and aliased coefficients.
So :

don't expect DWT coefficients of a sine to be constant
don't expect them correspond to the true frequency content of a signal

They are not meant for that. However, if a signal has constant frequency content, but phase shifts (thus non-stationary), then  wavelets (maybe continuous once) can be great at detecting their location or onsets.


Cancel Drift after numerical integration

Numerical integration may reduce noise under some conditions: the noise is zero-average, ergodic, its properties do not vary too much over time. Since numerical integration seems linear, pre-filtering the acceleration with a linear filter, or post-filtering after integration are likely to yield very similar results (up to border effects), since linear operations commute.

The separation of a useful signal from a complicated drift is a complicated topic, without further information of the signal, the noise and the drift. These are a few venues you can use, alone or combined, to improve your signal:

filter the acceleration with more efficient (nonlinear) filters. For signals with time dynamics, Kalman filters may embed integration drift (Kalman filter: Modeling integration drift). In econometrics, people use a family of Holt-Winters filters to smooth or extrapolate time series with noise, drift or seasonality.
use more advanced numerical integration methods, like Simpson that can filter more noise or taking into account the dynamics, see for instance in Direct numerical integration methods, the Wilson Theta method:

  the Wilson Theta method assumes that the acceleration of the system
  varies linearly between two instants of time

post-correct the drift (called trend, wander, baseline in other domains), with specific techniques based on sparsity of derivatives (BEADS: baseline placement) or using piecewise linear models like  LOWESS:
.


As you said, a linear combination of two signals $x_1[n]$ and $x_2[n]$ is simply

$$y[n]=w_1x_1[n]+w_2x_2[n]$$

with real-valued weights $w_1$ and $w_2$.


What information does the abs(normalized_image - normalized_quantised_image) carry?

Assuming your image is a $K$ bits depth image, and the quantized version is a $3$ bits version, pixels of the difference image, represent $(K-3)$ least significant bits of pixels of the original image. It means it completely is similar to original image in darker parts and completely masks the brighter parts.


An image moment is simply a number that characterizes the image, construed as the realization of a spatial random variable. If you've taken any probability class, you should remember the concepts of mean and variance, which are derived from the first and second moments of the random variable (the n'th moment of an r.v. is the expectation of its n'th power). Furthermore, the moments of a random variable collectively prescribe its distribution. In other words, you can reduce a probability distribution to a sequence of numbers, and this is useful when you want to compare distributions numerically.

In image processing, if you want to compare images, you might not want the comparison to be sensitive to minor things like rotation, translation, and scale (since the image remains fundamentally the same). Thus the motivation behind invariant moments you see in the Wikipedia article you cited.


Design a LTI system which returns DC value of the input signal

The DC value is simply the mean. Since the signal is periodic you only need to take the mean of one period. This can be simply done with

$$h(t)=\left\{\begin{matrix}
1/T,  0 < t < T\\ 0, {\rm elsewhere}
\end{matrix}\right. $$

  Clearly (1) and (2) are different.

(1) will not converge, so that doesn't work

(2) is the best way for a function that's not periodic. It will work for a periodic function as well but is needlessly complicated


What is the signal-noise model of laser induced fluorescence?

Second try!  This question is also at the Electrical Engineering stack, so I posted a short 'starter' answer there. A mod converted it to a comment and my tiny rep there prevented me making any further comments. This is what I posted there:

I suggest this link to cbeleites's answer to a question involving laser-excited Raman spectroscopy: https://stats.stackexchange.com/a/194247/247352 and this golden oldie paper: N.J. Dovichi, J.C. Martin, J.H. Jett, M. Trkula, R.A. Keller, "Laser-Induced Fluorescence of Flowing Samples as an Approach to Single-Molecule Detection in Liquids", Anal. Chem.,56 (1984) 348-354. There are thousands of laser-induced (aka, laser excited) fluorescence papers, but newer research papers tend not to dwell on the fundamentals.

My research career, prior to retirement, involved the use and study of laser-based spectrometries, including laser-induced fluorescence (LIF) of molecular and atomic species. So the reason for the above paragraph was that I simply cannot understand how the noise can decrease with increased number of photons. Likewise, the correlation of the signal and noise is puzzling. Hence the references above.

I have been trying to simulate what I think is happening and get this:

Maybe I am on the wrong track, but my model used 1) shot noise on the fluorescence peaks, 2) additive, white gaussian noise on the baseline and 3) some laser intensity fluctuations. It would help to know exactly how the PMT (or photodetector) output current leads to the observed fluorescence signals. Also, are there any unusual time constants that could result in unsuspected pulse stretching? The peaks in the OP's figure presumably have widths determined by the times it takes labelled cells to transit the region where the laser is focussed. This would be on the order of ms, so perhaps provides a clue about the correlation issue.

EDIT in response to OP queries. Without knowing any details, I assumed a fairly typical system: 1) CW laser excitation, 2) fluorophore-labeled cells, 3) a flow cytometer of some sort, 4) microscope objective, 5) PMT or other high gain photodetector, 6) preamplifier, 7) data acquisition, and so forth. As noted above, I assumed shot noise on the fluorescence emissions, white Gaussian noise (from the preamplifier) on the baseline, and about 0.5% Gaussian intensity fluctuations on the excitation laser. I ignored the variation in cell fluorophore loadings, because this has somewhat similar effect as laser intensity fluctuations. The simulation model (program, actually) is shown below:

Every simulation run (start = 0, end = 0.6, $\Delta t = 0.001$) results in a new temporal response plot. The numbers I used were just arrived at by trial and error, trying to get close to the OP's plot. Lowering the gains on the three peaks, via the triangular gain blocks that have 'x 1000' under them, results in lower SNRs. The next two plots show the results for gain = 300 and gain = 100.

I do not know if this is going to be of any help to the OP. Maybe it is over-simplified or I am missing something essential, but I did not have much to go on. I published my first laser induced (excited) molecular fluorescence paper in 1981 and my last in 2014, but I never encountered a molecular LIF system where the noise decreased with increasing number of emitted fluorescence photons. The correlation is also puzzling to me since I cannot think of how it might arise. Obviously, I am very curious to know what is going on!


MATLAB: Generating and Visualizing Mother Wavelet of different scale

Subsequent wavelets, spanning a specific range of scales are generated by shifting and scaling the Mother Wavelet.

Therefore, to produce "...a similar Daubechies with frequency twice larger..." you would have to "take every other sample" of the Mother Wavelet. For more information, please see this link and the Wavelet Toolbox User Guide on the section about the discrete wavelet transform.

It may also be useful to have a quick look at decimation and upsampling concepts in DSP.

Hope this helps.


How is the energy of $x_1\cdot x_2$ related to the energies of $x_1$ and $x_2$?

Knowing the energies of $x_1$ and $x_2$ is not sufficient for determining the energy of $x_3=x_1x_2$. What you can do is determine an upper bound for the energy of $x_3$ given the energies of $x_1$ and $x_2$ and their maximum values:

$$E_3=\sum_{k}\big|x_1[k]x_2[k]\big|^2\le\begin{cases}\max_k\big|x_1[k]\big|^2\sum_k\big|x_2[k]\big|^2=\max_k\big|x_1[k]\big|^2E_2\\\max_k\big|x_2[k]\big|^2\sum_k\big|x_1[k]\big|^2=\max_k\big|x_2[k]\big|^2E_1
\end{cases}\tag{1}$$

From $(1)$ an upper bound for $E_3$ is

$$E_3\le\min\left\{\max_k\big|x_1[k]\big|^2E_2,\max_k\big|x_2[k]\big|^2E_1\right\}\tag{2}$$

Unfortunately, this bound is generally not very tight. For your example you get

$$E_3\le\min\{19.178,17.787\}=17.787$$

but the actual value is $E_3=7.8165$.


Introduction to Clustering and Classification of Linear Sequence

This seems very much like a fit for Hidden Markov Models (HMM). They are pretty powerful in handling sequence based input such as time series.


Can I use standard computer vision techniques for images taken in the NIR spectral range?

As far as I know, it is a common practice to use regular image processing techniques in this field. In fact I have recently finished a project in this field.

  if it is possible to employ standard stereo matching algorithms like
  block-matching to images taken in the NIR spectrum and if so would it
  be possible (of course depending on the conditions) to obtain a dense
  depth map.

Yes, it is, I think, since NIR images as opposed to Far-Infra red Images has a lot of texture which you probability can (It is very very similar to RED band of the image), as an example people doing these:

repository.upenn.edu/cgi/viewcontent.cgi?article=1115&context=cis_reports

  if standard feature detection and description techniques like SIFT, HOG,etc produces similar results in the NIR spectrum.

I think so, since again it has lots of texture and wherever there is texture is corners, edges and all types of  feature points,which you can extract for further porcessing, check this out:

enter link description here

I think you are better to also take a look at reference books such as:

Machine Vision and Advanced Image Processing in Remote Sensing:

https://books.google.com/books?isbn=3642601057

Remote Sensing, Models, and Methods for Image Processing:

https://books.google.com/books?isbn=0126289816

Hope you've got your answer.


Power spectrum of 2D image: result interpretation

it does not make any sense to still use circles for calculating the histogram. I'd have to use ellipses, wouldn't I?

That's right. For radial averaging on a rectangular grid, the locus of equidistant points from the origin would look like an ellipse instead of a circle.

  there are four areas (area of the square minus the area of the circle) in each corner, from which I don't get any contribution

That's right, and there's nothing we can do about that. If you care about those high frequencies, you should re-acquire your image at a higher sampling rate.

  they used a crooked square, i.e. a rhombus, with each corner on one of the axis. ... Does this approach make any sense? What about the four triangle-shaped areas in this case?

In this case the locus represents points where the sum of the two spatial frequencies is a constant, i.e. instead of $\sqrt{f_1^2 + f_2^2} = \mbox{constant}$ we follow the points where $|f_1| + |f_2|=\mbox{constant}$. There's no reason why one is better or worse than the other. They are both valid methods of collapsing a 2D function of two spatial frequencies into a 1D function of one "combined" spatial frequency.

As with the circle/ellipse, there's nothing we can do about the bins outside the triangle area that never get analyzed.

When using a DFT, distance can be measured as the number of steps taken on a discrete 2D frequency grid. (See the discussion at the bottom of the page here: http://iubemcenter.indiana.edu/beyondNyquist.html). Constant frequency "ellipses" in this case are really "rectangles" of grid cells centered around the origin. This way all grid points get accounted for and there are no "dead" zones when computing the 1D collapsed PSD from the 2D DFT.


Does chirp have constant magnitude frequency response?

The book is not wrong, but it does present the concepts on LFM in a clunky manner and can be misleading. The book presents the analytical expression for the LFM spectrum, which is an approximation. It also plays with the plot views and most likely unwraps the phase angles, which is usually required to see the phases you expect.
Usually when you're introducing LFM, you'll show the modulated pulse itself as well as the phase progression in the time-domain. The analytical expressions in the time domain is all you need to observe the linear frequency and therefore quadratic phase nature of LFM. Doing this in the frequency-domain just tends to introduce more confusion. An example of the time-domain LFM pulse and its phase is shown below.

When deriving the expression for the Fourier transform of an LFM pulse, you do indeed yield a magnitude of 1 over the bandwidth of the pulse. This is intuitively satisfying because you have the same contribution from each frequency over the bandwidth.
Confusion does occur however when one goes to plot this if they expect a constant frequency response. With any practical LFM spectrum plot, even with very long pulse widths, one should expect a ripple effect which you have already identified. The quadratic nature of the phase is still captured in the DFT. The spectrum of the LFM pulse from above is shown below.

I haven't tried your code, but it might be that you just need to zoom in the proper areas and unwrap the phases to see what you want. Provided is the MATLAB code to produce the plots to help you in converting it to python.
%% LFM - Time and Frequency Domain

% Sampling
Fs = 50e6;

% Pulse parameters
tau = 50e-6;
bandWidth = 10e6;
alpha = bandWidth/tau;

% Define waveform
t = 0:1/Fs:tau - 1/Fs;
fmcwPulse = exp(1i*pi*alpha.*t.^2); % Complex transmitted LFM waveform

% Plot
figure;
subplot(2, 1, 1);
plot(real(fmcwPulse));
xlabel("Samples");
title("LFM Pulse - Real Part");

subplot(2, 1, 2);
plot(unwrap(angle(fmcwPulse)));
xlabel("Samples");
title("LFM Phase");

figure;
subplot(2, 1, 1);
plot(abs(fftshift(fft(fmcwPulse))));
title("LFM Spectrum Magnitude");

subplot(2, 1, 2);
plot(unwrap(angle(fftshift(fft(fmcwPulse)))));
title("LFM Spectrum Phase");

Update
Modifying the code above so that $\tau = 1 \space s$, which is relatively long, yields a spectrum closer to the ideal flat spectrum that one would expect analytically. The spectrum is shown below.


Cross correlation with short or non periodic signals?

Given two signals with the same finite duration and number of samples $N$, we can calculate the periodic  or cyclic crosscorrelation function which has $N$ values, or the aperiodic crosscorrelation function which has $2N-1$ values. For signals with $N$ and $M$ samples respectively ($N > M$), the aperiodic crosscorrelation function has $N+M-1$ values.

When the two signals are the same except that one is a cyclic shift of the other, the periodic crosscorrelation function has a maximum at perfect alignment. However, the maximum of the aperiodic crosscorrelation function is not necessarily at perfect alignment because the part that is "sticking out" and not overlapping with the other signal is missing, which lowers the peak.

When one signal is a truncated version of the other, and the longer signal is periodic, then the aperiodic crosscorrelation function will have multiple maxima all of the same height. Look at the array that
xcorr gives you as th eoutput, and check how the max function works in MATLAB when there are multiple maxima. Does max report the leftmost
maximum if there are multiple maxima of the same height?  (I am not a
MATLABi and cannot help you on this last point).


Ways to differentiate three respiratory signals

You're talking about ECG-derived respiratory (EDR) signals, but you didn't detail how you obtained them. Please detail your methods. Maybe then it will be easier to understand what each of your signals mean.

By 'experimental signal', I take it that you mean the reference/gold standard annotation of when respiration occurs, am I correct? But you did not detail if this was the onset or end or peak of the breathing cycle, as well as what kind of measure it was (tidal volume, airflow, other measure?)

Edit #1:
As for features, frequency is a good starting point. EDR signals may not evolve the same way as the reference signals but they should evolve at about the same period. If the breathing is controlled during a significant amount of samples, you can compare the frequency of the whole signals. If not controlled, then take segments of more than 2 seconds. You can apply FFT to get the frequency.

Phase will also be of interest, but it will difficult to comment further without knowing the algorithm for calculation of the EDR as well as the reference signal.


Derivation of fixed-point $\tt atan2$ with self-normalization

it's an approximation.

a better approximation is in this answer:

$$ \arctan(u) \approx \frac{u}{f(u^2)}  \quad \quad -1 \le u \le 1  $$

where

$$ \begin{align}
f(u^2)& = \sum\limits_{n=0}^{4} \ a_n \ u^{2n} \\
a_0  &  =  1.0 \\
a_1  &  =  0.33288950512027 \\
a_2  &  =  -0.08467922817644 \\
a_3  &  =  0.03252232640125 \\
a_4  &  =  -0.00749305860992 \\
\end{align} $$

not as fast, but quite accurate.  not bit-perfect accurate, but very accurate.

for atan2(), you need to put it in the correct tilted quadrant:

$$ \operatorname{atan2}(y,x) = \arg\{ x + j\,y \} =
\begin{cases}
\arctan\left(\frac{y}{x}\right) &\text{if } 0 < |y| \le x, \\
\frac{\pi}{2} - \arctan\left(\frac{x}{y}\right) &\text{if } 0 < |x| \le y, \\
-\frac{\pi}{2} - \arctan\left(\frac{x}{y}\right) &\text{if } y \le -|x| < 0, \\
\arctan\left(\frac{y}{x}\right) \pm \pi &\text{if } x \le -|y| < 0, \\
\text{undefined} &\text{if } x = 0 \text{ and } y = 0.
\end{cases} $$

this will require two fixed-point divisions, but you can arrange it so that the result is always no greater than 1 in magnitude (the numerator is always smaller in magnitude than the denominator).


z-transform causality properties: negative coefficents are zero ($x[-1]z^1=0$, $x[-2]z^2=0$, ...)

No. The system response is usually $H(z)$. If $H$ is causal, then it has no non-causal coefficients.  $X$ may very well be a signal that exists for all time (and so existed before time 0) and, therefore, have non-causal coefficients.  That means $Y$ will also have non-causal coefficients.

If $Y$ is taken as the system and $Y$ is causal, then your statement is true. But usually $Y$ is the output, $X$ is the input, and $H$ is the system response.


how much of DSP (digital signal processing) is related VLSI?

First, an example: designing FIR filters pops to mind, or for that matter any signal processing module. The FIR filter is a nice example, since it can easily be implemented in VLSI even at the lowest abstraction level - physical mask level (in Cadence for example).

If we recall an example of the FIR filter design:

We can see that we only need to design:

Delay
Multiplier
Adder

elements, which are quite easy to implement in VLSI, either at the high level (VHDL, Verilog), or lower level (Cadence, Magic IC, etc.).

So now we have a mixed problem - the DSP part is designing the FIR filter, which can be designed in various ways, and dealing with the transfer function/frequency response. The VLSI part is actually implementing the filter - creating the layout which can be printed on a real integrated circuit.

In general, the two fields are connected, but in industry usually one deals in only a single layer of abstraction. Either one deals with circuits/elements as black boxes, without caring so much about the physical implementation, or conversely focuses on the optimal design of a component, which would later be produced and used in projects where DSP is employed.


How to calculate the displacement from the measured acceleration data?

Generally, the observed effect occurs when there is an offset and/or linear drift present in your measured data. Double integration then leads to the quadratic (or higher order) effect observed in the data.

As @Stanley Pawlukiewicz states in his answer this can be due to the accelerometer registering gravity (not all accelerometers are capable of registering a constant acceleration field) but also due to temperature fluctuations, (electro)magnetic interference, strain at the base of the accelerometer, etc. Sometimes, it’s possible to “shield” your sensor for these effects but this can also be extremely challenging or expensive. In the latter case pre-processing your data can be an option too enable extracting meaningful results from your data.

The most straightforward approach to get rid of the quadratic drift in your double integrated acceleration signal is to detrend the data before integrating. If you substract the mean value before integrating the acceleration and velocity signals you should get reasonably “close” to the desired result (see matlab’s detrend function for other options).

Note that there also appears to be other noise and/or content in your measured acceleration signal; i.e. you will probably be able to get “close” to the desired result but will not get the exact same position signal. If that is critical you could use a more sophisticated signal processing approach or, if possible (and preferably) a better sensor (e.g. a position sensor).


What is the physical interpretation of the absolute value of a fourier transformed signal, $\left| F(t)\right|$?

This would be the so called magnitude spectrum, its unit is the same as the signal it's generated from, in the case of voltage its interpretation is "voltage  by time over frequency" as opposed to "power by time over frequency" of the squared magnitude spectrum.
In the ususal logarithmic representation, the two are proportional by a factor of 2 and essentially give you the same information.

The $\cos$ function in your example is a so called power signal, not an energy signal, as it is infinitly long it has infinite energy, so to assess its energy in a practical way, one refers to its power (energy per time).


I think the  most fundamental things that anybody trying to understand DSP should know are:

Nyqvist Sampling Theorem and Effects of Aliasing
Need for DFT and the use of FFT to achieve the same.
Basics of Digital Filters.

You have mentioned the audience as telecommunication engineers. So you could also add a lecture on the  need for modulation,some basic modulation and demodulation techniques and need for interpolation and decimation.


Octave error - STFT parameters when detecting pitch in guitar audio signals

The octave error is because you are using a spectral frequency estimator (STFT), not a psychoacoustic pitch detection or estimation method.  The spectral frequency peaks in the sound are often not the pitch frequency, but instead often one of the many harmonics.  As you found, the fundamental pitch frequency can even be completely missing.  This is especially true for guitar notes, voices, and sounds from other large stringed instruments.

The solution is not to play with sample rate and STFT FFT window parameters, but to look up more reliable pitch detection and estimation algorithms.


FFT "Brick Wall" Low Pass Filter Not Smoothing Output

Sorry for the post I just solved it.

The problem was I was using a full transform. The half-transform utilizes a symmetry condition and for some reason works infinitely better!


One approach would be to use the spectral flux of a sequence of STFT frames $X(n,k)$, namely

$$
\mathit{SF}(n) = \sum_{k = -N/2}^{N/2 - 1}{H(\left|X(n,k)\right| - \left|X(n-1,k)\right|)},
$$

where $n$ is the frame number, $N$ is the size (in samples) of each frame, and $H(x) = \frac{x + \left|x\right|}{2}$ is the half-wave rectifier function. Intuitively, it attempts to measure how much the spectrum changes from frame to frame.

The overall area of research is called onset detection. A good overview of techniques can be found in the following paper (which also describes spectral flux):

  S. Dixon. Onset detection revisited. In Proceedings of the 9th
  International Conference on Digital Audio Effects, volume 120, pages
  133-137, 2006.

An example implementation of spectral flux can be found in the NSGToolbox.


Rayleigh energy equation - applied to Fourier Transform of a sinusoid

An (infinite) sinusoidal signal does not really have a Fourier transform, since as you say it is not a finite-energy signal. It can be represented as a Fourier series, and there is a useful way to sort of write it as a "pseudo-Fourier transform" with Diract delta distributions, as well. However, if you try to take the magnitude of that spectrum and square it, it is not well-defined, since you would need the squared-magnitude of the Dirac deltas.


Time reversal of input signals in convolution

This question appears to involve confusion about the relationship between convolution and cross-correlation. To illustrate the relationship, consider the two rectangular pulses shown below:

The shorter pulse has width = 1 and the longer pulse has width = 2. For simplicity, assume unity heights for the pulses.
In my answer to a previous question, I used a figure from a textbook 1 to illustrate the convolution of $\mathrm {x_1(\lambda)}$ and $\mathrm {x_2(\lambda)}$, denoted by $\mathrm {x_1*x_2}$.
The convolution integral was given by
$$\mathrm {x_1*x_2 =\int_{-\infty}^{\infty}x_1(\lambda)x_2(t-\lambda)\,\mathrm d\lambda} \tag{1}$$
so $\mathrm {x_2(\lambda)}$ was reflected around the ordinate and is the function that gets shifted. This is clearly shown in the figure in my linked answer.
But it is also possible to reflect $\mathrm {x_1(\lambda)}$ about the ordinate and shift it. Then the convolution would be denoted by $\mathrm {x_2*x_1}$ and the convolution integral would be
$$\mathrm {x_2*x_1 =\int_{-\infty}^{\infty}x_1(t-\lambda)x_2(\lambda)\,\mathrm d\lambda} \tag{2}$$
This is illustrated the following figure that I drew in the same fashion as the figure in my linked answer:

Thus $\mathrm {x_1*x_2 = x_2*x_1}$, i.e., it does not matter which function is reversed for the convolution. Of course, if a function is symmetric about its center, then reversing becomes rather trivial. In this example, both rectangular pulses are symmetric about their centers, so the situation is especially simple.
In corroboration of the above, the following (through their equation 1.5-3) is directly quoted from Blinchikoff and Zverev [2, p. 12]:

Let p(t) and q(t) each be a piecewise continuous function for $-\infty < t < \infty$. The convolution of p and q, denoted by $p*q$, is defined to be the third function
$$\mathrm {v(t) = p*q =\int_{-\infty}^{\infty}p(t-\tau)g(\tau)\,\mathrm d\tau} \tag{1.5-1}$$
and the integral on the right is called the convolution integral. Furthermore, the convolution of q and p, $q*p$, is identical to p*q, that is,
$$q*p = p*q \tag{1.5-2}$$
This is easily shown by letting $\mathrm {x = t -\tau}$, $\mathrm {dx = -d\tau}$ in (1.5-1), yielding
$$\mathrm {p*q =\int_{-\infty}^{\infty}p(x)g(t-x)\,\mathrm dx} = q*p \tag{1.5-3}$$

Therefore, in answer to the OP's second question, it does not matter mathematically if f(t) or g(t) is the reversed function, but the calculation may be simplified, in practice, if either function is symmetric about its center.
For cross-correlation, all that is necessary is to shift one function by $\tau$, as per this photo image from [2, p. 305]:

So cross-correlation calculations are quite similar to convolution calculations, but note the important considerations discussed directly below Blinchikoff and Zverev's Table 7-2. I suggest that the OP think about their first question, i.e., what happens when both f(t) and g(t) are reversed and the overlap integrals are calculated.
1  C.D. McGillem, G.R. Cooper, "Continuous and Discrete Signal and System Analysis", 2nd Ed., Holt, Rinehart and Winston, NY, ©1984, pp. 58-59.
2  H.J. Blinchikoff, A.I. Zverev, "Filtering in the Time and Frequency Domains", Wiley-Interscience, John Wiley & Sons, NY, ©1976.


How is it determined that an adaptive filter has converged?

The adaptive filter is converged when the error is what they call Wide Sense Stationary, meaning the mean and variance of the error are unchanging over long time intervals.


Downsampling and low pass filtering in one step?

Is it possible to combine decimation and low pass filtering in one step? Not necessarily only for images but also for general signals.

Yes, that's what people usually do when they implement downsampling: since of the output of the anti-aliasing filter, you throw away N-1 samples, why even calculate these?
The trick is to decompose your filter into polyphase components, which enables you to run the resulting filter operation only once per output of the downsampling, instead of once per input. There's plenty of reference implementations - from GNU Radio's decimating FIR filters, to rescalers in image processing hardware.
Think of it this way:
The trick is to take your original filter $[h_0, h_1, h_2, h_3, \ldots, h_N, h_{N+1}, h_{N+2},\ldots,h_{2N}, h_{2N+1}, \ldots]$ and just split it up into filters where there's only one non-zero entry every $N$ coefficients. Choose the non-zero-value positions so that the first polyphase component filter gets $h_0, h_N, h_{2N},\ldots $, the second gets $h_1, h_{N+1}, h_{2N+1},\ldots$ and so on.
Add up the result of these filters, when you feed in the same input, to "undo" the splitting. This doesn't change anything, it's the same filter, just split into $N$ filters with many zeros in them, but with the non-zero elements in different positions.
After the addition, you decimate by $N$. Ok, you can do that before the addition, so now you have one input stream, fed into $N$ subfilters, each with a lot of zeros in them, each followed by a decimation by $N$.
Now you have a special kind of filter that only had every Nth filter tap occupied, so the first subfilter's coefficient vector is $[h_0, 0, \ldots, 0, h_N, 0, \ldots, 0, h_{2N}, 0 \ldots]$, and you'd decimate by $N$ afterwards, you could as well just swap decimation and filter, and just use the filter $[h_0,h_N,h_{2N},\ldots]$. The two things are identical in effect; this is called Noble Identity.
So, we can "pull the decimation up front" for that filter. You can, in fact, do that to all subfilters (you will have to add delay so that it works out mathematically for the non-zero-phase polyphase components, but the idea doesn't change. You have one input stream, going into $N$ different delays, decimate-by-$N$ decimators, subfilters, and a summation.
As it happens, this means that only one "branch" at a time actually gets input per input cycle.


The result you experience is reasonable. The IDFT performs exactly the inverse operation of the forward DFT. Hence, for every signal $x[n], n=0...L$ you have the identity:

$$x[n] = \text{IDFT}\{\text{DFT}\{x[n]\}\}$$

This property is valid for DFT and hence also for FFT (as it is a special implementation of DFT).

So, if you insert a zero-padded sequence $x[n]$ into the FFT and immediately perform IFFT you obtain the zeropadded sequence as the output.

The question is: Why would you want to do that? If you think of performing filtering in the frequency domain, dont do it.


Is it possible to design a digital filter that rejects 60Hz noise while keeping step response in the time domain below one mili-second

You can trade off steepness of the notch against time domain ringing. That's a fundamental trade off and there is no way around it. When using an IIR notch you can simply adjust the "Q" (quality factor). A Q of 3 will give a step response that's roughly 10ms long.


Calculating Savitsky-Golay filter coefficients - normalization

The SciPy library for scientific computing in Python contains functions for Savitzky-Golay filtering in its scipy.signal module since version 0.14.0, specifically you'll most likely be interested in the scipy.signal.savgol_coeffs function if you're only interested in the coefficients, or the scipy.signal.savgol_filter function which provides a nice interface for filtering alone.

When implementing some DSP method in Python, it's usually best to look into SciPy or scikits-* libraries rather than reinvent the wheel.


Selectively improve the frequency resolution of FFT (DFT)

You can compute any single FFT result bin using a complex Goertzel  algorithm or filter.

Realize however that if you want to compute a number of bins on the order of log N or greater (such as 10 extra bins with 10 Hz bandwidth and spacing), it may be faster to do the full longer FFT and throw away all the extra FFT result bins.

So just compute 2 FFTs of 2 different lengths and ignore all the unwanted results.


Showing linear-phase frequency response when only real part is known

You have two kinds information in the question:

The system is causal
The real part of the Frequency response is given.

Now, (repeating the steps from Wikipedia Entry on Causal Filter)

$h(t)$ is causal, i.e. $h(t)=0, t<0$.
Let $g(t)=\frac{1}{2} (h(t)+h^*(-t))$ which is non-causal, but has hermitian symmetry, hence its Fourier Transform is real.
We know (due to causality) $h(t)=2u(t)g(t)$ where $u(t)$ is the unit-step function
Now, computing Fourier-Transform we get
$$
H(\omega) = 2 \mathcal{F}(u(t))*G(\omega)
$$
where $*$ is the convolution operator.
Now, we can express this as
$$
H(\omega) = \left(\delta(\omega)-\frac{j}{\pi\omega}\right)*G(\omega)=G(\omega)-j\mathcal{H}(G(\omega))
$$
where $\mathcal{H}(G(\omega))$ is the Hilbert Transform of $G(\omega)$. Since we know that $G(\omega)$ is real, $G(\omega)$ equals the real part of $H(\omega)$. Then, $\mathcal{H}(G(\omega))$ is the imaginary part of $H(\omega)$.

So, to get the imaginary part, you need to calculate the Hilbert Transform of $1+\cos(\omega)$ which is given by

$$
\mathcal{H}(1+\cos(\omega))=\sin(\omega).
$$

Accordingly, the phase of your system is given by

$$
\phi(\omega) = \tan^{-1}\left(\frac{-\sin(\omega)}{1+\cos(\omega)}\right)
$$

Now, if we calculate the derivative regarding omega, we need to get a constant for a linear phase filter:

$$\frac{\partial \tan ^{-1}\left(\frac{-\sin (x)}{\cos (x)+1}\right)}{\partial x}=-\frac{\frac{\cos (x)}{\cos (x)+1}+\frac{\sin ^2(x)}{(\cos (x)+1)^2}}{\frac{\sin ^2(x)}{(\cos (x)+1)^2}+1}=-\frac{1}{2}$$

So, this shows that you have indeed a linear-phase system.

Your way with the group delay is also correct, once you have the imaginary part of your $H(\omega)$ obtained from its Hilbert Transform.


Can I segment characters from this image with CSV?

Your results (above) looks well for character segmentation. But if you have some problems with calculation of thresholds for total image, you can try slightly different approach. "Column sum vector" method works well for grey images too, but not only for black/white images. So there is no need to calculate threshold for total image, but you must choose some value for segmentation in vector. May be it will slightly simpler to do segmentation in this case.

You can use "line sum vector" for vertical segmentation   too.


Re question (1):

Most FFTs, or the way atan2() is commonly used with complex FFT results to compute a phase, produce a zero phase (strictly real) output for strictly symmetric input.  This means you want cosine waves, periodic in aperture, not sine waves starting at zero at the window edge, for FFT results with zero phase.


why don't we use Broadband technology to transmit digital signals ,and use Baseband technology to transmit analog signals?

"Broadband technology" is not inherently coupled with "analog signal transmission", as "baseband technology" is not coupled with "digital signal transmission". Both wire cables and optical fiber cables transmit signals in the form of electromagnetic waves, in RF frequency range for a wire medium and in IR/light frequency range for an optical fiber medium.
For digital transmission, the source data are transformed to changes in amplitude/frequency/phase parameters of a generated output signal and then the pulse shaping and filtering operations are performed on this signal. Resulting voltage/current or amplitude/phase waveforms have no sharp rise/fall jumps like those you see on the logic analyzer display. In this way, the resulting waveforms look similar to an analog signal.
In a baseband transmission the signal is sent "as is"; in a broadband transmission a modulator device creates a passband signal from a baseband signal and a carrier. The passband signal frequency is shifted by modulation process and occupies the frequency band around the carrier frequency with the width of the baseband signal spectrum. If the entire channel is used for a baseband transmission, we can use a medium's full bandwidth for the baseband signal. On the other hand, passband techniques allow frequency-division multiplexing, and we can share the channel for multiple concurrent transmissions. The number of subchannels is defined by a medium's full bandwidth divided by a baseband signal spectrum width.
Cable television systems use "broadband technology" for distribution of multiple analog video channels. Digital networks also can use frequency-division multiplexing (the "broadband technology", as your reference calls it). In 1980s, 10BROAD36 Ethernet standard was developed in an attempt to use analog TV cables for Ethernet networking. The 10BROAD36 communication line was even able to share a cable with TV channels. The attempt failed for a number of reasons. What's more, the frequency multiplexing is not a single (and so unavoidable) solution for sharing communication lines in packet-switching networks of Internet. The sharing of a communication line to multiple data streams is performed with packet switching. Still, frequency-division multiplexing is used in wired (Discrete Multitone Transmission) and wireless (WiFi, OFDM) digital transmissions.
As for analog transmissions in rare survived analog telecom systems, I do not know of analog signals with spectrum that wide as to occupy the entire bandwith of transmission lines. So these systems can safely use frequency multiplexing to share cables for multiple transmissions ("broadband technology", as your reference calls it).


Delay-and-sum beamforming: implementing steering in C

The code is calculating the beam pattern of a linear array which is the response at a single frequency.You might note that no actual sensor data is being summed or delayed.

The pattern is defined over all the angles from -90 to 90 and will have a main lobe.

There are essentially two ways of steering a beam, by multiplication by phase at a single frequency or by actually delaying the signal in the time domain, typically using a N point interpolation, where N is small.

For 4 elements,  500 hundred look directions is rather too much.


Correlation : How is correlation analogous to dot product?

So using your formula for the cross-correlation:

$$R_{xy}(n)=\sum_m x(m)y(n+m)$$

we see that

$$\begin{align}
R_{xy}[n]&=\sum_m x[m]y[n+m]\\
&= \sum_m x[m] \tilde y_n[m]\\
&= x \cdot \tilde y_n\\
\text{for}\\
\tilde y_n[k] &= y[n+k]\text{ ,}\\
\end{align}$$
which is a dot product.

So for any "shift" in the cross-correlation, you have the dot product of the first operand with the $n$-shifted second operand.

Note, by the way, that your formulas are only correct for real-valued $x$, i.e. in general, the correct formula for dot product and cross correlation are:

$$\begin{align}
x\cdot y &= \sum\limits_{i=0}^{N-1} x^*[i]y[i]\\
R_{xy}[n] &= \sum\limits_{i=0}^{N-1} x^*[i]y[n+i]\\
\text{with}\\
x,y &\in \mathbb C^N\\
{z}^*&= \Re\{z\} - \Im\{z\}
\end{align}$$


Generic response of an IIR filter from its transfer function


Noise Equivalent Bandwith of IIR and FIR Filters?

Hm, isn't it the same as for any filter?

$$
B_\text{eq} = \frac{1}{|H(e^{j\omega_0})|^2} \int_0^\infty |H(e^{j\omega})|^2 d\omega
$$

where $H(e^{j\omega})$ is the frequency response, and $\omega_0$ is the max-abs frequency. For FIR filters, you can just use e.g. freqz in MATLAB/Octave. For IIR filters, you'd have to do analysis or take measurements (or see if you can get away with truncating)


FFT for distance-limited sum of pairwise products

Your sequences are all correlations between $x_i$ and $x_j$, where you assume that $x_i[n]=0$ for $n<0$ or $n\geq N$. Correlation is very related to convolution and it has a similar expression in frequency domain:

$$ c[n]=\sum_l x_1[l]x_2[l+n] \Leftrightarrow C[k]=X_1^*[k]X_2[k]$$

(capital letters are Fourier transforms of the small ones).

Hence, you can efficiently calculate the correlation using this theorem.

One caveat is that you are actually only looking for linear correlation, but using DFT you calculate circular correlation. So, the solution is to zero-pad the sequences. Have a look at the following Python code:

N = 5
x1 = np.random.randn(N)
x2 = np.random.randn(N)

x1z = np.hstack([x1, np.zeros(N)])
x2z = np.hstack([x2, np.zeros(N)])

# calculate according to your expression
d12 = []
for l in range(N):
    d12.append(sum(x1[i]*x2[i+l] for i in range(N-l)))

# calculate using zero-padding of the sequences with no
# dependence of the range of i
d12z = []
for l in range(N):
    d12z.append(sum(x1z[i]*x2z[i+l] for i in range(N)))

# calculate using circular correlation, using the zero-padded versions
d12c = []
for l in range(N):
    d12c.append(sum(x1z[i]*x2z[(i+l)%len(x2z)] for i in range(len(x1z))))
# check they all match
print ("x1:", x1)
print ("x2:", x2)
print ("d12_1:", d12)
print ("d12_2:", d12z)
print ("d12_3:", d12c)

# now, exploit the correlation theorem to calculate the same stuff:
C = (np.fft.ifft(np.fft.fft(x1z)*np.fft.fft(x2z).conj()).real)
# take only the first N values, as the others correspond to parts
# where the sequences do not fully overlap.
print ("FFT:", C[-np.arange(N)])

A sample output could then be:

x1: [ 1.11998979 -0.81341718  1.47117346  0.59699824 -0.9233909 ]
x2: [-0.31884542 -0.2214878  -0.02678967 -0.25943789  1.36355875]
d12_1: [-1.6303355800547898, 0.20609114482882795, 2.1870585244152387, -1.3997099032157734, 1.5271718773447336]
d12_2: [-1.6303355800547898, 0.20609114482882795, 2.1870585244152387, -1.3997099032157734, 1.5271718773447336]
d12_3: [-1.6303355800547898, 0.20609114482882795, 2.1870585244152387, -1.3997099032157734, 1.5271718773447336]
FFT: [-1.63033558  0.20609114  2.18705852 -1.3997099   1.52717188]


Since you are using Matlab, I suggest using SPGL1 as your reconstruction solver. It can handle complex numbers and you can specify your matrix operator as a function.

This example demonstrates how to specify your matrix operator as a function instead of a matrix (see opA).

In your case, I would define your operator as a function that:

takes a vector x as input,
reshapes that vector to a matrix X of your desired shape,
applies fft to the matrix: Z = fft(X) (Matlab's fft transforms each column separately if you give it a matrix),
applies fft to the transpose of the result: W = fft(Z'),
reshapes W back into a vector w and returns that

If your operation additionally includes some sub-sampling (which it usually does in order to be compressed sensing) then that can be implemented as indexing into w for the entries that are sampled.

The above steps imply that x is in the domain of the solution you are seeking and therefore sparse. Thus, the sparse vector here is in the time domain and you are sampling in the Fourier domain. If, on the other hand you want to model a signal that is sparse in frequency and you sample in the time domain, you should replace fft by ifft above.


Usually, such figures show the error rate on a log scale, and $E_b/N_0$ in dB. So you would specify $E_b/N_0$ in dB, then convert it to a ratio, and use that value to evaluate the formula for the bit error rate (as you did). Then you plot the logarithm of the error rate as a function of $E_b/N_0$ in dB. So the main problem is that you used a linear scale for the BER instead of a logarithmic scale.

If you do things right, you'll end up with something like this:


Is there an equivalent of Parseval's theorem for wavelets?

Yes indeed! In theory as long as the wavelet is orthogonal, the sum of the squares of all the coefficients should be equal to the energy of the signal. In practice, one should be careful that:

the decomposition is not "expansive", i.e. the number of samples and of coefficients is the same.
wavelet filter coefficients are not re-scaled, as happens in some applications (like lifting wavelets, to keep integer computations).
wavelets are orthogonal (this is not the case in JPEG2000 compression).

You can verify this indirectly, looking at approximation coefficients. At each level, their number of samples is halved, and their amplitudes have around a $1.4$ scale factor, which is just $\sqrt{2}$. This feature is used for instance to estimate the Gaussian noise power from wavelet coefficients:

$$ \hat{\sigma} = \textrm{median} (w_i)/0.6745$$

A little further, there is a notion that generalizes (orthonormal) bases: frames. A set of functions $(\phi_i)_{i\in \mathcal{I}}$ ($\mathcal{I}$ is a finite or infinite index set) is a frame if for  all vectors $x$:

$$ C_\flat\|x\|^2 \le \sum_{i\in \mathcal{I}} |<x,\phi_i>|^2\le  C_\sharp\|x\|^2$$

with $0<C_\flat,C_\sharp < \infty$. This is a more general Parseval-Plancherel-like result used for general wavelets.

In other words, it "approximately preserves energy" by projection (inner product). If the constants $ C_\flat$ and  $C_\sharp $ are equal, the frame is said to be tight. Orthonormal bases are non-redundant sets of vectors with $ C_\flat=C_\sharp = 1 $.

For those using Matlab, you should care about the native border extension, which is obtained by dwtmode('status'). Some add tails to the data to help inversion with little border artifacts. With a periodic mode dwtmode('per') and a number of samples that can be divided by $2^L$ where $L$ is the wavelet level, you can get a good match in energy, with tiny differences:


How to Validate Bilateral Filter Implementation?

I'm assuming implementation which is defined by two parameters - Color Range STD and Spatial STD.

Few Properties:

If you set the Color Range STD to be high, the result should be very very similar to Gaussian Blur with the Spatial STD.
If you set the Color Range STD to be very near zero, the output image should be with no change (No smoothing at all).
Good image to test things on is Black White rectangles (Like Checkers).


How to generate random samples of Gaussian distribution directly in the frequency domain?

You can, but... you'll need to keep symmetry if your original time-domain signal is real-valued.

If a signal $x$ is real-valued, then its DFT $X$ will exhibit complex-conjugate symmetry:
$$
X[k] = X^*[N-k].
$$

So you can generate $N$ Gaussian pseudo-random noise samples, $g[n]$, and place them in the frequency domain noise vector, $\epsilon$ as:
$$
\epsilon[k] =  g[k] + j g[k+N/2]
$$
for  $k \in \{ 0, 1, \ldots, N/2-1\}$ and
$$
\epsilon[k] = \epsilon^*[N-k]
$$
for $k \in \{ N/2, N/2+1, \ldots, N-1 \}$ where $\epsilon^*$ is the complex conjugate of $\epsilon$ and is equal to $\Re[{\epsilon}] - \Im[{\epsilon}]$ (i.e. the same real part and the negative of the imaginary part).


Authoritative DSP book/reference material on Psychoacoustics

The Fastl-Zwicker book is usually the introduction to this topic.  The more I learn about it the more fascinating it becomes.  If you pick it up used be sure to get one with the CD examples in it, as they are completely fascinating.


Are 2nd Order Edge Detectors More Susceptible to Noise?

If you assume the Edge Detection is SNR driven operation, one could find a Mathematical justification for this.

First, the variance of Additive White Noise with Variance $ {\sigma}_{n}^{2} $ at the output of a Linear System given by $ g $ is $ {\sigma}_{n}^{2} {\left\| g \right\|}_{2}^{2} $.

Let's look on an Image filtered by a derivative approximating filter $ g = \left[ -1, 1 \right] $.
This increase the Noise Variance by a factor of 2.
For 2nd order the filter is $ {g}_{2} = g \ast g = \left[ 1, -2, 1 \right] $ hence the factor is 6 (It comes from the reason the second time we filter the noise is correlated).

Now, the question does the signal (Edge) is amplified by the same factor?
Moreover, on areas where there is no edge, the probability to detect noise as edge increases.
To see intuitively, just operate 1st order method on the gradient of the image and see how much False edges you have now.


How to measure noise, consonance, dissonance of a sound?

If there are spectral peaks and you want to calculate the dissonance of this sound, you can calculate the sensory dissonance. See the Essentia explanation and the references mentioned there.

If you rather want to measure how noisy a sound is, spectral flatness may be the audio feature of your choice. E.g. see this Wikipedia page.


The Remez exchange method isn't limited to DFT bin constraints so to answer your question, no FIR filters do not have the limitation you hypothesized.

If you can find METEOR, a program written in part by Tom Parks, you can be fairly free in your filter specifications provided they are feesable

and there is

http://gmeteor.sourceforge.net/


Finding the inverse Inverse $\mathcal Z$-transform

Multiply the top and bottom of the fraction by $z^{-60k}$ in each term, and you are left with each term being:

$$ z^{-60k}\frac{1}{z-a} = z^{-60k-1}\frac{1}{1-az^{-1}}
$$

The first factor delays the output by $60k +1$ samples, and the remaining factor is a simple one pole transfer function, giving:

$$  x[n] = \sum_{k=1}^{P} a^{n - 60k - 1} u[n - 60k - 1]
$$


Just so we're on the same page: given some sequence $x[n]$ for $n=0, \dots, N-1$, the $k$'th DFT coefficient is (forgetting the scaling factor)

$$X[k] = \sum_{n=0}^{N-1} x[n] \exp(-j 2\pi k n/N)$$

We can get all the coefficients by writing this as $ X = Dx$, where $D$ is the $N \times N$ DFT matrix and $D_{kn} = e^{-j 2\pi k n /N}$. For the N = 2 case,

$$D = \begin{pmatrix}
e^{-j 2\pi (0)(0) /2} &  e^{-j 2\pi (0)(1)/2}\\
e^{-j 2\pi (1)(0) /2} &  e^{-j 2\pi (1)(1)/2}
\end{pmatrix} =
\begin{pmatrix}
e^{-j 0} &  e^{-j 0}\\
e^{-j 0} &  e^{-j \pi}
\end{pmatrix} = \begin{pmatrix}
1 &  1\\
1 &  -1
\end{pmatrix},$$
for which $D^* D = D^T D = D^2 = 2I$, so (minus the scaling factor), unitary.


Color Artifacts in Fourier Transformed Image

Some things you should keep in mind:

Your image is real, keep your DFT Symmetric.
Numerical issues might cause that after the Inverse DFT imaginary parts are lefts.
Use either abs or real to get rid of them.
One should understand what it means to "Zero" out by mask DFT values.
It is not the "Ideal LPF". It causes Ringing / Halos due to Gibbs Phenomenon. One should use better designed LPF.

I hope it helps.


What is the difference between an impulse response and a kernel?

The impulse response of an operation is the same as the kernel of that operation, but only if the operation can be represented by a kernel (the operation is linear and translation invariant).

For example, consider the following input (the impulse):

  [[0., 0., 0., 0., 0.],
   [0., 0., 0., 0., 0.],
   [0., 0., 1., 0., 0.],
   [0., 0., 0., 0., 0.],
   [0., 0., 0., 0., 0.]]

If you give it as input to a laplace filter, you get:

  [[ 0.,  0.,  0.,  0.,  0.],
   [ 0.,  0.,  1.,  0.,  0.],
   [ 0.,  1., -4.,  1.,  0.],
   [ 0.,  0.,  1.,  0.,  0.],
   [ 0.,  0.,  0.,  0.,  0.]]

This corresponds to the kernel of the laplace filter.

But if you give the same input to a 90 degrees rotation around the centre, you get:

  [[0., 0., 0., 0., 0.],
   [0., 0., 0., 0., 0.],
   [0., 0., 1., 0., 0.],
   [0., 0., 0., 0., 0.],
   [0., 0., 0., 0., 0.]]

(the same as the input)

Obviously, this is not a 'rotation kernel', it is just the identity. This happens because there is no such thing as a rotation kernel, because rotation operations are not translation independent (besides the trivial case of rotating by 0 degrees).

If the operation was not linear, the kernel would not match the impulse response either (for example, for a thresholding operation).


Autocorrelation function and correlation integral

The definition of auto-correlation depends on the type of signal. For random processes, the auto-correlation function is defined by the expectation given in Eq. $(1)$ of your question.

For deterministic signals, there are two definitions, depending on whether the signal is an energy signal (i.e., has finite energy), or a power signal (i.e., has finite power but infinite energy). In the first case, the auto-correlation is defined by Eq. $(2)$ in your question. In the latter case, the integral $(2)$ doesn't exist, and the auto-correlation is defined by the limit given in Eq. $(3)$.


How to connect $\omega$ in time domain to physical frequency?

There are several "frequency" variables in signal processing, the most common all being related by constant factors. I will use the notation from Digital Signal Processing by Proakis & Manolakis:

$F$: "Ordinary" frequency in units of Hz.
$\Omega$: Angular frqeuency, in units of radians/second; $\Omega = 2\pi F$
$f$: Normalized or "digital frequency". $f = F/F_s$, where $F_s$ is the sample rate in Hz. This is a formally unitless quantity, but sometimes its units are listed as being cycles/sample.  Sampling produces aliasing, which means that the spectrum of a sampled signal will be periodic with period 1 cycle/sample. For this reason DSP engineers only look at the spectrum of the signal in the principal period of $-1/2 \leq f < 1/2$.
$\omega$: Normalized radian frequency. $\omega = \Omega/F_s = 2\pi F/F_s$. Sometimes its units are listed as being radians/sample. Because of aliasing, it is only necessary to study the spectrum of a signal from $-\pi \leq \omega < \pi$ in the digital domain.

Other notions of the concept of "frequency" exist (such as instantaneous frequency, and a notion of an "instantaneous spectrum"), but these are more advanced concepts that you probably don't need to worry about yet.


How do DAWs align the phase when working with plugins?

As the author of a few AUs/VSTs an AU hosting app, I'll chime in here:

There is usually no phase alignment of signals passed through effects processing (AU/VST)
The "latency" of the effect is normally reported by the effect to the "hosting" environment on initialization. This way the hosting app can keep all channels in sync (i.e "wait" until all effects have returned the same time-range of audio before mixing, etc). Core Audio SDK states that AUs should implement this property kAudioUnitProperty_Latency if their latency is > 0 seconds.
Sometimes effects do not report their latency to the host app, so there are tools like Voxengo to work around that problem.

If the hosting environment doesn't have an automatic way to handshake & compensate for latency, then its compensation is a manual affair. For example Ableton Live has a feature (17.5 Device Delay Compensation) to help work around these timing issues.


Find the center of positive plateaus of rectangular function

Assume $x(t)$ is your original continuous function. The rectangular function $\pi(t)$ is generated by thresholding, something similar to

$$\pi(t)=1_{x(t)>\gamma}$$

where $\gamma$ is the threshold. Consider one of the rectangular segments of $\pi(t)$ which is between $t=t_a$ and $t=t_b$, where there is an upward edge at $t_a$ and a downward edge at $t_b$. The center point is at $$t_m=t_a+\frac{1}{2}(t_b-t_a)$$
So the problem reduces to finding the upward and downward edges of each rectangular segment. This can be done by looking at the first derivative (or if the functions are discrete, the first difference) of $\pi(t)$. The derivative has positive and negative impulses that correspond to upward and downward edges. In vector notation:

$$\mathbf{t}_a=\{t:\pi'(t)>0\}$$

$$\mathbf{t}_b=\{t:\pi'(t)<0\}$$
$$\mathbf{t}_m=\mathbf{t}_a+\frac{1}{2}(\mathbf{t}_b-\mathbf{t}_a)$$


yes ok thanks for the reply. Specifically, what I intend to do is an algorithm of digital image processing for use in retinal implants. Clarifying a bit my problem I know that there is a mathematical operation applicable to digital images that is called Wavelet Transform and that is in some way to compress the image can also lower its resolution through downsampling. My question is how can I from Matlab (with wavemenu GUI for example) compress an image using the Wavelet Haar transform and lower its resolution for example 640x480 to 32x32. In the following article is explained in a general way in the last step of table 1
https://www.researchgate.net/publication/323193463_Comparative_analysis_of_image_processing_algorithms_for_visual_prosthesis


Why is the backprojection of the Middlesbury depth map "cone shaped"?

The reason that pixels get wider with depth is that the depth map works in projection space, so each pixel is basically shaped like a frustum, not like a box.  Not sure where the noise is coming from but the depth map doesn't look perfectly clean either.


What is state of the art in multivariate 2D image display?

If "kind of works" but slightly less ugly is acceptable, than any two not-too-near-to-co-linear RGB color space vectors will do, and you or the user can just play with them to get a preferred combination of differentiation and non-ugliness (say, from a suitable artistic palate coordinated with other artwork color schemes in or around the data presentation, etc.) for the particular data sets.  Perhaps let the user choose one color vector, and auto rotate the other color vector to some roughly orthogonal point in the color space if it's too close to co-linear.

If you invert the data vector (or use a subtractive color vector) and thus fade to white, the lighter pastel version of some color vectors may blend better, but this depends on the mean and distribution of the data.


Looks like a convolution. S is the convolution of y with another signal of w+1 consecutive ones. So, S is the result of a moving sum filter. But the moving sum filter completely removes certain frequencies in the signal. For example, if you use w=3 and have the following signal y:

 y: ... 0  1  0 -1  0  1  0 -1  0  1  0  -1 ...

Then S[i] for any i (the moving sum over 4 neighboring samples) will always be zero, no matter how you shift or scale y. And that makes the problem of recovering y from s something between impossible and hard.

In case your data is finite, it is possible to write your equations as a matrix-vector product and solve it with the usual algorithms. You can ask tools like Matlab or GNU Octave about how noisy the inverse would be. This is something you should check for since the problem is ill-conditioned. Your equations, assuming y and s are finite can be written like this:

$$ A \cdot y = s $$

Here, A would contain mostly zeros and w+1 consecutive ones per row, at least for first rows. At the end, you have less ones because y is supposed to be finite. Then, asking for the condition number is a matter of typing

cond(A)

in Matlab or Octave. Example for 6 samples and w=2

octave:1> A = toeplitz([1 0 0 0 0 0],[1 1 1 0 0 0])
A =

   1   1   1   0   0   0
   0   1   1   1   0   0
   0   0   1   1   1   0
   0   0   0   1   1   1
   0   0   0   0   1   1
   0   0   0   0   0   1

octave:2> cond(A)
ans = 7.9053

Let's actually test it:

octave:3> y = [1 9 3 2 4 3]';
octave:4> s = A*y
s =

   13
   14
    9
    9
    7
    3

And now the inversion:

octave:5> y2 = A\s
y2 =

   1
   9
   3
   2
   4
   3

This inversion is actually easy to implement because the matrix exhibits a triangular structure. You can recover y sample by sample starting at the end and moving to the front.

This example was free of noise. Keep in mind that the condition number tells you how the relative error in s will be amplified (worst case) if you recover y from s. So, for example, if s contains a relative error of 3%, the recovered y may have a relative error as large as 3% * 7.9053 = 23.716%. Unfortunately, the condition number increases with the number of samples in this case. So, if you have to recover not six but 1000 samples and/or your w parameter is large, you're out of luck because the inversion will amplify the noise in s a lot, too.


finding the audio pulse compression peak in extra-low SNR scenario?

So this answer is coming from experience with radar systems which use complex signals (in-phase and quadrature signals so that phase information can be recovered), be aware that you may have to adapt it to fit your exact needs:

With radar systems, this is a pretty common problem. Sometimes, for a variety of reasons, the target that is present can’t rise above the noise, even with regular pulse compression gain. In this case, most systems turn to some sort of multi-pulse solution; a common one is called coherent integration. Basically, you send out and receive the exact same pulse N many times.

After doing this, you’ll have N many versions of the same target. Assuming no target migration has happened, you simply add all of the complex signals together. The idea here is that the target signal will cohere across all of them, such that a processing gain is achieved, but the noise will not cohere with itself and it’s power level will remain fairly constant.


Number of DFT (FFT) Points Required for a Specific Frequency Resolution for an Oversampled Signal

The resolution in the DFT is given by: $ \frac{{F}_{s}}{N} $.
Hence you need 10e6 / 100 = 100,000 samples to get the resolution you want.
You may bring the signal to the baseband (demodulation) and then you'll need a lower frequency of sampling to achieve what you want.
Since the effective bandwidth of the signal is (Two Sided) 50 [KHz], it can be sampled at $ {F}_{S} = 50 [KHz] $ and only 500 samples would be required for the asked resolution.


How to get an interpolation weight from a mathematical definition

I'll give a concrete example. Let's say you have a signal that was sampled every second. The frequency is $f = \frac{1}{1} = 1 Hz$.
Time: $T = \begin{bmatrix}0 & 1 & 2 & 3\end{bmatrix}$
Value: $X = \begin{bmatrix}1 & 2 & 3 & 4\end{bmatrix}$
We want to increase the frequency by $2$ i.e. $2 Hz$ (sample every 0.5 seconds).
Time: $\tilde{T} = \begin{bmatrix}0 & 0.5 & 1 & 1.5 & 2 & 2.5 & 3\end{bmatrix}$
Value: $\tilde{X} = \begin{bmatrix}1 & x_1 & 2 & x_2 & 3 & x_3 & 4\end{bmatrix}$
$x_1, x_2, x_3$ is determined by your interpolation functions $h_{nn}(\delta)$ and $h_{lin}(\delta)$. Here $\delta$ is time and both functions define intervals. Each value is given by $S(\delta) = \sum_{i=0}^{n-1} X_i \cdot h_{nn}(\delta-i)$. For nearest neighbor change the interval to $0 \leq \delta < 1$. Then $S(0) = S(0.5) = X_0 \cdot 1 + X_1 \cdot 0 + \cdots = X_0$. See also Cardinal B-splines.
While downsampling by $M$ requires strided convolution $$y[n] = \sum_k x[nM - k]h[k]$$
Upsampling needs fractionally strided convolution which is also called transposed convolution (see stackexchange)
$$y[j + nM] = \sum_k x[n-k]h[j+kM] \text{ and } j = 0, \dots, M-1$$
Transposed convolution with kernel size 3, stride 2 and padding 1 is equivalent to inserting 1 zero between inputs, pad by 1 and stride 1.
The kernel is $\begin{bmatrix}1 & 1 & 0\end{bmatrix}$ or $\begin{bmatrix}0 & 1 & 1\end{bmatrix}$ (either cross correlation or convolution) for nearest-neighbor interpolation (to double frequency):
from torch.nn import ConvTranspose1d
import torch
import numpy as np

def interpolate_nn(X):
  X = torch.from_numpy(X)
  with torch.no_grad():
    op = ConvTranspose1d(in_channels=1, out_channels=1,
                         kernel_size=3, stride=2,
                         bias=False, dilation=1, padding=1)
    op.weight.data = torch.tensor([0, 1, 1]).view(1, 1, -1).float()

    return op(X.view(1, 1, -1).float()).numpy().flatten()

X = np.array([1, 2, 3, 4])
print(interpolate_nn(X))

The result is [1. 1. 2. 2. 3. 3. 4.]
For linear interpolation use $\begin{bmatrix}0.5 & 1 & 0.5\end{bmatrix}$. The result is [1.  1.5 2.  2.5 3.  3.5 4. ]
Compare it with your $h_{lin}(\delta)$:
$\begin{align*}
S(0) &= X_0h_{lin}(0 - 0) + X_1h_{lin}(0 - 1) + \cdots = X_0(1 - |0|) = X_0\\
S(0.5) &= X_0h_{lin}(0.5 - 0) + X_1h_{lin}(0.5 - 1) + \cdots = 0.5X_0 + 0.5X_1\\
S(1) &= X_0h_{lin}(1 - 0) + X_1h_{lin}(1 - 1) + \cdots = 1X_1\\
\vdots
\end{align*}$


In a very simple form: Gaussian-like noises are frequent, and derivatives of squared functions are easy to minimize.

In very simple words: error, directly, cannot be used as a criterion. You can minimize a signed function to $-\infty$, which is not really useful here. As least, you need a cost function with a lower bound.

Power (if by that you mean the squared error) is an increasing function of the error, and positive. So if you minimize power, you can expect to minimize, somehow, a measure of the error (more precisely its absolute value).

The first problem is that the absolute value is not differentiable around zero, and can be difficult to optimize (one may look for a zero in the derivative). And when can you very easily find a zero in a derivative? When it is linear, because since Laplace, Gauss, Seidel such systems can be solved efficiently (matrix inversion). What are the functions that yield linear derivatives? Squares, parabolas...

Basically, for a while, squared functions were the only ones one could solve easily. The standard mean is just the minimum of the sum of differences squared. Interestingly, in linear cases, unrelated disturbances can add with a Gaussian distributions. And Gaussian distributions are well-related to squared errors.

But you are interested, recent advances in optimization allow to optimize more easily other types of functions, including absolute value loss criteria   or penalties.


So, you describe what you want to do relatively explicitly:

You care about the amplitude, not the sign of the signal.
you need to detect periods where that doesn't change much.
You want to quantize the result.
You want to Run-Length Encode (RLE) the result

So, to achieve each point:

abssignal = numpy.abs(signal)
Apply a low-pass filter. Simply design one using scipy.signal.
Quantization as ${0,\ldots,N}$: Divide by max(filtered_signal), multiplying by $N$ and then casting to integer (numpy.astype).
Apply RLE


Cross correlation vs FFT for finding phase between 2 signals

I don't see it as one or the other. You can use the FFT to more efficiently compute a cross-correlation (versus computing the cross-correlation in the time-domain).

There is lots of good information on this thread, about computing the "phase shift", from the offset of the peak in the cross-correlation function.


To answer your specifiy questions:

  Is the cv2.filter2D() function able to convolve two kernels?

Yes, but by default, it actually computes the correlation, not the convolution. However, the only difference is that the filter kernel (your second kernel) needs to be flipped. See OpenCV documentation for filter2D.

  Where does my logic/understanding begin to fail? (i.e. what topics/concepts am I completely not understanding?)

I think one more issue is the size of the resulting kernel.
The filter2D function returns an "output image of the same size and the same number of channels as src". However, the combined kernel increases is size (the combined size should be f1-size + f2-size - 1). So before you apply filter2D, you should add a black border to the first kernel.


50 KBaud means the symbol rate is 50 KHz (and in the OP's case the bit rate with 1 bit per symbol). The units of Hz is 1/sec. So the duration of each bit (the bit period) is simply

$$t_p = \frac{1}{50,000} = 20 \mu\text{s}$$


What real object could be modeled by this transfer function?

You transfer function models the behavior of a system. That may be a mechanical system, an electric system, etc. In your case you have a second order low pass filter.

If you calculate the step response of your system, you will notice that the system cannot follow the step. It will either grow fast, with a high slope and then start to decrease its slope when approaching the step, to never ever reach it (at least not in a finite lapse).

On the other hand, the system could grow faster to the step, reach it, then have an overshoot and start descending, reach again the step but continue descending, grow again, etc. This oscillation goes mathematically forever, or at least until its amplitude is negligible.

All this will depend on the position of the system poles.

Lets consider that the system models the behavior of, lets say, an elevator trying to reach the first floor. Here the mass of the elevator, the spring coefficients, etc. all can be used to model the system, and the step (input) is the wish to reach the first floor. The actual response of the system will be the step response, which will model the actual movement of the elevator reaching the first floor.


How to understand the sum of all "fourier frequencies"?

Let me put a practical answer with the following Matlab / Octave Code :

L = 2*1000;     % signal sample count
n = 0:L-1;      % discrete-time index

Fs = 44100;     % sampling frequency
am = [1, 0.5, 2, 0.5, 0.3, 0.6, 0.1, 0.2];                % magnitudes of 8 components
fm = [882, 2646, 4410, 6615, 8820, 10000, 13230, 15876];  % frequencies

% Time domain signal :
% x = a1*sin(2*pi*f1*n) + a2*sin(2*pi*f2*n) + ... + am*sin(2*pi*fm*n)
x = sin((2*pi/Fs)* n'*fm)*(am');

% Frequency domain signal :
X = abs( fft( x, L ) );

% Display :
figure,subplot( 2,1,1 );
plot( n/Fs , x );
title(' signal as a sum of sinusoidal frequencies ' );
xlabel(' time [s]');
subplot( 2,1,2 );
plot(linspace( 0 , Fs/2 , L/2 ) , (2/L)*X(1:L/2) );
title(' Fourier (frequency) spectrum magnitude of the        signal' );
xlabel(' Frequencies from 0 to Fs/2 [Hz]' );
ylabel(' FFT Magnitude ');

The output will be :

The output speaks for itself. The upper plot is the time-domain signal as a superposition of (M=8) different sinusoidal signals each with a particular frequency and amplitude. The lower plot shows its frequency spectrum magnitude which displays a peak per frequency component in $x[n]$. The peak's amplitude is proportional to the time-domain component amplitude.


Generally speaking, the BER only depends on the SNR. If the received signal power $P_\mathrm{S}$ is reduced due to path loss, while the noise power $P_\mathrm{N}$ remains constant, then the SNR will be reduced, because the SNR is defined as
$$
\gamma = \frac{P_\mathrm{S}}{P_\mathrm{N}}
$$
Consequently, the BER will be increased and the answer to your first question is: yes. However, this should not influence your BER vs SNR curves. Thus the answer to your second question is: no.

If your simulation behaves differently, i.e. if introducing a path loss "shifts" your BER curves, then the SNR is probably not calculated correctly. I suggest that you estimate the received signal power $\hat P_\mathrm{S}$ at the receiver (just before the sampler). Then the power of the additive noise should be calculated by
$$
P_\mathrm{N}=\frac{\hat P_\mathrm{S}}{\gamma}.
$$
Then add noise with mean power $P_\mathrm{N}$ to the received signal. With this setup, your BER/SNR curves should be independent of path loss.


I was going to post this as a comment, but then it became so long that I thought, this constitutes really as an answer in fact.

Your question reminds me of the lab-work in my MSc course on Adaptive Filters. We used "Wiener filter" to remove some unwanted background noise (wind and tire noise) from the recorded input to make the speech clearer. The Wiener filter is quite simple to begin with.

You can look at the wiener filter from the scipy package of python - scipy.signal.wiener.

Another advantage of Wiener filter could be that you don't have to do any additional computational step of taking the Fourier transform of the input. That is, Wiener filter can be applied to the time domain signal directly.

A short example of how to use Wiener would be:

from scipy.signal import wiener
import matplotlib.pylab as plt
plt.plot(t,signal,'k')
plt.plot(t,wiener(signal,mysize=55),'r',linewidth=3)
plt.show()

Going by the assumption that the "breathing sound" in your audio inputs would be quite uncorrelated to the actually speech signal, just like the wind and tire noise that we removed to make the speech clear in our audio data, I expect the filter to provided a reasonable output.

(Note: There were then additional processing steps added sometimes to improve the audio output like making the parameters of the filter adaptive. You have to play around with the parameters a bit, as in one case, I remember the Wiener filter produced an output like it was coming from a hollow pipe, if you chose/modified the parameters in  a particularly wrong way. But as a start, just applying the Wiener filter to your signal in Python would be a good first step. And you can improve from their onward, step by step.)

And if you are curious as to how this filter works then you catch up more on the math here.


Simulink Floating-Point and Fixed-Point Modeling

The model is fixed point if it contains a single block with the data type setting "fixdt()". If you want to override the fixed-point settings you should go in your Simulink model to Analysis then open "Fixed-point Tool". There you should change "Data Type Override" from "Use local settings" to something else.


Fourier Transform negative amplitude meaning

The fft returns complex values, to get the amplitude you need to take the abs( ). The real and imaginary portion tell you about the signals phase. Remember the fft is changing the basis by projecting your signal onto a complex sinusoid:
$$e^{i \omega t} = \cos(\omega t) + i \sin(\omega t)$$

and thus your signal is now a set of complex sinusoids which have some phase and amplitude. Think about the phase of a vector $v = [a, \, i\cdot b]$ on the complex plane and what this would mean.


Codec Output has "Wavy" Amplitude on Chirp Input

These could be ripples in the frequency response of an (ill-configured) anti-aliasing filter involved in the acquisition or restitution of the signal. It could be that the code running on your MCU and initializing the codec sends invalid configuration information to the codec; causing the codec to work with filters coefficients defined for a wrong sample rate. One thing worth noting is that the STM32F4 I2S master can only output an MCK at 256x FS. Check that your codec can work at 48kHz with such an MCK frequency.

Another thing I can think of is your software "pass-through" routine on the MCU dropping samples. Depending on the ratio between the block size and the period of the waveform, an amplitude modulation effect could appear.


How is GPS used to synchronize oscillators in a distributed system?

Part of the GPS signal is a 1 pulse per second clock that can be used in distributed systems for timing synchronization.

This patent has some diagrams that might help, e.g.

Regarding how accurate it is, I'm not sure. This paper has the following graph... which seems pretty good to me.


TL;DR: No, they are not necessarily the same.

Gory Details

Least squares is just an optimization technique. It is used in a variety of ways.

For filter design it is used to select that realizable filter $H_r(e^{j\omega})$ that most closely matches, in the least squares sense, the ideal required filter response $H_i(e^{j\omega})$:
$$
H_r(e^{j\omega}) = \arg \min \parallel H_r - H_i \parallel_2
$$
where $\parallel \cdot \parallel_2$ is the 2-norm or least-squares norm.

This sort of filter $H_r$ is not adaptive. That is, it doesn't change once it has been designed.

Adaptive filters may also use the least squares criterion, but in a different way: as part of the adaptation step.

Adaptive filters start off with initial filter coefficients $\vec{w}_o[0]$ and then use an update:
$$
\vec{w}_o[n] =\vec{w}_o[n-1] + \mu g[n-1]
$$
where $\mu$ is the step-size and $g$ is the gradient of the least squared error  surface in the direction of the minimum (from our current "location" of $\vec{w}_o[n-1]$).

Here, $g$ is determined by our error criterion: least squares. This means:
$$
\parallel \vec{w}_{\tt opt} - \vec{w}_o \parallel_2
$$
where $ \vec{w}_{\tt opt}$ is the unknown optimal (minimizing) solution.


Controllable realisation of $\frac{s^4+1}{4s^4+2s^3+2s+1}$ is both controllable and observable?

Controllability is a property of $A,B$ and observability is a property of $A,C$  matrices.

Transfer functions don't exhibit this property since if there is such a problem then you will not get a nonzero transfer function as there will be no possibility of control so nothing to transfer or nothing to observe. Hence, this is also a state space property.

The main theorem says that a state space realization is minimal if and only if it is both observable and controllable.

These representations allows one to quickly see, if exists, the unobservable and uncontrollable parts of the system. Such representations are obtained via the procedure called Kalman decomposition.

Here you have only constructed the controllable canonical form and already found out that the system is also observable. Thus the represantation is minimal.


How would Fourier and Cosine Transforms responds to summation of cosines with same frequency but different phases?

These integral transforms are linear, which means the transform of a sum is identical to the sum of the individual transforms, i.e.

$\mathcal F\left\{\cos(a)+\cos(b)\right\}= \mathcal F\left\{\cos(a)\right\}+\mathcal F\left\{\cos(b)\right\}$,

and also, that when multiplied with a constant, that constant can be "inside" or "outside" the transform:

$\mathcal F\left\{\alpha g(t) \right\}=\alpha\mathcal F\left\{g(t) \right\}\, \forall \alpha \text{ const.}$

No magic happening here; since (Euler's formula! It's important!):

$$\begin{align}
\mathcal F\left\{\cos(a)\right\}&=\mathcal F\left\{\frac 12\left(e^{i a}+e^{-ia}\right)\right\}\\
&=\frac12\mathcal F\left\{e^{i a}+e^{-ia}\right\}\quad\quad\quad\quad\quad||\,a=2\pi f t + \varphi\\
&=\frac12\mathcal F\left\{e^{i(2\pi ft + \varphi)}+e^{-i(2\pi ft + \varphi)}\right\}\\
&=\frac12\mathcal F\left\{e^{i 2\pi ft}\cdot e^{i\varphi}+e^{-i 2\pi ft }\cdot e^{-i\varphi}\right\}\\
&=\frac12\mathcal F\left\{e^{i 2\pi ft}\cdot e^{i\varphi}\right\}+\frac12\mathcal F\left\{e^{-i 2\pi ft }\cdot e^{-i\varphi}\right\}\\
&=\frac{e^{i\varphi}}2\mathcal F\left\{e^{i 2\pi ft}\right\}+\frac{e^{-i\varphi}}2\mathcal F\left\{e^{-i 2\pi ft }\right\}\tag{*}\\
\end{align}$$

we see that for any fixed phase $\varphi$ of a cosine, for example $\frac\pi 3$ or $\frac\pi 4$, the Fourier Transform "consists" of two dirac impulses $\delta$ with positions that only depend on the cosine's frequency $f$ and with amplitudes that depend only on the phase offset.

Now, just insert your different $\varphi$ in $(*)$ and add up those sums of Fourier Transforms to see the very basic point:

What you're modelling here, mathematically, is the interference of two waves at the same frequency. Tada! Physical significance added.

Now: good luck with the rest of your homework.


Algo Name: Detecting Elongated Shape in 2D Image by Convolution with Directed Lines

What you are essentially doing is a matched filter. However, thanks to Hough transform, your filter (line) is oriented and therefore I would call it an oriented matched filter. For generating the Bresenham line and sampling the pixels you might want the use the OpenCV line iterator. The simple usage would be similar to:

cv::LineIterator it(image, pt1, pt2, 8);
std::vector<cv::Vec3b> buf(it.count);
std::vector<cv::Point> points(it.count);

for(int i = 0; i < it.count; i++, ++it)
{
    buf[i] = *(const cv::Vec3b)*it;
    points[i] = it.pos();
}

The hole detection part can also be implemented using convolutions as follows:

Convolve the image with a horizontal convolution kernel composed of the oriented line. Use filter2D.
Threhsold the result to identify the areas with holes (use threshold).
Perform a blob-filtering and select the desired hole.

Even though I think that this approach is really valid and would work fine in many cases, the typical blob-person would approach the problem differently: First identify the large connected components in the image and filter out the rest. Within that single large component, perform a morphological hole filling. This is a pretty standard pipeline, and with that method, you will be able to find more tools (blob analysis) to help you. For example OpenCV methods such as findContours, drawContours or morphologyEx come handy.

Finally, if you know the shape of the hole, you might benefit from e.g. ellipse detection to directly find it.

Hope these help.


Scale and Rotation invariant feature descriptors

As far as alternatives to SIFT/SURF go, the question you linked provides very good answers.

There were two more questions I could read out:

"how could I build a useful (e.g. rotation invariant) feature descriptor"?
"regarding the statement from the linked question, how does he accomplish free rotational invariance?"

Building feature descriptors

This is a valid research topic. Good feature descriptors are not something just anyone can build in an afternoon. People publish articles when they successfully model feature descriptors with desirable properties. This is a reason that currently only a handful of state-of-the-art descriptors are used, and that is also what I advise you to do: find a feature descriptors that's good for your needs.

Achieving rotational invariance "for free"

You can determine the dominant gradient, or orientation at an image patch (your feature area). Then, rotate the image patch so that the gradient is always looking in the same direction, e.g $0$ (upwards). E.g. if you had a |black->gray->white| and a |white->gray->black| image, their dominant gradients would be pointing left ($-90$) and right ($90$) respectively, and when you rotate them by this amount, you get the same images.

This way you will always calculate the descriptor on an image patch with the same dominant orientation (the rotated patch), and thus you achieve rotational invariance.


Determine g(t) for minimal symbol error probability, given h(t), K=1 and Minimum distance receiver

Your computation using the energy of $h(t)$, $E_h$, to find the value of $c$, is correct. But note that you are only forcing the zero ISI condition for $k=0$. You still need to check the zero ISI criterion whenever $k\neq0$, but there is no mention of the value of $T$, and we can't derive what happens in those cases until we know the value of $T$.
For a general approach, I'd recommend going to the basics and not taking shortcuts. One key problem in your definition is that your convolution integral is not taking into account the time instant $t$, and therefore, you only compute one value of the convolution (for $t=0$), instead of all the values that are relevant (namely, $kT$ for all $k$).

The convolution of $h(t)$ and $g(t)$ is defined as
$$p(t) = (g* h)(t) = \int_{-\infty}^{\infty}g(\tau)h(t-\tau)\,d\tau = \int_{-\infty}^{\infty}h(\tau)g(t-\tau)\,d\tau$$
and we need to satisfy the Nyquist first criterion for zero ISI,
$$p(t)\Big|_{t=kT}=(g* h)(kT)=\left\{\begin{array}{rl}1,&k=0\\0,&k\neq0\end{array}\right.$$
Lastly, the matched filter is defined as $g(t)=ch^*(T-t)$ in its causal version, and $g(t)=ch^*(-t)$ in its non-causal version. For simplicity, let's use the noncausal version:
$$\begin{align}
p(t)& = (g* h)(t) = c\int_{-\infty}^{\infty}h(\tau)h^*(\tau-t)\,d\tau\\
\end{align}$$
In the particular case of $t=0$ (equivalently, $k=0$), we have
$$\begin{align}
p(0)& = (g* h)(0) = 1\\
&=c\int_{-\infty}^{\infty}h(\tau)h^*(\tau)\,d\tau= c\int_{-\infty}^{\infty}|h(\tau)|^2\,d\tau\\
&= c\Bigg(\int_{0}^{10^{-3}}|2|^2\,d\tau+\int_{10^{-3}}^{2\cdot10^{-3}}|1|^2\,d\tau\Bigg)=c\big(4\cdot10^{-3}+10^{-3}\big)=5c\cdot10^{-3}\\
&\Rightarrow c=10^3/5 = 200
\end{align}$$
Now, for verifying that the system satisfies the zero ISI criterion, we need the value of $T$. But you can follow a similar procedure.


Natural order would be to save it on frame by frame basis - that's the way you calculate them, otherwise you would need some gymnastics to do so. What's more, each 12 coefficients will be one vector in your 12D space - frame after frame.

Also it is easier to manage that in memory. In MATLAB like environments you simply use the matrices to store each vector in separate column. For C++, you can use array and array of pointers:

float **frames;
float  *data;
frames = new float*[num_frames];
data   = new float*[12 * num_frames];
memset(..)

After that you just set all frames pointers to first coefficient in each frame and you have easy way to access your features by using frames.

With interleaved features it wouldn't be so easy.


The DC coefficient, of the 2D-DCT (discrete cosine transform) of an 8 x 8 image block, represents the average value of the samples within the 8 x 8 block.

For your question, if you compute 8 different 1D-DCTs of 1 x 8 rows or 8 x 1 columns of this 8 x 8 block, then the DC coefficient of 2D-DCT of this 8 x 8 block will be the average of individual DC coefficients of 8 x 1 (or 1 x 8) 1D-DCTs.


Retaining filter state for processing data in chunks

You can use the filter function of matlab with initial conditions with the following syntax: (excerpted from Matlab 6.0 command window help)

[Y,Zf] = FILTER(B,A,X,Zi) gives access to initial and final conditions, Zi
and Zf, of the delays.  Zi is a vector of length MAX(LENGTH(A),LENGTH(B))-1
or an array of such vectors, one for each column of X.

Hope this helps.


If you add an accelerometer to the project, a Kalman filter can give a good estimation of vertical speed.

With only a barometric sensor, I don't think it's possible to reduce the lag below 1 second.

import numpy as np
import matplotlib.pyplot as plt
import random

from filterpy.kalman import KalmanFilter
from filterpy.common import Q_discrete_white_noise

random.seed(65537)

# Standard deviation of simulated sensor data
std_a = 0.075
std_h = 0.42

# Earth gravity
g_n = 9.80665

# Sensor sample period, s
dt = 0.05

kf1 = KalmanFilter(dim_x=2, dim_z=1)
kf1.H = np.array([[1, 0]])
kf1.F = np.array([[1, dt], [0, 1]])

kf2 = KalmanFilter(dim_x=3, dim_z=2)
kf2.H = np.array([[1, 0, 0], [0, 0, 1]])
kf2.F = np.array([[1, dt, dt * dt * 0.5], [0, 1, dt], [0, 0, 1]])

# initial process covariance
kf1.P = np.array([[1, 0], [0, 0.01]])
kf2.P = np.array([[std_h * std_h, 0, 0], [0, 0.0001, 0], [0, 0, std_a * std_a]])

# Process noise matrix
std = 0.0039
var = std * std
kf1.Q = Q_discrete_white_noise(dim=2, dt=dt, var=0.0625)
kf2.Q = Q_discrete_white_noise(dim=3, dt=dt, var=var)

# Measurement covariance
kf1.R *= np.array([[std_h * std_h]])
kf2.R *= np.array([[std_h * std_h, 0], [0, std_a * std_a]])

n = 300
r_n = 1.0 / n

t = np.zeros(n)
h_sim = np.zeros(n)
v_sim = np.zeros(n)
a_sim = np.zeros(n)
measured_h = np.zeros(n)
measured_a = np.zeros(n)
v_est_kf1 = np.zeros(n)
v_est_kf2 = np.zeros(n)

v0 = 0
v = v0
h = 0

for i in range(n):
    t[i] = (i * dt)
    a = 1.0 / 32 * g_n * np.sin(4 * np.pi * i * r_n)
    a_sim[i] = (a)
    v += a * dt
    v_sim[i] = v
    h += v * dt + (a * dt * dt) / 2
    h_sim[i] = h
    measured_a[i] = a + random.gauss(0, std_a)
    measured_h[i] = h + random.gauss(0, std_h)

# Compute the speed estimations

for i in range(n):
    kf1.predict()
    kf1.update(np.array([measured_h[i]]))
    v_est_kf1[i] = kf1.x[1]

    kf2.predict()
    kf2.update(np.array([[measured_h[i]], [measured_a[i]]]))
    v_est_kf2[i] = kf2.x[1]

# Plot the results

plt.figure(1, figsize=(8, 12), dpi=80)

plt.subplot(311)
plt.axis([0, n * dt, -0.75, 0.75])
plt.plot(t, measured_a, 'y+')
plt.plot(t, a_sim, 'r')
plt.title('Vertical acceleration - gravity, m/s^2')
plt.legend(('Measured', 'True (simulation)'), loc='best')

plt.subplot(312)
plt.axis([0, n * dt, -1, 7])
plt.plot(t, measured_h, 'c+')
plt.plot(t, h_sim, 'b')
plt.title('Altitude, m')
plt.legend(('Measured', 'True (simulation)'), loc='best')

plt.subplot(313)
plt.axis([0, n * dt, -0.5, 1.0])
plt.plot(t, v_sim, 'b')
plt.plot(t, v_est_kf1, 'r')
plt.plot(t, v_est_kf2, 'c')
plt.title('Vertical speed, m/s')
plt.legend(('True (simulation)', 'Kalman filter (1 sensor)', 'Kalman filter (2 sensors)'), loc='best')

plt.savefig('kalman_vsi.svg')
plt.show()

After some experiments with simulated sensor data, I wrote an Android app that uses smartphone sensors. Source code on GitHub


Precision delay is, in my opinion, best understood in terms of the Nyquist-Shannon sampling and reconstruction theorem.

If a continuous-time (a.k.a. "analog", but that is somewhat imprecise) signal $x(t)$ is bandlimited to below the Nyquist frequency $\frac{1}{2T}$ (where $T$ is the sampling period and $f_\mathrm{s} \triangleq \frac1T$ is the sample rate), then the original continuous-time signal $x(t)$ can be reconstructed from the discrete samples $x[n]$ as:

$$\begin{align}
x(t) &= \sum_{n=-\infty}^{\infty} x(nT) \, \operatorname{sinc}\left(\tfrac{t - nT}{T}\right) \\
&= \sum_{n=-\infty}^{\infty} x[n] \, \operatorname{sinc}\left(\tfrac{t - nT}{T}\right) \\
\end{align}$$

where

$$ \operatorname{sinc}(u) \triangleq \begin{cases}
\frac{\sin(\pi u)}{\pi u},  & \text{if } u \ne 0 \\
1, & \text{if } u = 0 \\
\end{cases} $$

All terms of the summation are bandlimited to a maximum frequency of $\frac{1}{2T}$, so the summation is bandlimited to the same bandlimit as the original $x(t)$.  The samples are

$$ x[n] \triangleq x(t) \Bigg|_{t = nT} \triangleq x(nT)$$

Because $x(t)$ can be evaluated at any $t$, not only multiples of $T$, this is the basis in which we interpolate to get a precision delay.

Suppose we wanna delay by an amount of time $\tau$, then the output of the delay is

$$\begin{align}
y(t) &= x(t-\tau) \\
&= \sum_{n=-\infty}^{\infty} x[n] \, \operatorname{sinc}\left(\tfrac{t-\tau - nT}{T}\right) \\
\end{align}$$

So that is explicitly where you get the delayed signal of $y(t) = x(t-\tau)$ from the samples $x[n]$.  The first problem we have is that we cannot add up an infinite number of terms.  So our first practical approximation is truncating terms.  But that is the same as applying a rectangular window (usually considered to be the worst one), so instead, we'll apply a better window.

$$\begin{align}
y(t) &= x(t-\tau) \\
&= \sum_{n=-\infty}^{\infty} x[n] \, \operatorname{sinc}\left(\tfrac{t-\tau - nT}{T}\right) \\
&= \sum_{n=-\infty}^{\infty} x[n-n_\tau] \, \operatorname{sinc}\left(\tfrac{t-\tau - (n-n_\tau)T}{T}\right) \\
&= \sum_{n=-\infty}^{\infty} x[n-n_\tau] \, \operatorname{sinc}\left(\tfrac{t-(\tau-n_\tau T)}{T}-n\right) \\
\end{align}$$

... for any integer $n_\tau$.  Let's choose it as the integer part of the dimensionless $\frac{\tau}{T}$, which is the delay in sample units.

$$ n_\tau \triangleq \left\lfloor \frac{\tau}{T} \right\rfloor $$

The fractional part of the delay (in sample units) is

$$ \frac{\tau}{T} - \left\lfloor \frac{\tau}{T} \right\rfloor $$

(i gotta return to this a little bit later.)