Audio texture analysis for visualization

Listening to that 'gritty' track you posted, two qualities stood out to me: firstly the bitcrushing / 8-bit effect, and secondly the overdriven quality of the guitar. Overall the effect sounded like there was all this 'noise' in the background, and the pitches of the instruments were competing with this.

Some thoughts based on that observation:

If you're already doing frequency detection, maybe you could you measure the proportion of energy at/near these frequencies, as a proportion of the total energy. This is almost like an estimate of SNR. So based on that, you could say that a higher 'SNR' is gritty sounding.
Alternatively, you could do some sort of frequency domain correlation with an ideal noise spectrum of some sort, e.g. white or pink noise.
Based on this page, both bitcrushing and overdrive create square waves in the waveform. So you could try and measure how square the waves are, perhaps by correlating the waveform with a superposition of square waves at your detected frequencies.


Your notes gave a counter example, you just proved it in the general case.
Your note used your x[n] as delta and showed that
y[n-1] unequal to H{x[n-1]} which is exactly what you showed.
In the notes example H{x}=0 for all n, as discrit delta equals zero foll all n but n=0, which means H is zero every where for all input of the kind delta[n-n0]  (becouae of the n factor which zeros y at its only non zero point)
Than you have zero output for delta input. But as you showed time varied output isn't zeros every where.


From your question it is unclear what is the dynamic range exactly, or is it needed to normalize the signal or not? Anyway, I assumed it to be normalized (since somehow this happens in ADCs). To quantize the normalized signal, you can use the following equation :
$$x_q[n]=\frac{\bigg\lfloor\left(\frac{x[n]}{\textrm{Dynamic Range}}\right)\cdot2^\textrm{Number of Bits}\bigg\rfloor}{2^{\textrm{Number of Bits}-1}}$$
where, $\textrm{Dynamic Range}$ here is $2$, $\textrm{Number of Bits}$ is $8$.

Quantization error in each sample would simply be the difference of the original signal and the quantized version of it.

$$Q_e=x[n]-x_q[n]$$

I can't provide the complete solution, however the following code gives you the idea of application of above discussions.

fs=8000;
Ts=1/fs;
f1=1500;
f2=2000;
f3=3000;
t=(0:512)*Ts;
x=2*cos(2*pi*f1*t+13124)+3*sin(2*pi*f2*t+pi/8)+7*cos(2*pi*f3*t);
x=x./max(abs(x)); %Normalization

Range=2; % Dynamic Range ; from -1 to 1 I assume
Q=8; %Number of Bits

Quantized_Signal=floor(x./Range*2^Q)./2^(Q-1);


Using Pilot tones to estimate carrier frequency offset in ofdm

I am familiar with 802.11a, so I'll answer your questions with regard to how 802.11a does it.

What are the nature of those Pilot tones?

There are 64 sub-carriers in 802.11a.  Four of them are used as the pilot channels.  The pilots themselves are simply known data- if I remember correctly they are either (1 + i) or -(1 + i).

how are they being generated?

The transmitter just puts the appropriate data into the pilot sub-channels and then ifft's the entire symbol (i.e. with all of the other sub-channels).

the authors of the paper wrote equation (18) without any derivation?

This isn't really a question.

can I get the system model for this technique?

What do you mean by a system model?  We don't normally provide code.


Understanding the units of wavelet time & frequency resolution

Confusion seems rooted in notation:

In $s = f_0 / (f \cdot \text{sampling rate})$, $f_0$ is the center frequency of mother $\psi$, and $f$ is the center frequency of the scaled $\psi$, i.e. $\psi_s$.
In $\sigma_{\hat\psi}=$ ..., $\omega_0$ is the (radian) center frequency of the wavelet of interest, or $\psi_s$, rather than $\omega_0=2\pi f_0$.
However, you can compute sigma of scaled from sigma of mother via: $\sigma_{\hat\psi_s} = \sigma_{\hat\psi} / s$ (see under 4.54).
$\text{sampling rate}$ is only relevant in discretized computations and must get canceled to yield physical units (Hz, sec, etc), so if you seek only physical use $s = f_0 / f_s$ (adjusted $f\rightarrow f_s$ to avoid confusion; recall that $f$ (or $\omega$) is what we integrate over, rather than plug in)
Do not involve $\text{sampling rate}$ in defining $\psi, \hat\psi$, scaled, or computing $\sigma$, or effectively most formulas you come across unless they explicitly include it, it's only to be used after computation, e.g. you first find $\sigma_{\hat\psi}$ without $\text{sampling rate}$, and then rescale as in $(4)$ here.
You commented "$\sigma_\omega$ Hz", but it's actually "rad-Hz"; to get Hz do
$f=\omega / (2\pi)$.

Let me know if this makes sense (if not, again with a specific example).


The description of the MRI signal from the magnetization vector

You always measure the (vector) sum over all magnetization vectors in your object, $S(t) = \sum\limits_{(x,y)} M(x,y,t)$. This sum gives you an amplitude and a phase. That is your measured signal.

From the gradients that you simulate, you know the $k$-space position as a function of time, i.e. you have $k(t)$. Now you sample your signal at certain points and sort the signals $S(t_n)$ into the respective $k(t_n)$-space. A FFT should now convert your $k$-space into the image again.


What is the difference between transport delay and dead time?

According to this document (p. 4) and this one (p. 13), transport delay and dead time denote the same thing: latency.


Digital has one major challenge for control application: latency.

Any digital system needs to include A/D conversion, processing and D/A conversion. The minimum amount of latency is 1/SampleRate even if you don't do anything else but a straight wire. You also typically can't use off-the-shelf sigma-delta-converters since these have significant amount of conversion latency.

Most digital control systems tend to use heavy oversampling and fast one-shot converters to work around this problem. Once you have the basic infrastructure in place, the digital domain offers a large variety of algorithms that are highly configurable and effect. However the cost of entry is high.


Deconvolution of an Image Acquired by a Square Uniform Detector

Your model is exactly a Convolution with Uniform Kernel where the output is what is called the Valid Part of the Convolution.

In MATLAB lingo it will be using conv2(mA, mK, 'valid').

So the way to solve it will be using a matrix form of the convolution and solving the linear system of equations.

Let's use the Lenna Image as input (Size was reduced for faster calculations):

We have a uniform kernel for the sensor model.
The output of the convolution with uniform kernel is given by:

The output from the sensor is both blurred and smaller (Less 2 rows and 2 columns as it is 3x3 kernel) just as in your model. This is the model of Valid Convolution.

In Matrix form what we have is:

$$ \boldsymbol{b} = K \boldsymbol{a} $$

Where $ \boldsymbol{b} $ is the column stack vector of the output image, $ \boldsymbol{a} $ is the column stack vector of the input image and $ K $ is the convolution operator (Valid Convolution) in matrix form. In the code it is done in the function CreateConvMtx2D().

So, now all we need is to restore the image by solving the Matrix Equation.
Yet the issue is the equation is Underdetermined System and the matrix has high condition number which suggest not to solve this equation directly.

The solution is to use some kind of regularization of the least squares form of the problem:

$$ \arg \min_{\boldsymbol{a}} \frac{1}{2} {\left\| K \boldsymbol{a} - \boldsymbol{b}  \right\|}_{2}^{2} + \lambda r \left( \boldsymbol{a} \right) $$

Where $ r \left( \boldsymbol{a} \right) $ is the regularization term.  In the optimal case the regularization should match the prior knowledge on the problem. For instance, in Image Processing we can assume a Piece Wise Smooth / Constant Model which matches the Total Variation regularization.

Since we have no knowledge here, we will use the classic regularization to handle the Condition Number - Tikhonov Regularization:

$$ \arg \min_{\boldsymbol{a}} \frac{1}{2} {\left\| K \boldsymbol{a} - \boldsymbol{b}  \right\|}_{2}^{2} + \frac{\lambda}{2}  {\left\| \boldsymbol{a} \right\|}_{2}^{2} =
{\left( {K}^{T} K + \lambda I \right)}^{-1} {K}^{T} \boldsymbol{b} $$

The output is given by (For $ \lambda = 0.005 $):

We can see that near the edge we have some artifacts which are due to the fact the system is Underdetermined and we have less equations to describe those pixels.
One can use the $ \lambda $ parameter to balance between how sharp the output is (Yet with artifacts) to how smooth it is, basically governing the level inversion of the system.

I advise playing with the parameter to find the best balance for your case but more than that, find a better regularization. Since the information you're after looks smooth you can use something in that direction.

The full MATLAB code is available on my StackExchange Signal Processing Q63449 GitHub Repository.

Enjoy...


Consider the function $\sin(2\pi f t)$. When $2\pi f t$ goes from $0$ to $2\pi$ you get oscillation of the sine wave, $2\pi$ to $4\pi$, another, and so on. So every time the argument changes by $2\pi$ you get one oscillation.

Now lets plot $2\pi f(t) t$ where $f(t) = 10 + 10 \sin(2 \pi t)$. (I've modified it a bit to show the effect more):

As $t$ increases, this value changes faster and faster, meaning $\sin(2\pi f(t) t)$ will have higher and higher frequency.

Another way of thinking about it - the frequency is the derivative of $f(t) t$. For the example this is $20 \pi (\sin(2 \pi t)+2 \pi \mathbf{t} \cos(2 \pi t)+1)$. Note the $t$ multiplier in there - $t$ increases thus the frequency increases.

So what you really want is
$\tfrac{d}{dt} f(t) t = 40 + 10 \sin(t)$ which gives $f(t) t = 40t - 10\cos(t) + c$ and for your example try $ \sin(2\pi(40t - 10 \cos(t)))$


Invertibility of Time-Dependent Fourier Transform

I have the 2nd edition, and there the case $m=0$ is considered, so the equation reads (Eq. $(10.21)$ on p. 716)

$$x[n] = \frac{1}{2\pi w[0]}\int_0^{2\pi} X[n, \lambda)d\lambda$$

(which of course assumes $w[0]\neq 0$). This makes sense because $e^{j\lambda m}$ equals $1$ for $m=0$.


Gradient of Total Variation of Magnitude of Complex Function for Denoising

The problem with $\left|f\right|$ is that since is not analytic the standard definition of complex derivative does not apply. A solution is to use Wirtinger derivatives:

http://en.wikipedia.org/wiki/Wirtinger_derivatives

A detailed account of Wirtinger calculus for signal processing problems is

http://arxiv.org/abs/0906.4835

Another (probably simpler) option is to treat the complex image as a two-channel (real,imag) image and use the definition of derivative for vector fields. This paper has a very clear explanation on how to do this:

Lee, H.-C.; Cok, D.R.; "Detecting boundaries in a vector field" (IEEE Transactions on Signal Processing, vol.39, no.5, pp.1181-1194, May 1991)


Statistic that is responsive to changes in time series, yet not too volatile

Hi: I don't think that there are general statements that can be made because the median component used in the second approach is pretty ad-hoc.

You probably already know this but the $\alpha$ that you defined for the EWMA is such that your EWMA has the same weighted age as just taking a moving average of length "window". Therefore, if you an moving average of length window and an EWMA with parameter $alpha = \frac{2}{window + 1)$, they will look quite similar.

The way you look at things in terms of "reaction time" till convergence or catching up is interesting.

The standard approach for this ( and more systematized but I'm not critcizing your creative approach ) is to just use an ewma with parameter $\alpha$.
The different values of $\alpha$ chosen will equate to the half life of an impulse response for $S_t$. So, just by plotting the impulse response, you can see what $\alpha$ you prefer.

So, for example, suppose you used the EWMA with parameter $\alpha$. ( $\alpha$ is between zero and one for invertibility purposes ) Therefore,
your exponentially smoothed value, $\hat{S}_{t}$, at every time step, $t$ would be:

$\hat{S}_t = \alpha \times \hat{S}_{t-1} + (1-\alpha) \times X_{t} $.

Now,  we assume that $S_{0} =1 $ and $X_{t} = 0 $ for $t$ greater than zero and that we start start smoothing $S$ at $t=1$. So, since  the halflife corresponds to the period when $\hat{S}_t$ becomes half of what it was originally, this corresponds to

$\hat{S}_{halflife} = \frac{1}{2}$.

Therefore, we have $\frac{1}{2} = \alpha^{halflife} $.

Taking the log of both sides gives $halflife = \frac{log(1/2)}{log(\alpha)}$.

So, the point is that, given the impulse response at time $t=0$, then different $\alpha$ correspond to different halflives
for the response and you can play around with the halflife to see what works best for your problem.

Finally, if this was obvious or you were already familiar with it, my apologies. Also, IMHO, the best textbook exposition of this is Brown's "Smoothing, Forecasting and Prediction of Discrete Time Series". There are used versions on Amazon that are pretty cheap.


How to find the output of a system with the unit-step response and the input?

This can indeed be done. Write the input function $x(t)$ as

$$x(t)=x(-\infty)+\int_{-\infty}^t x^{\prime}(\tau)d\tau=
x(-\infty)+\int_{-\infty}^{\infty} x^{\prime}(\tau)u(t-\tau)d\tau\tag{1}=\\
=x(-\infty)+(x^{\prime}*u)(t)$$

where $x^{\prime}(t)$ is the derivative of $x(t)$, and $u(t)$ is the unit step function. From (1), the system's response $y(t)$ to the input $x(t)$ can be written in terms of the step response $a(t)$:

$$y(t)=a(\infty)x(-\infty)+\int_{-\infty}^{\infty} x^{\prime}(\tau)a(t-\tau)d\tau\tag{2}$$

So, apart from the first constant term, the output signal $y(t)$ is given by the convolution of the derivative of the input signal with the system's step response.

Equation (2) is valid because the system's response to the constant $x(-\infty)$ is $H(0)x(-\infty)$, where $H(0)$ is the system's frequency response at DC:

$$H(0)=\int_{-\infty}^{\infty}h(t)dt=a(\infty)$$

The integral on the right-hand side of (2) is the convolution

$$((x^{\prime}*u)*h)(t)=(x^{\prime}*(u*h))(t)=(x^{\prime}*a)(t)$$

EDIT: The discrete-time case is completely analogous. Write the input sequence as

$$x[n]=x[-\infty]+\sum_{k=-\infty}^n(x[k]-x[k-1])=x[-\infty]+\sum_{k=-\infty}^{\infty}(x[k]-x[k-1])u[n-k]$$

This results in the following expression for the output sequence in terms of the step response:

$$y[n]=x[-\infty]a[\infty]+\sum_{k=-\infty}^{\infty}(x[k]-x[k-1])a[n-k]$$

Again, apart from the constant first term, the output sequence is obtained by convolving the first order difference of the input sequence with the system's step response.


Spectral plot shows more notes than there really should be

Since the notes from a guitar are no pure sinusoids, you should expect to see some harmonics, even when analyzing the dry signal without effects. E.g., the note E is the perfect fifth of the note A, i.e., it is the second harmonic.

If you use distortion or modulation effects (chorus, flanger, and phaser) you get even more additional frequencies due to the non-linearity (distortion) or due to time-varying filtering (modulation effects).

So in sum there are three sources of additional frequency components different from the fundamental pitches:

harmonics (due to non-sinusoidal signal)
non-linearities (e.g., overdrive/distortion)
time-varying filtering (e.g., modulation effects)

And in your example you get all three of them.


Impulse response of the inverse system to the backward difference

You're on the right track. Let's call the impulse response of the second system (equalizing the first one) $h_2(n)$. The trick is to assume causality of the second system, i.e. $h_2(n)=0$ for $n<0$. By doing so, the solution for $h_2(n)$ becomes unique.

The total impulse response of the cascaded system is

$$h(n)=h_2(n)-h_2(n-1)\tag{1}$$

and this must equal $\delta(n)$. So we get

$$h_2(0)=1\\
h_2(n)=h_2(n-1),\quad n\ge 1$$

Which gives

$$h_2(1)=h_2(0)=1\\
h_2(2)=h_2(1)=1\\
h_2(3)=h_2(2)=1\\\vdots$$

So $h_2(n)=u(n)$, where $u(n)$ is the discrete-time step function.

Note that you can also find an anti-causal filter $h_2(n)$ satisfying the requirement $h(n)=\delta(n)$. This anti-causal solution is $h_2(n)=-u(-n-1)$. You can easily convince yourself that for this choice of $h_2(n)$ the total impulse response $h(n)$ given by (1) also equals $\delta(n)$.

Note also that all the solutions you gave at the end of the question have the form $h_2(n)=u(n)+c$ where $c$ is some constant. It is indeed true that all these systems are solutions to the problem because the frequency response of the first system has a zero at DC, compensating for any DC component of $h_2(n)$. However, all of these systems with $c\neq 0$ are neither causal nor stable.


How is a "good" or "bad" singer defined in signal processing terms?

First order, some DSP analysis might be able to detect central pitch, pitch variation (vibrato), timbre (overtone richness), and note onset, and say whether these are grossly out-of-tune, thin, or off-time, by how much and how often.

However, second-order, good musicians play with subtile variations in pitch (equal vs. just intonation, barbershop quartet harmonization as the key changes, etc.) and in timing (jazz, blues, "swing" rhythms, etc.).  How to "score" the latter is probably more a topic for computer AI than DSP.


Can realizations of a filtered Gaussian white noise process be represented as a Fourier transform?

I don't think it makes sense to talk about the Fourier transform of a (realization of a) random process, because it generally doesn't exist if no further restrictions are imposed. This problem is also discussed here.

What can be done is analyze the statistics of the discrete Fourier transform (DFT) and the discrete-time Fourier transform (DTFT) of a windowed version of a discrete-time stationary random process. This has been done in this paper.


Finding symmetric regions/patterns in image

As I've said in the comments, medical image registration is a topic with lots of research available, and I'm not an expert. From what I've read, the basic idea commonly used is to define a mapping between two images (in your case an image and its mirror image), then define energy terms for smoothness and for image similarity if the mapping is applied, and finally optimize this mapping using standard (or sometimes application-specific) optimization techniques.

I've hacked together a quick algorithm in Mathematica to demonstrate this. This is not an algorithm you should use in a medical application, only a demonstration of the basic ideas.

First, I load your image, mirror it and split these images into little blocks:

src = ColorConvert[Import["http://i.stack.imgur.com/jf709.jpg"],
   "Grayscale"];
mirror = ImageReflect[src, Left -> Right];
blockSize = 30;
partsS = ImagePartition[src, {blockSize, blockSize}];
partsM = ImagePartition[mirror, {blockSize, blockSize}];
GraphicsGrid[partsS]

Normally, we would do an approximate rigid registration (using e.g. keypoints or image moments), but your image is almost centered, so I'll skip this.

If we look at one block and it's mirror-image counterpart:

{partsS[[6, 10]], partsM[[6, 10]]}

We can see that they're similar, but shifted. The amount and direction of shift is what we're trying to find out.

To quantify the match similarity, I can use the squared euclidean distance:

ListPlot3D[
  ImageData[
   ImageCorrelate[partsM[[6, 10]], partsS[[6, 10]],
    SquaredEuclideanDistance]]]

sadly, using this data is the optimization directly was harder than I thought, so I used a 2nd order approximation instead:

fitTerms = {1, x, x^2, y, y^2, x*y};

fit = Fit[
   Flatten[MapIndexed[{#2[[1]] - blockSize/2, #2[[2]] -
        blockSize/2, #1} &,
     ImageData[
      ImageCorrelate[partsM[[6, 10]], partsS[[6, 10]],
       SquaredEuclideanDistance]], {2}], 1], fitTerms, {x, y}];

Plot3D[fit, {x, -25, 25}, {y, -25, 25}]

The function is not the same as the actual correlation function, but it's close enough for a first step. Let's calculate this for every pair of blocks:

distancesFit = MapThread[
   Function[{part, template},
    Fit[Flatten[
      MapIndexed[{#2[[2]] - blockSize/2, #2[[1]] - blockSize/2, #1} &,
        ImageData[
        ImageCorrelate[part, template,
         SquaredEuclideanDistance]], {2}], 1],
     fitTerms, {x, y}]], {partsM, partsS}, 2];

This gives us our first energy term for the optimization:

variablesX = Array[dx, Dimensions[partsS]];
variablesY = Array[dy, Dimensions[partsS]];

matchEnergyFit =
  Total[MapThread[#1 /. {x -> #2, y -> #3} &, {distancesFit,
     variablesX, variablesY}, 2], 3];

variablesX/Y contains the offsets for each block, and matchEnergyFit approximates the squared euclidean difference between the original image and the mirrored image with the offsets applied.

Optimizing this energy alone would give poor results (if it converged at all). We also want the offsets to be smooth, where the block similarity tells as nothing about the offset (e.g. along a straight line or in the white background).

So we set up a second energy term for smoothness:

smoothnessEnergy = Total[Flatten[
    {
     Table[
      variablesX[[i, j - 1]] - 2 variablesX[[i, j]] +
       variablesX[[i, j + 1]], {i, 1, Length[partsS]}, {j, 2,
       Length[partsS[[1]]] - 1}],
     Table[
      variablesX[[i - 1, j]] - 2 variablesX[[i, j]] +
       variablesX[[i + 1, j]], {i, 2, Length[partsS] - 1}, {j, 1,
       Length[partsS[[1]]]}],
     Table[
      variablesY[[i, j - 1]] - 2 variablesY[[i, j]] +
       variablesY[[i, j + 1]], {i, 1, Length[partsS]}, {j, 2,
       Length[partsS[[1]]] - 1}],
     Table[
      variablesY[[i - 1, j]] - 2 variablesY[[i, j]] +
       variablesY[[i + 1, j]], {i, 2, Length[partsS] - 1}, {j, 1,
       Length[partsS[[1]]]}]
     }^2]];

Fortunately, constrained optimization is built-in in Mathematica:

allVariables = Flatten[{variablesX, variablesY}];
constraints = -blockSize/3. < # < blockSize/3. & /@ allVariables;
initialValues = {#, 0} & /@ allVariables;
solution =
  FindMinimum[{matchEnergyFit + 0.1 smoothnessEnergy, constraints},
   initialValues];

Let's look at the result:

grid = Table[{(j - 0.5)*blockSize - dx[i, j], (i - 0.5)*blockSize -
      dy[i, j]}, {i, Length[partsS]}, {j, Length[partsS[[1]]]}] /.
   solution[[2]];
Show[src, Graphics[
  {Red,
   Line /@ grid,
   Line /@ Transpose[grid]
   }]]

The 0.1 factor before smoothnessEnergy is the relative weight the smoothness energy gets in relation to the image match energy term. These are results for different weights:

Possible improvements:

Like I said, perform a rigid registration first. With a white background, simple image moments-based registration should work fine.
This is only one step. You can use the offsets you found in one step and improve them in a second step, maybe with a smaller search window or smaller block sizes
I've read articles where they do this without blocks at all, but optimize an offset per pixel.
Try different smoothness functions


Are there any signal processing libraries that implement TD-PSOLA?

When I'm looking for source code, I generally search the open repositories rather than Google.

Try this one in C++.

However, there don't seem to be many on github, and sourceforge has nothing of interest.


Compressed Sensing vs Common Compression Approaches

Compression works by reducing redundancy and there are a number of ways that this can be done.

Very briefly, consider the pair of functions $C = f(S_{original}, \Theta_{f}), S_{recovered} = g(C, \Theta_g)$, where $f,g$ are the compressor and decompressor functions respectively, $\Theta$ is the set of parameters specific to $f,g$ and $C, S_{original}, S_{recovered}$ are the compressed, original and decompressed sequences respectively. If the compression is loseless then the error between the original and recovered sequences should be zero ($S_{recovered}-S_{original}=0$).

However, we might accept some small error so that $S_{recovered}-S_{original}<=e$. In this case, the compression is lossy.

What $f$ does is to reduce the redundancy in $S_{original}$ so that it can be represented using less data so that the length of $C$ is less than the length of $S_{original}$ ($|C| < |S_{original}|$).

Run-length encoding achieves this by looking for sequences of repeated values and substituing them for a pair of "value, count". The redundancy here is most obvious.

Huffman encoding, does this by creating a carefully constructed "dictionary" of codewords based on knowledge about the source of a signal. In this, frequently occuring values are mapped to short codewords and infrequently occuring values are mapped to longer codewords. In this way, instead of representing all values with the same amount of data (for example 8 bits), frequent values are represented with less data (less than 8 bits) and on average lead to describing a long $S_{original}$ with a short $C$. More frequently appearing characters than others, again, imply redundancy.

Lempel-Ziv-Welch encoding reduces redundancy by looking for repeating sequences. Once a repeating sequence is found it is placed in memory and assigned a codeword. The next time the same sequence is found it is substituted with its codeword. In this way, long sequences (of more than one values) are substituted with just one value (an index into the dictionary).

Other methods rely on preserving just a part of the frequency spectrum of a signal or image therefore reducing the overall ammount of data that needs to be described.

The key thing to remember here is redundancy. Repeating values, slow varying signals (which again implies nearby samples having similar values) or even components with negligible contribution. If a signal contains redundancy then it is predictable and if it is predictable then it can be captured in a model. (As a side note, please note that this redundancy, or amount of information in a signal can be quantified by Entropy).

So, compressed sensing is in a way similar to lossy compression in the sense that we are trying to recover an $S_{original}$ from a (shorter) $C$ accepting some small error $e$ in the recovery with some "decompressor" $g$.

BUT, to do this we need to make certain assumptions about the dynamics of $S_{original}$. For example, it is varying slowly, therefore it is not necessary to sample the full spectrum of frequencies that could generate some potential $S_{original}$

Another way to think about this is to look at a picture or a drawing and cut (or black)out a rectangle from it. That is, throw away some data or, don't acquire a part of data from the beginning. If the picture is predictable enough then it is easy to reconstruct the missing part. WHY is it predictable? Because it is a repeating pattern and we assume that the missing bit behaves similarly. If the picture is less predictable then we can still attempt to reconstruct a rectangle of it from the surrounding information but perhaps it will not be as accurate.

If we were to take this to the extreme (as it is done in practice anyway), then we could say that we don't really need to acquire all of the pixels of the image. We can acquire some of the data and by exploiting the redundancy of the signal and the fact that it is not going to be too unpredictable (i.e. we can model it) reconstruct the missing bits...almost as closely to reality.

This is exactly what iterative reconstruction techniques (a family of $g$ "decompressors") are doing by progressively trying to minimise the error between an $S_{recovered}$ at some iteration $n$ and the subset of known values contained in the actually acquired $C$ and stopping this process once the error drops below an "acceptable" $e$.

For more information please see this link which is a relatively simple iterative reconstruction scheme that shows how to "update" the unknown values from the known values. For even more information please this link.

Hope this helps.


Different methods exist for estimating PSDs. They can be broadly classified as parametric or non-parametric methods. Indeed you should learn which method they used to obtain the associated PSD.

An FFT based PSD estimation is a non-parametric method such as Periodogram and its variants. This method will produce random looking non-smooth output result which can be smoothed, but this has limits and has side effects too.

On the other hand, parametric estimators can produce very smooth and high resolution looking PSD estimates. There are many parametric PSD estimators. The authors probably used a model based parametric estimator to obtain the plotted PSD, that's why they are smooth looking and that's what I understand from your word curved ?


Filter odd or even harmonics with notch or inverse notch filter

What you are looking for are what we, in the audio space, call comb filters.  Comb filters may or may not have a feedback path, just like FIR and IIR filters. In fact, there is a generalized theory about designing comb filters based simply on designing any digital filter.
In fact, you can think of every LTI digital filter as being a comb filter where the teeth of the comb are spaced apart by the sampling frequency.  To space the teeth apart by a frequency spacing much less than the sampling frequency, you need a delay element that is longer than a single-sample delay ($z^{-1}$).  If a delay with precision of an integer number of samples is good enough to tune your comb filter, good for you.  But if it isn't, you need to deal with interpolation and precision delay.


Does the frequency response assume the zero input response to be 0?

A frequency response describes a linear time-invariant (LTI) system. It is the Fourier transform of the system's impulse response, and only LTI systems are fully characterized by their impulse response.1

A system with non-zero initial conditions is not LTI because its output has a component which only depends on the initial conditions and which is independent of the input signal. Consequently, a system with non-zero initial conditions cannot fully be described by a frequency response. If you assume zero initial conditions, then such a system (if otherwise LTI) can be described by a frequency response.

So to answer your question, a system can indeed only be fully described by a frequency response if its zero-input response (ZIR) is zero, because if its ZIR were not zero it wouldn't be an LTI system.

Also take a look at these two answers to related questions: answer 1, answer 2.

1Also linear time-varying (LTV) systems can be characterized by an impulse response, but it is a two-dimensional function, and LTV systems don't have a frequency response in the conventional sense.


Method to rescale signals to mean length

If you really want to scale data, you can look at the scale transform, detailed in The scale representation, L. Cohen, 1993. There also exists a time-scaling method given in A. Papoulis, "Signal Analysis", McGraw-Hill, p. 290, 1977; one Matlab/Scilab/Octave implementation is provided in the DiscreteTFDs toolbox by  Jeffrey C. O'Neill.For the latter, you just have to provide a scaling factor proportional to the inverse of the signal's length.

However, in the same line as @Florian, and looking at the signals, I would be more prone to 1) find the maxima at fractional steps through interpolation (eg Quadratic Interpolation of Spectral Peaks) 2) realign with fractional filters (Splitting the unit delay: tools for fractional filter design) 3) truncate signals.


Local Histogram Processing clarification

If you move the 3x3 window to the right by 1 pixel, a new column gets introduced to the right of the previous window and the leftmost column gets deleted.

If the histogram component for the k$^{th}$ intensity level is $p_{k}=n_k/n$ where n is the total number of pixels in the neighborhood and $n_k$ is the number of pixels with intensity k.

For new window, $p'_{k} = (n_k+n_{r_k}-n_{l_k})/n$
where $n_{r_k}$ is number of occurrence of the $k^{th}$ intensity level pixels in the added right column and $n_{l_k}$ is number of occurrence of the $k^{th}$ intensity level pixels in the deleted left column.

$$p'_{k} = n_k/n+(n_{r_k}-n_{l_k})/n$$
       $$p'_{k} = p_k+(n_{r_k}-n_{l_k})/n$$


Deconvolution of shifted gaussian function in the frequency range

It is naive, but you didn't tell why you don't use $s(t) = A(t) \cdot \exp(-i \omega_0 t)$, and calculate $A(t) = s(t) \cdot \exp(i \omega_9 t)$...
I named your signal $s(t)$ here.
Other not so naive things.
You mention number of points, so I assume you are discretizing the signal. And your Fourier transform is a finite length DFT. In that case the Fourier transform will have alias, if you sample rate is enough to capture all the information in $A(t)$ then you can simply calculate $S(\omega) = DFT(s(k T_s))$, and the $A(\omega) = S(\omega + \omega_0~ \textrm{mod}~ F_s)$. Where $T_s$ is the sampling period, and $F_s$ is the sampling frequency, thus $s(k T_s)$ is the k-th sample of your signal.


how do you even measure the frequency of a digital signal?

This may be the source of your confusion. In digital communications, data is transmitted using an analog, continuous-time waveform. This is necessary because only this kind of waveform can exist in the physical world. Like all waveforms, it has a bandwidth, and its Fourier Transform determines what frequencies are present in it.

In digital communications, the waveform used to communicate has a specific form: it is a sequence of pulses, and the amplitude of the pulses encodes the transmitted bits.

The simplest waveform is one that encodes a bit $1$ with a pulse of positive amplitude, and a bit $0$ with a pulse of equal but negative amplitude. You may think that, for random bits, the negative pulses cancel out with the positive pulses and the average DC level of the signal is 0. This is in general not true over short time periods; even random bits will have sequences with many more $1$'s than $0$'s, or viceversa. This differential causes the signal's DC to grow, which causes distortion at the receiver end.

Line encoding can help with this problem. For example, the Alternate Mark Inversion line code encodes a bit $0$ to a pulse of amplitude zero, and a bit $1$ to pulses of positive and negative amplitude, alternating. So, every positive pulse is followed by a negative pulse, and the signal's DC never builds up.


Impulse response (general form for linear systems as two-variables function $h(t,\tau)$) applied to Time-invariant systems

If the response to $x(t)$ is given by

$$y(t)=\int_{-\infty}^{\infty}h(t,\tau)x(\tau)d\tau\tag{1}$$

then the response to $x(t-T)$ is

$$\tilde{y}_T(t)=\int_{-\infty}^{\infty}h(t,\tau)x(\tau-T)d\tau=\int_{-\infty}^{\infty}h(t,\tau+T)x(\tau)d\tau\tag{2}$$

If the system is time-invariant we require

$$\tilde{y}_T(t)\stackrel{!}{=}y(t-T)=\int_{-\infty}^{\infty}h(t-T,\tau)x(\tau)d\tau\tag{3}$$

For $(2)$ and $(3)$ to be equal for any $x(t)$ we must have

$$h(t,\tau+T)\stackrel{!}{=}h(t-T,\tau)\tag{4}$$

which is equivalent to

$$h(t+T,\tau+T)\stackrel{!}{=}h(t,\tau)\tag{5}$$

From $(5)$ it is clear that for a linear and time-invariant system the value of $h(t,\tau)$ only depends on the difference of its arguments, and not on their individual values. Consequently, the impulse response can be rewritten as a function of $t-\tau$:

$$h(t,\tau)=\tilde{h}(t-\tau)\tag{5}$$


Let us say you detected noise level power at -100dBm. You keep additional 10dB threshold to detect energy. Means if you detect >= -90dBm power, you try to detect preamble of Wifi packet. If you detect -95dBm power, you ignore it assuming it may not be a valid Wifi packet. This really works assuming you are very far from base station/access point(AP). Now suppose say you want a higher threshold so that you don't want too many false alarms (may be in order to save battery). So you raise the threshold to -60dBm. Here, if you are too far from AP, you may not trigger at all. If you are close to say within 10 meters, you start to detect energy more than -60dBm and try to detect Wifi packets.

So, you can decide on threshold on how frequently you want to try to detect Wifi signal.

If your threshold is too low, you may more frequently try to detect packets. Advantage - you can detect packets if you are far away from AP. Disadvantage- you may consume more power trying to detect packets.

If you increase threshold, you will detect packets less frequently. Advantage - you can filter out signals which are weak and can use higher modulation rates. Disadvantage - you cannot move too far away from the AP.


If the impulse response is $h[n], n\geq 0$, then the step response $s[m]$ at any time $m≥0$ is
$$s[m]=h[0]+h[1]+h[2]+\cdots +h[m].$$
Note that if the filter is a finite impulse response (FIR) filter so that
$h[m] = 0$ for $m > M$, then the step response at $M$ is
$$s[M] = h[0]+h[1]+h[2]+\cdots +h[M]$$
and remains at this value forever afterwards, that is, $s[m]=s[M]$ for all $m \geq M$.

What other information do you want to get? The transfer function of the filter? That's easy:
$$H(z)=h[0]+h[1]z^{−1}+h[2]z^{−2}\cdots $$
which works out to be polynomial of degree $M$ in $z^{-1}$ for a FIR filter.


Connection with system analysis and laplace&Z transform

You are right that the (bilateral) Laplace transform can be interpreted as the Fourier transform of $e^{-\sigma t}f(t)$. However, I think that the significance of the Laplace transform only becomes clear when $s=\sigma+j\omega$ is viewed as a complex variable because then we can study the analytic properties of the system function. E.g., electrical networks are characterized by the poles and zeros (in the complex $s$-plane) of their system function.

The practical importance of the Laplace transform for the analysis of continuous systems lies in the fact that it transforms linear differential equations with constant coefficients into algebraic equations, which are obviously much easier to solve. This is achieved by the fact that in the Laplace transform domain, integration is equivalent to multiplication by $1/s$, and differentiation becomes multiplication with $s$.

Any linear time-invariant (LTI) system, such as the RLC-circuit you mentioned, will result in linear differential equations with constant coefficients. This is why any LTI system can be analyzed by the Laplace transform. Note that in practice usually the unilateral transform is used (where the integral starts at $t=0$ instead of $t=-\infty$), because practical systems are causal, and because in this way initial conditions (such as charged capacitors at $t=0$) can be easily taken into account. The unilateral Laplace transform is thus an appropriate tool to solve initial value problems.

Everything said so far carries over to the $\mathcal{Z}$-transform if you replace 'continuous system' by 'discrete system', and 'differential equation' by 'difference equation'. Linear time-invariant discrete systems consisting of adders, multipliers and delay elements are treated by the $\mathcal{Z}$-transform in the same way as LTI electrical networks are treated by the Laplace transform. And indeed, IIR filters are analogous to RLC circuits in the sense that both have an infinitely long impulse response and that both can be characterized by the poles and zeros of their system functions.


Two-sided impulse response and its inverse in real life

The analog front end will always contain some sort of low pass filter. This can be chosen to be either linear phase or minimum phase. Linear phase preserves the phase response of the original signal but the time domain response is non-causal, i.e. the impulse response is symmetric around the center tap.

Minimum phase has a one sided impulse response, with the strongest tap up front, but it incurs some phase distortion. Ideally we want both, but we have to pick one.

The same consideration implies to the inverse: the inverse of a linear phase filter is linear phase and the inverse of a minimum filter is minimum phase. However there is an added complexity: not all filters are invertible: If the magnitude of the original transfer function is very small or zero at some frequencies, the inverse is not well defined. You get "divide by zero" or "divide by noise" problems.


What Is the Cause for Poor Results with Adaptive Thresholding?

Some guidelines for using Thresholding:

Stretch the image to use the whole Dynamic Range (DR).
Apply some Denoising (Very very gentle).
Median with small radius would be a good idea.
Unless you hand tweak the Threshold, Otsu's Method generally yields good results for this kind of tasks (Text on background).
If one use Adaptive Local Methods (Mean / Gaussian based) the Radius is another hyper parameter to tweak.


Can cross-c.(x,y) be positive if signals x and y do not share frequencies?

The solution to the self-proposed problem was correct, as it stands.


I hope I won't repeat what you already know...
It's not the DFT/FFT algorithm "who" views (or "assumes", or "interprets") it being periodic, rather it's the "mathematician's mindset" which assumed periodic signals at the input of the algorithm, and designed the algorithm that way.
It's exactly on par with saying that "Newton's laws assume that space & time are absolute..."; laws of mechanics do not have a "brain" to assume such a thing, rather Newton himself made that assumption that the time & space were absolute, and then formulated his laws accordingly, which compute dynamical quantities based on that assumption.
Today we know (believe) that those assumptions are not valid (Special Relativity), yet if you use Newton's laws to compute high-speed phenomena, you may end up objects having speeds larger than speed of light; this inconsistency is due to fact that the "assumptions of the law" are not validated. Here the assumptions are reflected into the mathematical structure of the formulation of the law.
Moving to DSP; in essence, Discrete Fourier Transform (DFT) is by definition one period of the Discrete Fourier Series (DFS); and very naturally DFS works with periodic sequences. I don't know who defined it that way, but it's the case.
So, you cannot compute the DFS of a single rectangular pulse, for example. If you try to find the DFS/DFT coefficients of such a nonperiodic pulse, internally, the mathematical algorithm treats it as if it were periodically extended, and computes the coefficients that way; by fitting periodic sine waves to the periodically extending waveform. Then, when you want to synthesize back the input waveform with the given DFT/DFS coefficients, and evaluating the formula for all n, then you will get not a single pulse but a periodic extension of that pulse. So the zeros outside of the pulse duration are not naturally generated by the inverse DFS/DFT (synthesis) formula.
You must externally force your output to be zero for those out of duration samples. And that's typically what's done with the DFT block length definition.


Check out Julius O. Smith III's write up.

There is a Hilbert transform relationship between the magnitude response, $G(\omega)$, and the phase response, $\theta(\omega),$ of the associated minimum phase filter.

If

$$
H(\omega) = G(\omega) \exp(\jmath \theta(\omega))
$$
then
$$
\ln( H(\omega) ) = \ln(G(\omega)) + j\theta(\omega)
$$
and
$$
\theta(\omega) = - {\cal H}\left[ \ln(G(\omega)) \right]
$$
where ${\cal H}$ is the Hilbert transform.

Also check out Wikipedia's write up.


DFT sample point k < N has negative frequency

The DTFT frequency $\omega_k$ corresponding to the $N$-point DFT index $k$ is given by (frequency sampling relation)

$$ \omega_k = \frac{2 \pi}{N} k ~~~,~~~ k = 0,1,...,N-1.$$

The confusion might arise with the interpretation that the, on the continuous-time domain these frequencies correspond to:

$$  f_k = \frac{f_s}{N} k ~~~,~~~ k = 0,1,...,N-1 ?$$

But no. Because, after about $k = N/2$ which corresponds to $\omega = \pi$ and $f = f_s/2$, you will be sampling the negative portion of the shifted spectrum: $X(e^{j (\omega-2\pi)})$, which corresponds to negative discrete-time frequencies of $X(e^{j \omega})$ due to the periodicity of DTFT $X(e^{j\omega})$.

A more relevant (eventually the same) interpretation of the CTFT, DTFT, DFT frequency relation can be given by the following.

Let $x_c(t)$ be bandlimited to $-B \leq f <B$ and sampled at $F_s = 2B$ samples per second. The base period of the spectrum of $x[n]$ is between $-\pi \leq \omega < \pi$. Note that only the left member is included and the right member is discluded.

Now assume an $N$-point DFT $X[k]$ of $x[n]$ is computed. According to the frequency sampling interpretation of DFT, samples of $X[k]$ are given by :

$$ X[k] = X(e^{j \frac{2\pi}{N} k}) ~~~~ , k = ...-3,-2,-1,0,1,2,3... $$

Note very carefully, that eventhough most of the time we restrict our attention of the DFT index $k$ into the range of $0 \leq k < N$, indeed, the relation is valid for all $k$ which is a manifestation of DFT $X[k]$ and DTFT $X(e^{j \omega})$ being periodic by $N$ and $2\pi$.

Hence, instead of taking the range of $k$ in $[0,N-1]$, we can also take it in the range: $[-N/2 , N/2-1]$ (assuming $N$ even for the moment).

So, if you take $k = -N/2,-N/2+1,...,-1,0,1,2,...,N/2-1$ as the range of $k$ for DFT $X[k]$, then the mapping into the negative frequencies become obvious.

Note that, conventional FFT software always take the range of $k$ in $[0,N-1]$. Therefore, one has to relate them too which is simple: for $k$ in $[-N/2,-1]$ it corresponds to $k'$ in $[N/2,N-1]$ in the conventional DFT range, and for $k$ in $[0,N/2-1]$ it corresponds to the $k'$ in $[0,N/2-1]$ range in the conventional range.


As a pragmatic applied mathematician who has job to do, say, rather than to prove something via theorem or induction, my first approach would be to reconstruct the time series without IMF C2 and C3 and see what happens. I think that's the equivalent of using noise filters in the FFT or deleting outliers in a data set.

As an aside, I'm impressed that you could pick up and utilize the HHT so quickly. Kudos to you.


Trying to find an intuitive understanding of OQAM/FBMC

Since FBMC doesn't depend on the channel acting as circular convolution on the transmit signal like OFDM does, it makes no sense to use a cyclic prefix. (By the way, not all OFDM is cyclic-prefix OFDM, some OFDM-schemes simply use a silent guard interval, but these are the exception.)
To battle ISI without usage of a guard interval, one has to fall back to doing an equalizer per sub-channel in FBMC (something that OFDM reduces to multiplication with the inverse of the complex channel coefficient, but bought at the expense of rate loss due to necessity of the guard interval/cyclic prefix).
However, that equalizer is typically relatively benign in complexity, since the FBMC subchannels are designed to not see very many symbols in delay spread.
That's exactly what the PHYDYAS-style FBMC does: it uses the techniques used in modern single-carrier communications per carrier, and uses OQAM to combat inter-carrier-interference, whilst employing FFT methods for the multicarrier modulation.


Any way to reverse engineer professional audio algorithms or those in e.g. rack units?

Not really.

Anything that's a linear and time invariant is of course easy to characterize but most pro audio stuff doesn't fall into this category. It's typically dynamic, time variant, non-linear, signal dependent, etc. These types of algorithms are extremely difficult to quantify if you don't know the internal structure and the parameter space.

Even if you get a full characterization, you may still not learn something about the way it's done. A time-invariant reverb, for example, can be measured through it's impulse response. However, you can't really tell how it's generated (unless it's a convolutional reverb in the first place).


stretch factor in SOLA time stretch methods

No, the start index doesn't change the stretch factor because the overlapping windows do not have to absolutely fit the available length of the recording.

You can still effectively pad the original signal with zeros by returning a zero value for any request of samples outside the time interval of the original recording.


Special case: Band-limited in time domain and band-limited in frequency domain

The condition is that the function $f(t)=0$, if you are talking about continuous time signals and frequencies (some of your formulations are confusing to me). In other words, a (nonzero) time-limited signal cannot be also band-limited. In other words, a function and its (continuous) Fourier transform cannot both have finite support.

The core of one proof (using complex analysis) relies on Paley-Wiener theorems: the Fourier transform of a time-limited function would be an entire function (a complex-valued function, holomorphic over the whole complex plane), that would vanish on an open interval, and hence be zero (almost) everywhere.

You can find details in Fourier transforms of compactly supported functions, and the paper (for instance) The uncertainty principle for Fourier transforms on the real line. This is a reason for the development of time-frequency analysis.


Replace $x$ by $x - \tau$ . So option b i.e $3x-3\tau$.


Buffering block and filter length in real-time processing

Your data block can be as short as 1 sample (although longer is usually more efficient), as long as you save each full convolution response vector to overlap-add/sum into all the following output blocks.  Look up overlap-add/save fast convolution see how this can be done using an FFT/IFFT plus some add/save buffering.


SDR : Signal divided between the edges of the receiving window

That's aliasing in its most classic form.

The anti-alias filter in your receiver isn't perfect, so that frequencies slightly outside of what you want to sample appear on the other end.


Autocorrelation - Understanding reduced correlation at periodic time shifts using np.correlate (versus statistical autocorrelation)?

The difference here is that statistical autocorrelation assumes a stationary power signal (so basically an infinite periodic signal) and does a normalisation to $[-1,1]$ and that autocorrelation just takes the finite signal as it is and does the classical convolution of the signal with itself, leading to just one point of perfect overlap and decreasing amplitudes before and after.


Alternative for pcshow() function in matlab

pcshow() was implemented in 2015b and in R2015a the equivalent function is called showPointCloud.


This Octave / Matlab code gives you a 2nd ord notch filter at $\omega_n = \pi/6$

r = 0.99;          % notch radius (closer to 1 stiffer)
wn = pi/6;         % notch radian frequency...

% Create the 2nd order NOTCH filter coefficients b() and a()
b = r*conv([1,  -exp(j*wn)],[1,  -exp(-j*wn)]);
a = conv([1,  -r*exp(j*wn)],[1,  -r*exp(-j*wn)]);

figure,freqz(b,a,2048);

with the following result:


Why Sparse Priors Like Total Variation Opts to Concentrate Derivatives at a Small Number of Pixels?

I will divide my answer into 3 sections.

The Distribution of the Derivative of Images

Take a real world image, any image.
Apply the derivative operator on it (Namely apply the kernel $ \left[ 1, -1 \right] $ on it.
Display the histogram of the filtered image.

I took this image:

The histogram I got is this:

This distribution is very similar to Laplace Distribution.
Clearly, this distribution is "Sparse" namely most of the values are 0 (Or close) and very few are different.
In real world we assume most of the values which are small yet not zero are actually due to the noise.

Total Variation Optimization Problem

Look at the Optimization Function:

$$ \hat{x} = \arg \min_{x} {\| Hx -y \|}^{2} + \lambda {\| G x \|}_{1} $$

Where $ H $ is the Blur Operator and $ G $ is the Derivative Operator.

Now, you can look at it in the Sparse sense or you can look at it as the MAP solution given a Laplace prior for the Gradient.

As written above, it fits very well to real world images.

Sparse Optimization

Look at the Optimization Function:

$$ \hat{x} = \arg \min_{x} {\| Hx -y \|}^{2} + \lambda {\| G x \|}_{0} $$

Where $ H $ is the Blur Operator and $ G $ is the Derivative Operator.

This optimization problem clearly promote Sparse solution (With respect to the Derivative).
Yet, this is a very hard problem to solve.
Hence it was shown that under some circumstances the solution of the problem with $ {\ell}_{1} $ coincide with this solution.

Moreover, by looking at the Unit Sphere of different norms (And the Pseudo Norm $ {\ell}_{o} $, or more correctly Cardinality Function) you can see why the lower the norm the Sparse Solution it promotes.

All in all you have many point of view on the same problem.


2D convolution of image with filter as successive 1D convolutions

If a 2D $K_2$ filter kernel is of rank $0$ or $1$, it can be written as a separable product of $2$ 1D kernels   $K_1^r$ and $K_1^c$ on rows and columns. As such, it can implemented by 1D convolutions, as long as one properly reshape the 2D matrices into 1D ones, and take care about "out-of-range" values, to avoid wrap-around. For instance, you can pad in every direction by the size of the filter, and make sure the convolution does not add unwanted information.

Assuming that you know that you have a separable 2D filter, the following code does the job. A one-liner would be:

xRowFull = reshape(conv(reshape(reshape( conv(x(:),s1,'same'),nRow,nCol)',nRow*nCol,1),s2,'same'),nRow,nCol)';

And the code is:

% https://dsp.stackexchange.com/questions/62115/2d-convolution-of-image-with-filter-as-successive-1d-convolutions
%% Initialization
clear all
nRow = 16;
nCol = 16;
HalfSizeCentralImageKernel = 1;
x = zeros(nRow,nCol);
x(nRow/2-HalfSizeCentralImageKernel:nRow/2+HalfSizeCentralImageKernel,nCol/2-HalfSizeCentralImageKernel:nCol/2+HalfSizeCentralImageKernel)=rand(2*HalfSizeCentralImageKernel+1);

%% Original 2D version
s1=[1,0,-1]';
s2=[1 2 1];
y = s1*s2;

%% Step by step 2x1D version
xRowFlat1 = x(:);
xRowFlat1FiltCol = conv(xRowFlat1,s1,'same');
xRowFlat2 = (reshape(xRowFlat1FiltCol,nRow,nCol))';
xRowFlat2 = xRowFlat2(:);
xRowFlat2FiltRowFlat = conv(xRowFlat2,s2,'same');
xRowFlatFilt2Row = reshape(xRowFlat2FiltRowFlat,nRow,nCol)';

%% Compact vectorized 1D version
xRowFull = reshape(conv(reshape(reshape( conv(x(:),s1,'same'),nRow,nCol)',nRow*nCol,1),s2,'same'),nRow,nCol)';

%% Display
figure(1);
imagesc(x);

figure(2);
subplot(1,3,1)
imagesc([conv2(x,y,'same')]); xlabel('Original')
subplot(1,3,2)
imagesc(xRowFlatFilt2Row); xlabel('Separable, step by step')
subplot(1,3,3)
imagesc(xRowFull); xlabel('Separable, one-liner')

diff1=conv2(x,y,'same')-conv2(conv2(x,s1,'same'),s2,'same');
disp(['Max error 1: ',num2str(max(abs(diff1(:))))]);

diff2=conv2(x,y,'same')-xRowFlatFilt2Row;
disp(['Max error 2: ',num2str(max(abs(diff2(:))))]);

[First answer]

Here is a crude  Matlab code. Can you test it, and if OK, i'll send a one-liner (if I can).

nRow = 8;
nCol = 8;
HalfSizeCentralKernel = 1;
x = zeros(nRow,nCol);
x(nRow/2-HalfSizeCentralKernel:nRow/2+HalfSizeCentralKernel,nCol/2-HalfSizeCentralKernel:nCol/2+HalfSizeCentralKernel)=rand(2*HalfSizeCentralKernel+1);
figure(1);
imagesc(x);

% 2D version
s1=[1,0,-1]';
s2=[1 2 1];
y = s1*s2;
diff1=conv2(x,y,'same')-conv2(conv2(x,s1,'same'),s2,'same');
disp(['Max error 1: ',num2str(max(abs(diff1(:))))]);

% 1D version
xRowFlat1 = x(:);
xRowFlat1FiltCol = conv(xRowFlat1,s1,'same');
xRowFlat2 = (reshape(xRowFlat1FiltCol,nRow,nCol))';
xRowFlat2 = xRowFlat2(:);
xRowFlat2FiltRow = conv(xRowFlat2,s2,'same');
xRowFlatFilt2Row = reshape(xRowFlat2FiltRow,nRow,nCol)';

figure(2);
subplot(1,2,1)
imagesc([conv2(x,y,'same')])
subplot(1,2,2)
imagesc(xRowFlatFilt2Row)


Differences between OpenCV Canny and MatLab Canny?

As suggested above, the Matlab Canny edge detector calculates the gradient using a "derivative of a Gaussian filter" (as stated in the documentation). In other words, Matlab does a Gaussian blur of the image and then finds the gradient of that smoothed image... all using a single fancy filter. [If you want to know the details, just type in edit edge as Andrey suggested, and then scroll down to the smoothGradient() function.]

The blurring operation significantly reduces the amount of noise present in the image, eliminating many spurious edges and leaving behind the good stuff.

Unfortunately, the OpenCV Canny function doesn't let you change the filter kernel it uses via the function parameters. However. You can generate the same results by first blurring the input image, and then passing this blurred image into the Canny function.

This significantly cleans up the resulting edge map. To blur the input image, I personally use OpenCV's GaussianBlur() function with sigmaX=2. This mimics the default sigma in Matlab. The best blurring kernel size can vary from case to case, but in Matlab it is calculated using filterLength = 8*ceil(sigma);, so for a sigma of 2 that would mean a kernel size of (16,16)

Since both the Gaussian Blur and Sobel filters are linear, passing a blurred input image to the OpenCV Canny() function is mathematically equivalent to what Matlab does because of the principle of superposition, as demonstrated in this pseudocode (note: * is the convolution operator):

// The Matlab method: the sobel and blur operations are combined into
// a single filter, and that filter is then convolved with the image
matlabFancyFilter = (sobel * blur);
gradient = matlabFancyFilter * image;

// Equivalent method: image is first convolved with the blur filter, and
// then convolved with the sobel filter.
gradient = sobel * (blur * image); // image is filtered twice

There is an OpenCV Canny tutorial here that demonstrates how to do this using C++. I'm a python guy, so here is what I do:

smoothedInput = cv2.GaussianBlur(image, (7,7), 2);
edges = cv2.Canny(smoothedInput, 25, 50);


My attempt:

Consider the following diagram showing the beam with length L, sensors A and B on each side of the beam and an arbitrary position x representing the location of the discontinuity. First consider the case of a signal source emanating from x equally toward A and B, with zero loss in the beam. The overall distance of the beam itself is such that the delay from A to B does not cause the phase of a signal at it's sampling rate to exceed $2\pi$ (otherwise alias solutions will exist). Smaller phases at this frequency would be allowable, but this will allow for the best precision.

Thus similar to the transfer function of any cable the transfer function from x to A over the normalized radian frequency of $0$ to $2\pi$ has a magnitude of $1$ and a linear phase that extends from $0$ to $-\frac{2\pi}{x/L}$. Similarly the transfer function from x to B would have a magnitude of $1$ and a phase that extends from $0$ to $-\frac{2\pi}{(L-x)/L}$

These are given as:

$$H_{xA}(\omega) = e^{-j\omega L/x}$$

$$H_{xB}(\omega) = e^{-j\omega L/(L-x)}$$

And the ratio of these two is given as:

$$H(\omega) = H_{xA}(\omega)/H_{xB}(\omega) = e^{-j\omega L (L-x)/x}$$

With $\omega$ as the normalized radian frequency in units of radians/sample ($\omega$ extends from $0$ to $2\pi$ with $2\pi$ representing the sampling rate assuming the beam is sized such that the phase delay of the sampling rate was $2\pi$ (or this can be modified accordingly for shorter delays).

Consider a single tone given as $e^{j \omega_n}$ that is originating at $x$. The signal received at A compared to the signal received at B is given by $H(\omega_n)$ to be $e^{-j\omega_n L (L-x)/x}$.

This would appear at A as $e^{-j\omega_n L/x}$ and at B as $e^{-j\omega_n L/(L-x)}$.

Specifically if we followed the guidelines above on length versus sampling rate, the phase of this comparison signal is unique to $\omega_n$ and $x$ for a given $L$ as $\omega_n L/(L-x)$.

Given that and to the extent we can model the discontinuity as a source that emanates equally in all directions, I would recommend the following approach:

Measure A as $a(t)$ and compute the DFT to be $A(\omega)$

Measure B as $b(t)$ and compute the DFT to be $B(\omega)$

Compute the comparison signal $R(\omega) = A(\omega)/ B(\omega)$

Compare the phase of $R\omega$ to the phase of $H\omega$ and solve for x for each bin:

$$arg(R(\omega)) = \omega L/(L-x)$$

$$x(\omega) = L\bigg(  1 - \frac{\omega} {arg(R(\omega))} \bigg)$$

Where $\omega = 2\pi k/N$ for DFT bins $k = 0$ to $N-1$.

This approach can be extended for lengths with phase shifts less than $2\pi$ at $f_s$ as well as actual transmission loss factors (since lossless was assumed here).

Note that two or more discontinuities that react to the same frequency will not be resolvable: their relative phase contributions will add and so will therefore appear to be coming from another location (such as images in a reflection).

We can calculate this for your specific example:

$$f_s = 9000 \text{ Hz}$$
$$L = 74 \text{ m}$$
$$v = 12000 \text{ m/s}$$

The wavelength in meters of the sampling frequency is $v/f_s$ = $12000/9000 = 1.33 \text{ m}$

This means a the maximum length of the beam is 1.33 meters for a unique solution over the full sampling rate given this sampling rate given. However all frequencies up to point where their delay would be $2\pi$ would be resolvable and this frequency would be:

$$f_{max} = \frac{v}{L} = \frac{12000}{74} \approx 162.2 \text{ Hz}$$

Given this, there is no reason to sample higher than twice this frequency but it doesn't hurt to do so. We would expect to see repeating images of everything from DC to 162.2 Hz at the higher frequencies.


Ok you are mixing the order of the transforms. When you decompose the argument $3t-6$ into $3\cdot(t-2)$ , then you should first scale $x(t)$ by $3$ and then shift the scaled result by $2$, that's where you have chosen the wrong order. Let's show the correct order:

Let $x(t)$ be the original signal whose CTFT is $X(\omega)$, then define the following signals in two steps to get the final signal $z(t)= x(3t-6)$:

$$ y(t) = x(3t) \tag{1}$$ and
$$z(t) = y(t-2) \implies z(t) =  x(3(t-2)) = x(3t-6) \tag{2}$$

Notice how we replaced $y(t-2)$ by using its definition from the eq. (1).

Now applying CTFT to the above rows yields:

$$ Y(\omega) = \frac{1}{3} X(\omega/3) \tag{3} $$ and
$$ Z(\omega) = Y(\omega) e^{-j2\omega} \tag{4} $$

finally insert $Y(\omega)$ from eq. (3) into eq. (4) to get:

$$ \boxed{Z(\omega) = Y(\omega) e^{-j2\omega} = \frac{1}{3} X(\omega/3) e^{-j2\omega} }$$

By the way, you could have also used the alternate approach: first shift $x(t)$ by $6$: $$y(t) = x(t-6) \implies Y(\omega) = X(\omega)e^{-j6\omega}$$ and then scale $y(t)$ by $3$ to obtain $$z(t) = y(3t) = x(3t-6) \implies Z(\omega) = \frac{1}{3}Y(\omega/3) = \frac{1}{3}X(\omega/3)e^{-j2\omega}$$


Cross correlation should work.  I think the problem is the waveform that you are using.  A square wave has bad auto-correlation properties.  If it is a periodic square wave it will have multiple peaks.  It sounds like you are just using a single pulse which is better, but it will still have a gradual roll-off which is a problem.  Instead, use a Barker code, which will give you a nice sharp peak and very low side peaks.  I would BPSK modulate the barker code and then cross correlate what you receive with the "ideal" signal.


List of Fourier coefficientes to simulate musical instruments

look up Additive Synthesis somewhere.  those Fourier coefficients are not constants.  the musical note is not normally a perfectly periodic function in time.  it is what i like to call "quasi-periodic".  so this is not exactly true (given some period $P$):

$$ x(t) = x(t+P) \quad \text{for all } t $$

but in the neighborhood of some given time $t_0$ this is approximately true

$$ x(t) \approx x(t+P(t_0)) \quad \text{for } t \approx t_0 $$

so the exact Fourier series for $x(t)$

$$ x(t) = \sum\limits_{n=-\infty}^{+\infty} c_n \ e^{j 2 \pi n f_0 t} \quad f_0 \triangleq \frac{1}{P} $$

is not precisely true, but more like

$$ x(t) = \sum\limits_{n=-\infty}^{+\infty} c_n(t) \ e^{j \theta_n(t)} $$

where

$$ \theta_n(t) = \int\limits_{0}^{t} 2 \pi f_n(u) du + \theta_n(0) \quad f_n(t) \triangleq \frac{n}{P(t)} $$

so what you're looking for is a collection of envelopes, rather than a collection of constant coefficients, for some musical instrument.

in a previous life, living on an older computer (that might not power up now) and maybe even on some 3.5 inch floppy disks, i have some of that additive synth data.  actually, i am not sure where to get it all digested.  within the context of Wavetable Synthesis, i know how i would go about extracting it from "samples" of notes (like what you would get from a sample library like from Synthogy ).


Designing a root-locus compensator to attain design requirements

I don't have access to matlab, but I have played with this in R (see code below) and can confirm that if you're trying to transform your $G(s)$ into:
$$
H(s) = \frac{K G(s)} {1 + KG(s)}
$$
using only a scalar $K$, then you cannot get to the target pole locations (based on the values of $\xi$ and $\omega_n$).

In the diagram:

The black circles are the starting root (pole) positions.
The green circles are the target pole positions.
The red circles are the root locus for varying values of $K$.
The large black circle is the unit circle.

R Code Below

#30341

library("signal")

xi <- 0.5
omega_n <- 0.3

target_den <- c(1, 2*xi*omega_n, omega_n*omega_n)

target_roots <- roots(target_den)

num <- c(0, 0.5, 0.5)
den <- c(1, -2, 1)

start_roots <- roots(den)

par(pty="s")
plot(Re(start_roots),Im(start_roots),xlim=c(-5,5),ylim=c(-5,5), col="black", lwd=5)
points(Re(target_roots), Im(target_roots), col="green", lwd = 3)
for (k in seq(-100.0, 100.0, 0.1))
{
  rts <- roots(k*num + den)
  points(Re(rts), Im(rts), col="red")
}

theta <- seq(-1,1,.01)*pi
lines(sin(theta), cos(theta), col="black")


The problem you are most likely having is that you need to search for files in the correct format.  DICOM is what you need to look for.
A few hits:

One example
A library from a DICOM viewer software producer
A library of images from Siemens

You will also need a DICOM viewer.
Image J  is free (non-commercial, no restrictions on use) and can be used under Windows, Mac, and Linux - it does require Java.
There are other viewers available, free as well as commercial.  Just search for DICOM viewers.


When generating the sine wave, the code uses a sampling frequency to establish discrete points in the time domain. It is not continuous. This means the wave is broken up into bins that are exactly 1/fs in size/width.

So if you want a true integral of the wave, you have to take the average density of all the points in a bin by dividing by the sampling frequency.

Hope that helped.


The Number of Sine and Cosine Waves in an $ N $ Point DFT

There is no contradiction.

Both cases you have information from $ -\frac{ Fs }{ 2 } $ to $ \frac{ Fs }{ 2 } $.

The representation using Cosine ans Sine functions is using Real Functions.
Those functions are symmetric as Real Valued functions.
Namely if you a coefficient for a Sine in a given frequency it is for sure have energy on the negative frequency as well (See the DFT / Fourier Transform of Sine and Cosine).
Namely using those functions all needed it to go from $ 0 $ to $ \frac{ Fs }{ 2 } $ as the negative value is implicit.

On the other end, for the Complex Exponent the energy is one sided.
Hence you must go from $ -\frac{ Fs }{ 2 } $ to $ \frac{ Fs }{ 2 } $.


Yes, it is possible to implement a Hilbert transform in the analog domain, although it is difficult to get the desired phase shift over a significant bandwidth when you take component uncertainties and imperfections into account. There is some discussion of the topic in this document.

I'm not sure why you think you would use an integrator instead of a delay. The analog counterpart of a discrete-time delay is a continuous-time delay. This can be implemented in various ways, although as noted before, achieving precise delays in analog electronics is quite difficult. Some methods that I've seen before for inducing delays include introducing capacitive delays (e.g. with an RC filter) or using surface acoustic wave (SAW) devices.


How to find phase difference between 2 signals, one has an unknown function

I don't fully understand your question, or what your code is supposed to do. However, this line is wrong:

fs =  14*frequency;

You should use the same sampling frequency as was used to acquire the data. Assuming uniform sampling, you can do something like

dt = t(2) - t(1);
fs = 1/dt;

Also, be aware that estimating the actual phase from the FFT is not straigthforward. Here is a quick and dirty example showing one way to do it. Start with this code:

fs = 1000;
t = 0:1/fs:0.1;
f = 100;    % signal frequency
phase = 0;  % signal phase
y = cos(2*pi*f*t+phase);
Y = fftshift(fft(y,2^nextpow2(length(y))));
Ym = abs(Y);
[m, idx] = max(Ym);
Ym(idx) = 0;
[m, idx] = max(Ym);  % index of f
Ya = angle(Y);
ph_e = Ya(idx);
fprintf('Actual phase is: %1.5f\n', phase);
fprintf('Estimated phase is: %1.5f\n', ph_e);

If you run this code you'll get:

Actual phase is: 0.00000
Estimated phase is: -0.49087

The phase estimate depends on many factors, including the start and end times of your time vector. The thing is that, in this particular setup, the FFT interprets a phase of zero as -0.49087. Then, modify the third to last line in the script:

ph_e = Ya(idx) + 0.49087;

Now the script will properly identify the phase. Try it by setting the phase in line 4 to, say, pi/3, then run the script:

Actual phase is: 1.04720
Estimated phase is: 1.03347

Which is pretty good accuracy.

I hope this helps you to get going.


I point you to this https://stackoverflow.com/questions/12723699/changing-image-bit-depth-using-matlab it can be done. This particular post is only for png, but if I'm not mistaken all image types have the depth parameter so they can likely all be changed in a similar manner

EDIT

For your particular problem, since your images are already grayscale you would do something like

%change this as needed
desired_bit_depth = 2;

my_pixel_depth = 2 ^ desired_bit_depth;
%converts grayscale to an indexed image at the appropriate depth
[ind_im, reduced_colormap ] = gray2ind(my_gray_image, my_pixel_depth)

%displays your new image
figure()
subplot(1,2,1);imshow(my_gray_image); title('original')
subplot(1,2,2);imshow(ind_im,reduced_colormap); title(sprintf('reduced to %d bits',desired_bit_depth))

% save indexed png
imwrite(ind_im, reduced_colormap, 'test.png', 'bitdepth', desired_bit_depth);

EDIT #2: While reducing the "bits" to display the image, the "colormap" is made by taking the average of each range
I kept the original answer because it still works and may suit someone elses needs in the future.

For this solution, rather than using evenly spaced values for our "colormap' we actually use the gray value that is the average over the specified range

ex

to convert to 1 bit image we have two ranges 0,128 and 128,256

replace all values between 0,128 with the average value between 0,128

replace all values between 128,256 with the average value between 128,256

function result_im = ChangeBitDepthGrayImage(gray_im, desired_bit_depth)
    if (desired_bit_depth < 1)
        disp('converting to binary(1 bit) image');
        desired_bit_depth = 1;
    end

    if (desired_bit_depth > 8)
        disp('converting to 8 bit image');
        desired_bit_depth = 8;
    end

    %assuming we start with 8 bit image 256 levels
    num_levels = 2 ^ desired_bit_depth;

    %figures out how big each range should be, we use +1 because if we
    %divide the data into N levels, there should be N+1 boundaries
    limits = linspace(0,256,num_levels + 1);

    result_im = uint8(zeros(size(gray_im)));

    for i = 1:num_levels
        lower_lim = limits(i);
        upper_lim = limits(i+1);

        %creates a binary mask of values between the limits, the output is
        %0 or 1, but we need to make it uint8 for the next step
        temp_mask = uint8((gray_im >= lower_lim) & (gray_im < upper_lim));

        %multiplies image by mask, this isolates only pixels in the given
        %range
        image_only_in_range = temp_mask .* gray_im;

        %finds the mean of that small part of the image. this weird notation is
        %taking the average of nonzero elements
        avg_val_for_range = round(mean(image_only_in_range(image_only_in_range~=0)));

        %replaces all pixels in that range with the average val
        result_im = result_im +(avg_val_for_range * temp_mask);
    end

    %i just picked some random figure
    figure(32)
    subplot(1,2,1);imshow(gray_im);title('original image');
    subplot(1,2,2);imshow(result_im);title(sprintf('modified to %dbit image',desired_bit_depth));
end


signed 2's complement of a negative number for Integer = 3 and Fraction = 9

If your routine works correctly for converting positive numbers, you could just use it to convert a negative number by first adding the value of the sign bit (MSB) to the negative number (which must make it positive if it is not out of range), then convert the resulting positive number, and finally add the sign bit.
Example: convert the number $-0.7$ to Q3.9 format. The MSB has a value of $-2^2=-4$, so we first convert the number $-0.7+4=3.3$, which results in
011.010011010 ($= 3.30078125$)
Now you just add the sign bit and you obtain the representation of the original number:
111.010011010 ($=-0.69921875$)


Convex Optimization in Signal and Image Processing

There's a whole area of signal processing dedicated to optimal filtering. In pretty much every case I've seen the filtering problem is formulated with a convex cost function.

Here's a freely available book on the subject - Sophocles J. Orfanidis - Optimum
Signal Processing.


Obtaining the FFT of a data set collected through an accelerometer

From you code it appears you do a low pass butterworth filter:

[b1,a1]= butter(order_of_butterworth_filter,wn,'low');
accelerationd2_filtered = filtfilt(b1,a1,accelerationd2);

Then another low pass moving average filter which you minus from the previous.

accelerationd2_smoothned = smooth(accelerationd2_filtered,36.5);
accelerationd2_DC_shifted = accelerationd2_filtered - accelerationd2_smoothned;

Yes, this second part will remove DC, but I think this is the first thing you should do. That is, use a HIGH pass butterworth filter (with a few more orders perhaps)

cut_off_frequency = 5;
wn = cut_off_frequency/sampling_frequency;
order_of_butterworth_filter = 4;
[b1,a1]= butter(order_of_butterworth_filter,wn,'high');
accelerationd2_filtered = filtfilt(b1,a1,accelerationd2);

Then apply the FFT. If there is two much high frequency noise, you can apply a low pass butterworth filter to the above output.

As for you analysis, the FFT uses the entire signal to compute. So any frequency magnitudes represent the average over the whole time series. That is, there is no time information in there. Since you have a large sample of animal data, I imagine you would want to see how things change over time, so I suggest instead using the STFT or some other wavelet transform.

In MATLAB you might get good results out of the box using spectrogram, no pre-filtering required.

window = 1000;
overlap = round(0.9*window);
s = spectrogram(accelerationd2,window,noverlap)


Why $y[n] = x[-n]$ is not time-invariant?

A time-invariant system is one that, when you shift the input signal, the output is shifted by the same amount.

A system that reverses the signal cannot be time-invariant because when you shift the input, the output is shifted the other way. $k$ and $-k$ are not the same amount.

$$
y[n-k] = x[k-n] = x[-n \mathbf{+} k]
$$


Multipath resistance of spread spectrum systems

The key idea in spread-spectrum communications is to use much more bandwidth than the minimum required for the signal being sent.

Multipath effects occur when the sent signal $s(t)$ and a delayed, attenuated version of it interfere at the receiver:
$$
r(t) = s(t) + \alpha \cdot s(t-\tau)
$$

If $s(t)$ is really narrow band (a sinusoid), and $\alpha = 1$, then
$$
r(t) = \sin(\omega t) + \sin(\omega(t+\tau))
$$
and if $\omega\tau = \pi$ you can get perfect nulling:
$$
r(t) = \sin(\omega t) + \sin(\omega t + \pi) = 0
$$

If, however, $s(t)$ is not a simple sinusoid, then perfect nulling can be eliminated:
$$
r(t) = [\sin(\omega t) + \sin(\omega_2 t) ] + [\sin(\omega(t+\tau)) + \sin(\omega_2 (t+\tau)) ]
$$
even if  $\omega\tau = \pi$, if $\omega_2$ is chosen well, then $\omega_2\tau \ne \pi$ so there is still some signal to receive.


Remove noise on edge (compression deffects)

The Total Variation (TV, a.k.a. Rudin-Osher-Fatemi functional) denoising wold help here.
You can apply it on every color layer or look for a vectorized version if th eresults ar enot good enough for you.

Good references of tsable and easy to implememnt algorithms for TV denoising are for example:

Chambolle's algorithm
the Chambolle-Pock primal dual algorithm

Both references have good performance and are easy to implement from the papers.


The way to measure chromatic tuners' precision

Precision and accuracy are 2 different things.

1 cent is 1/100th of the frequency difference between 2 adjacent semitone pitches in equal temperament tuning.  The frequency ratio is exp(ln(2)/1200) or roughly 1.00059

Precision of 1 cent might mean the tuner reports frequency error from the nearest semitone in 1 cent units, whether accurate or not.

To test accuracy, you might want to synthesize test signals that "sound like" your desired instrument (in software, or using a calibrated programmable synthesizer).  Feed those test signal to your tuner (or your tuner's pitch estimation algorithm in a software environment as unit testing).  If the results matches the synthesized pitch to within 0.5 cents, then your accuracy is 1 cent at that frequency.  Vary the test signal (perhaps randomized) to different semitone notes and cent (or cent fraction) offsets over your tuners entire range to find the worst case to evaluate (get a statistical estimate of) your tuners accuracy.

You can add Gaussian random white noise or recorded typical background room noises to your generated test signals to test measurement precision (repeatability) at realistic S/N ratios as well.

One "gold standard" is to measure the WWV audio frequencies (500 and 600 Hz), and use those to calibrate your synth and test your tuner.


The expectation in power spectral density

I think a better definition of the power spectrum is the following:

  The power spectrum of $x(t)$ is the Fourier transform of the autocorrelation function of $x(t)$, where $x(t)$ can be either a deterministic power signal, or a wide-sense stationary (WSS) random process. The definition of the autocorrelation function depends on the model for $x(t)$.

If $x(t)$ is modeled as a WSS random process, then the autocorrelation function is defined by

$$R_x(\tau)=E\big\{x^*(t)x(t+\tau)\big\}\tag{1}$$

For deterministic power signals, the autocorrelation function is given by

$$R_x(\tau)=\lim_{T\to\infty}\frac{1}{T}\int_{-T/2}^{T/2}x^*(t)x(t+\tau)dt\tag{2}$$

In this answer it is shown that the following definition of the power spectrum for a WSS random process $x(t)$

$$S_x(\omega)=\lim_{T\rightarrow\infty}E\left\{ \frac{1}{T}\left| \int_{-T/2}^{T/2}x(t)e^{-j\omega t}dt \right|^2 \right\}\tag{3}$$

is equivalent to the definition of the power spectrum as the Fourier transform of $(1)$.

For deterministic power signals, the corresponding definition of the power spectrum is

$$S_x(\omega)=\lim_{T\rightarrow\infty}\frac{1}{T}\left| \int_{-T/2}^{T/2}x(t)e^{-j\omega t}dt \right|^2 \tag{4}$$

which can also be shown to be equivalent to the Fourier transform of $(2)$.

The definitions of autocorrelation and power spectrum of deterministic power signals is described in Chapter 12 of

Papoulis, A., The Fourier Integral and its Applications, McGraw Hill, 1962.

A good reference on random processes and the corresponding definitions of autocorrelation and power spectra is

Papoulis, A. and S.U. Pillai, Probability, random
variables, and stochastic processes, Boston: McGraw-Hill, 2002.


Psycho acoustic masking is a VERY complicated subject. The masking thresholds (energy required to have one sound be audible in the presence of another) is a very complex function of

Spectral content and overlap
Spatial presentation (there are many books just on binaural masking alone)
Time shapes of the waveforms
Time alignment and envelopes of the wave forms

Even the more advanced audio codecs (such HE AAC) use very simplified models because a more accurate one is very complicated.

The best choice really depends on your application and requirements, monoaural masking of simultaneous sinusoids is probably the easiest case. You may find some helpful stuff here https://engineering.purdue.edu/~malcolm/interval/1998-010/


A real-valued negative number has a phase of $\pm\pi$, because $e^{\pm j\pi}=-1$. So for all frequencies for which the real-valued spectrum is negative, the corresponding phase angle equals $\pm \pi$:

$$S(j\omega)=\begin{cases}|S(j\omega)|,&&S(j\omega)\ge 0\\
|S(j\omega)|e^{\pm j\pi},&&S(j\omega)<0\end{cases}$$


Determine filter type (HPF or LPF) by its coefficients

I am not sure if my answer will correctly suits your question as it seems a little confusing. Anyway, I understand that you want to eliminate the initial transition of your filter. What I propose is to initialize your filter state so that the output equals the input at $t=0$

I had to solve that problem in the past. What I did was to assume that all delays were initialized to the same value. Doing that, you simplify your problem a lot. I don't have a lot of theory to back this technique, but that seems to work quite good on my side.

Let's start with Direct Form II

If we consider that all delays are set to the same value $K$, we find this output equation.

$$ y = (-a_1k-a_2k+x)b_0+b_1k+b_2k $$
A little algebra :
$$ y = b_0k(-a_1-a_2-...)+b_0x+k(b_1+b_2+...)$$
$$ y = -kb_0\sum_{n=1}^{N_a}{a_n} + b_0x+k\sum_{n=1}^{N_b}{b_n}$$

Now, since we know we want the output to be equal to the input, we can say $y=x$

$$ x = -kb_0\sum_{n=1}^{N_a}{a_n} + b_0x+k\sum_{n=1}^{N_b}{b_n}$$

Then

$$ x-b_0x = -kb_0\sum_{n=1}^{N_a}{a_n}+k\sum_{n=1}^{N_b}{b_n}$$
$$ x(1-b_0) = -kb_0\sum_{n=1}^{N_a}{a_n}+k\sum_{n=1}^{N_b}{b_n}$$
$$ \frac{x(1-b_0)}{-b_0\sum_{n=1}^{N_a}{a_n}+\sum_{n=1}^{N_b}{b_n}} = k $$

So there you go. Initialize all your delays with $K$ and your filter will start with an output equal to the input.

You can do this for other forms as well.

Hope that helps


How can I properly detect a ball moving in a film capture?

First, for removing shadows try to normilize your RGB image prior to converting it to HSV color model. For more information on this regard see :
http://aishack.in/tutorials/normalized-rgb/

In case the ball is always moving, Optical flow would be a good option I guess. Simply calculate optical flow and mark the region with highest flow as the ball.

The other option would be emplying a tracking algorithm, which, I feel MeanShift would be a good option here. However, a tracking algorithm requires you to pinpoint the object of the interest in first few frames. If the object is fixed (it's always the same ball) then tracking would diffinitly work.


What steps come after finding the gradients of the SIFT key points?

First of all, there are two distinct parts to SIFT.  The first part is interest point detection algorithm (aka key-point detection), which finds local extrema of the multi-scale difference-of-gaussians function.

The second part is computing the feature descriptor, which is a vector describing the image patch around each key point. SIFT computes this descriptor by taking a 16x16 block of pixels centered at a key point, dividing it into 4x4 cells, and computing an 8-bin histogram of gradient orientations within each cell.  This results in a 4*4*8=128-element vector, which is the descriptor.

Finally, the SIFT algorithm normalizes the feature descriptor for rotation.  It finds the dominant gradient direction in the patch around the key point and shfits the gradient histograms such that the first orientation bin corresponds to that direction. This makes the descriptors invariant to in-plane rotation.


How can I obtain the response signal for this question?

HINT:

The (zero-state) step response is just the cumulative sum of the impulse response $h[n]$:

$$y_{ZS}[n]=u[n]\sum_{k=0}^nh[k]$$

This follows in a straightforward manner from the convolution of a unit step $u[n]$ with the impulse response:

$$y_{ZS}[n]=(h\star u)[n]$$


Why is the Kalman gain constant when position and velocity change?

The Kalman Gain $K_k$ for a discrete, linear system is computed using the state covariance matrix $P_k$, measurement matrix $H_k$ and measurement noise covariance matrix $R_k$:

$$
K_k = P_kH_k^T(H_kP_kH_k^T+R_k)^{-1}
$$

If the measurement matrix and the measurement noise covariance matrix are constant and the state covariance matrix converges to a steady-state value $P_\infty$, the Kalman Gain also converges to a steady-state value $K_\infty$, regardless of actual measurements or state.


The essence of your problem is one of Voice Activity Detection (VAD). This is a well researched topic and you should investigate the available algorithms - there are many papers available online, here for example

You should apply a VAD algorithm to the signal and replace non-speech areas with silence, you can then look at the specific edge cases which are concerning you, such as animal sounds have being incorrectly considered speech.

For the case of animal sounds you may want to investigate Formant Detection (for example this paper) since as you pointed out there is a distinction between them. You may also want to consider developing a heuristic based on your own empirical data about the differences. You can take many signal features such as the fundamental frequency, the signal power/amplitude, peakedness etc... and weight these appropriately to give you a confidence rating about the likelyhood of it being animal rather than human speech. This is a fuzzy approach, you would need to iterate a lot on the weightings and metrics.


If you want to merely compare the periods of the functions $f_1$ and $f_2$, you could use the following "metric" (it might not technically be one if it disobeys the triangle inequality, but someone else will have to check for that):

$\Delta (f_1, f_2) \equiv \dfrac{\left| T_1 - T_2 \right|}{\max {T_1, T_2}}$

Then the similarity can simply be defined as $s(f_1, f_2) \equiv 1-\Delta(f_1, f_2)$


Yes, they both give you the same SNR gain. I suppose there might be a dynamic range issue, as you say, if a LOT of samples are added, but I think the main reason that the literature speaks of "averaging" instead of "adding" is because for systems that allow a variable integration rate, we can be assured that other parts of the system work the same, like AGC.


You know the geometric series

$$\sum_{n=0}^{\infty}r^n=\frac{1}{1-r},\quad |r|<1\tag{1}$$

Due to convergence of (1) (for $|r|<1$) you can take the derivative with respect to $r$ by taking the element-wise derivatives of the left-hand side. Equating the derivatives of both sides of (1) gives

$$\sum_{n=0}^{\infty}nr^{n-1}=\sum_{n=1}^{\infty}nr^{n-1}=\sum_{n=0}^{\infty}(n+1)r^n=\frac{1}{(1-r)^2},\quad |r|<1\tag{2}$$

With $r=az^{-1}$ you get the following $\mathcal{Z}$-transform relation:

$$(n+1)a^nu[n]\Longleftrightarrow \frac{1}{(1-az^{-1})^2}=\frac{z^2}{(z-a)^2}\tag{3}$$

where $u[n]$ is the unit step sequence.
With the time-shifting property of the $\mathcal{Z}$-transform (see here) you get the following pair:

$$(n-1)a^{n-2}u[n-2]\Longleftrightarrow \frac{1}{(z-a)^2}\tag{4}$$

With (4) you immediately get

$$\frac{5}{(z-2)^2}\Longleftrightarrow 5(n-1)2^{n-2}u[n-2]$$

Of course, all of this is only valid inside the region of convergence, i.e. for $|z|>a=2$.


Take a look on Generalized/constrained Procrustes Problems.

It should be sufficient to update singular values in $\Sigma$.

The $\det(SVD)=-1$ case is discussed in more detail in 1987 paper of Arun, Huang and Blostein: "Least-squares fitting of two 3-D point sets".


I don't think the exam question is stated properly, unless there is additional context that the student is supposed to know.

First, the bandwidth depends on the encoding. The question seems to assume binary encoding, but PCM does not specify this.

Second, the bandwidth depends on the pulse shape. Usually one assumes rectangular unless otherwise specified, but this information should bave been included in the question.

In general, the PCM signal can be written as a pulse train:

$$s_{\text{PCM}}(t) = \sum_{k=-\infty}^\infty a_k p(t-kT_p),$$

where $T_p$ is the pulse rate and $a_k$ is the encoded PCM data. In general, the bandwidth of $s_{\text{PCM}}(t)$ is equal to the bandwidth of $p(t)$.

As an example, let's say we want to transmit the bits $1010$ using a bipolar, non-return to zero binary encoding with $a_k \in \lbrace 1, -1 \rbrace$. The pulse shape is rectangular with duration $T_p = 1/(nf_s)$. Assuming $n$ starts at $n=1$, the PCM signal can be written as:

$$s_{\text{PCM}}(t) = p(t-T_p)-p(t-2T_p)+p(t-3T_p)-p(t-4T_p).$$

The problem is, what is the bandwidth of this signal? Theoretically, it's equal to the bandwidth of the rectangular pulse $p(t)$ with duration $T_p$. However, there are several ways to measure that bandwidth.

For instance, if you use the "absolute bandwidth" definition, then the bandwidth is infinite. There are ten or more engineering definitions of bandwidth, and each gives you a different number. Note, however, that transmitting $s_{\text{PCM}}(t)$ over any channel with less than infinite bandwidth will result in some amount of distortion.

Now let's say that $p(t)$ is a sinc pulse of appropriate width. The expression for the signal doesn't change. Theoretically, now the signal has infinite time duration, but its bandwidth is finite and equal to $1/(2T_p)=nf_s/2$.


Detecting changes between voice and music

This is a well-studied problem, dating back from the mid 90s (DARPA/NIST broadcast transcription challenges). Search for "speech/music segmentation" or "audio segmentation" and you'll find thousands of research papers.

There are two broad approaches to solve this problem:

Supervised classification

Train a speech/music classifier, using a standard machine learning approach. You can use MFCCs as input features, along with other basic feature like zero-crossing rate, amplitude modulation at 4Hz, etc. Recently it became common to throw in as many features as possible, and using feature detection techniques to identify the most discriminant ones.

Any classification algorithm will do - support vector machines, gaussian mixture models, decision trees. Once the classification is done, you'll have misclassified frames (for example a tiny acapella segment in a song will be classified as speech; or a FX or jingle between speech will stand-out). This requires post-processing, the most common approach is to apply mode filtering (voting) on the sequence of classifier outputs. The classification/temporal smoothing are sometimes rolled into one through the use of hidden markov models for both classification and temporal smoothing.

Ref: Content-based audio classiﬁcation and segmentation by using support vector machines, Lu et al.

Unsupervised segment change detection

Consider a 10s window sliding over the signal. Compute audio features on the first half, on the second half, and use a statistical test to decide which hypothesis is the most likely: the two sets of audio features are drawn from the same distribution, or are drawn from two different distributions. The output of the test will tell you how likely it is that the middle of the window corresponds to the boundary between a speech and a music segment. Select the points with the highest scores as the segment boundaries.

The same audio features as for the supervised approach (MFCC, ZCR, amplitude modulation at 4 Hz...) can be used.

"Textbook" criterion for the statistical test: bayesian information criterion (BIC).

Ref: Unsupervised Audio Stream Segmentation And Clustering Via The Bayesian Information Criterion, Zhou & Hansen (for an introduction to BIC).

Combined Supervised and Unsupervised approaches for automatic segmentation of radiophonic audio streams, Richard, Ramona & Essid (for more exotic change detection tests).


Where to get a proper dataset collection for bikes and cars

Bike and Cars database

http://www.emt.tugraz.at/~pinz/data/GRAZ_02/

http://homepages.inf.ed.ac.uk/rbf/CVonline/Imagedbase.htm

Hope this helps .. Download the bike.zip and cars.zip files and extract it.. If it is correct, click the tick mark to accept the answer.


does audio feedback consist of harmonic content?

The contradicting statements may be a result of looking at two different aspects of the situation.  As Hilmar pointed out in an earlier answer, any clipping (be it electronic or mechanical) will indeed produce harmonics.

However, feedback distortion (not necessarily reaching to the point of the loud feedback oscillation being produced) definitely does not consist of harmonics, or even any spurious frequencies added!  Feedback distortion comes from a frequency response with a very weird structure (very unnatural to the ear, often described as adding a metallic quality to the original sound).

The reason is the following: since the sound takes a fixed amount of time to travel from the speaker to the microphones (simplified analysis assuming single speaker and single microphone, no echos/reflections, and no attenuation of sound w.r.t. distance traveled), then the frequencies at which the feedback is positive, thus producing peaks in the frequency response, are the frequencies $f$ such that $$d = k \frac{\,c\,}{f}$$ where $d$ is the distance between the speaker and the microphone, $c$ is the speed of sound, and $k$ is a positive integer.  In other words, positive feedback (and thus larger gain) occurs at frequencies where the distance between the speaker and microphone is an exact multiple of the corresponding wavelength.

Now, you see, that means positive feedback occurs at frequencies $$f_k = k \frac{\,c\,}{d}~~~~k = 1, 2, 3, \cdots$$

There, you can see that the distortion is not harmonic in nature --- first of all, the above, assuming gains are low enough that you don't produce clipping or oscillation, is linear distortion (strangely enough, seeing how the feedback distortion is a very unnatural-sounding artifact); and second, the spectral locations of the "spurious" frequencies are fixed and equi-spaced, independent of the input frequency. In fact, harmonic distortion applies to an input tone (a sinusoid at a fixed frequency), whereas feedback distortion simply introduces a "weird/unnatural" frequency response, so it applies to signals with a spectrum covering ranges of frequencies.

A feedback oscillation may be subject to a combination of both effects: it may still carry a perfectly audible "feedback distortion" quality (especially if the gain is just enough to produce the oscillation), while it may or may not produce clipping that adds harmonic distortion to the mix.


I don't quite know what you mean by different viewpoints, however, there are a lot of different similarity metrics you could use to compare images. You could try:

extracting features and comparing the features Like Scale-invariant feature transform (SIFT) features
you could try dimensionality reduction and comparing that representation, like taking the discrete cosine transform (DCT).
I know you said you tried mean squared error but you could also try comparing the structural similarity (SSIM).

Edit:

Now that I better understand your question, I believe I found a resource that could be of assistance. Here is a very aptly named paper that is titled "A Novel Algorithm for View and Illumination
Invariant Image Matching".


In z-transform, if z means "delay", why do we talk about the value of z?

It is a bit to wrap your head around.

Mathematically, if you denote the z transform as $X(z) = \mathcal{Z}\left \lbrace x_k\right \rbrace$, then when you take the transform of $x_k$ after it's been delayed by one sample you get $\mathcal{Z}\left \lbrace x_{k-1}\right \rbrace = z^{-1}X(z)$.  Except in strange corner cases* this always works, and you don't need to know $X(z)$ to know it works.

So $z^{-1}$ (and thus $z$) are frequency-domain variables, just like the $k$ that gets lost in the transform from $x_k$ to $X(z)$ is a time-domain variable.  But $z^{-1}$ acts exactly like a delay operator, so in most circumstances you can treat it as such.

* Which I can't even call out -- you simply don't run into them in practice.


How can I detect buttons in a GUI of an app with pattern detection using python?

First, have a look at the squares.py sample provided by OpenCV. It should handle a fair number of button types with some tweaking.

Here is the output I got (with some tweaking) for your Calculator example:

I made the following tweaks to the squares application:

Change this code (starting on line 84):

if(result.total == 4 and
   abs(cv.ContourArea(result)) > 1000 and
   cv.CheckContourConvexity(result)):

To this:

if(result.total == 4 and
   abs(cv.ContourArea(result)) < 1300 and
   abs(cv.ContourArea(result)) > 300 and
   cv.CheckContourConvexity(result)):

Also, since you don't have to worry about scale or rotation variations, check out using matchTemplate. Also, its corresponding tutorial (in C++) is here.

See my other answer for another example of how template matching works. Also, you might find this answer helpful for detecting the top X matches with matchTemplate.

Hope that helps!


Best (Perceptually / Objectively) Super Resolution Methods Out There?

In our days the Deep Neural Network methods certainly are generating best results.
Due to the intense research going on in this field the best method is a moving target hence one can not pin point to one.
One generation before them the best methods were based on Dictionary Learning.
For example you can use the K-SVD for Single Image Super Resolution.
Those approaches / methods works pretty well and are even improved using the EPLL approach.
You should be aware those methods are really slow for real world images.
Probably you should invest in the DNN world.
Pick any one of the latest on ArXiv Sanity.


Upper Bound for the DFT (FFT) Coefficients of a Bounded Signal

You can obtain a bound on the magnitude of the DFT of $x[n]$ ($|x[n]|\le 1$) as follows:

$$\big|X[k]\big|=\left|\sum_{n=0}^{N-1}x[n]e^{-jkn2\pi/N}\right|\le\sum_{n=0}^{N-1}|x[n]|\le N\max_n|x[n]|=N\tag{1}$$

Note that for signals $x[n]=e^{j2\pi nl/N}$, $l\in\mathbb{Z}$, the bound is tight. However, for most other signals with $|x[n]|\le 1$ you will find that the maximum of $|X[k]|$ is substantially smaller than the bound given by $(1)$.