Time recovery algorithm and a symbol with samples

Preamble: This answer is about timing recovery in a sense of symbol synchronization, i.e. finding the proper sampling phase of a baseband signal.
Based on the stated requirement of only 8 samples per symbol, I will assume that you are employing a fully digital approach to timing recovery. That means that you have no control over the times of sampling by ADC. The fully digital approach looks like this (this nice picture is from slides here):

The "digital processor" is your FPGA-based implementation of the timing recovery algorithm, including the selected Gardner algorithm to find a phase error. "Sampler" - is an ADC with fixed sampling rate, 8 samples per symbol. "Analog processor" - this is an analog front-end, which is implemented in your board outside of the FPGA. And "signal in" is a baseband signal from your target communication channel.
In order to verify an implementation of your algorithm, you will want to model it, e.g. using a Matlab program or a Python script (whatever is more convenient). With the real signal samples from ADC you will be able to see the details of how your algorithm works in a long run, which otherwise is difficult to simulate in RTL and may be hard to debug when programmed on board. With a synthetic signal you can vary the parameters of the signal to see, how the algorithm reacts to it, and also you can plot a sigmoid for the TED, which can be very helpful.
About TED:
Timing error detector (TED) is intended to detect a phase error in sampling. Based on the sign of its output, the timing recovery scheme will incrementally shift the sampling phase by a little step value, until the timing error becomes zero (or close to zero). This step value is usually selected to be 1-2% of a symbol period. Note, that with a fully digital timing recovery approach, when there is no actual control of the sampling phase, an interpolator is used instead, to fill the gaps between the signal samples. As a result the phase shift is made for the recreated version signal and not the real one.
Sampling phase error is defined by a phase difference between the transmitter, which generates symbols, and base sampling frequency of ADC (actual sampling frequency/8). Ideally the algorithm should be able to converge for arbitrary initial value of that error.
The answer to the question:
In order to adequately model the initial phase error, the model must include not only a digital processor, but also a sampler. This means, that the input signal for the model should be continuous or oversampled. The degree of oversampling will define the resolution ability of your tests, i.e. how many different initial values you can provide. The goal can be 50-100 points per symbol period.
If you intended to use a real signal as an input of your model, you would use a logic analyzer to capture the signal at the input of the ADC with high enough frequency.
Alternatively, you may want to use a synthetic signal, that models a real one. Then you have to create its mathematical model. Usually you will start with random symbol values (as you did), then apply the modulation scheme (you will skip this if the symbol is a signal level), then apply a pulse shape (if the target transmitter does), then apply the distortions/pulse shaping relevant to your intended target communicating channel and analog front end (e.g. a low-pass filter can model the distortion of the signal in a twisted pair).
As you can see, the particular way of creating a synthetic signal will strongly depend on the type of the target communication channel.
Edit: Here is an illustration of what the resulting input data $x(k)$ for the TED may look like. It shows two different cases. You can see that they differ in the phase of ADC sampling. Also the examples show that the selection of 2 samples out of 8, which is necessary for the Gardner TED, is random.

The black lines are the exemplary eye diagram of the signal, the green line is the optimal sampling point and the blue lines are the ADC samples.


FFT on non-periodic signal and signal power at each harmonic - is there a more robust approach?

I'll try to address your questions one at a time.

  Is it suitable to calculate the fundamental frequency from the length of the signal?

Your signal is not periodic, so speaking of "fundamental frequency" is not correct. The fundamental frequency corresponds to the lowest frequency of a periodic waveform, which is not the case here. So, regarding your question, no: you cannot define a fundamental frequency from the length of a signal - two things not related at all.

  I have read that it is not suitable to use an FFT on a non-periodic signal. Is there a way to adapt the FFT method to make it more suitable e.g. some sort of windowing function?

FFT can indeed be used for non-periodic signals. In case you are not too familiar with DFT, I'll try to explain how it works in a nutshell.

You have a discrete signal. You window it (you only take $N$ of its values). You perform the DFT on those values. What does "perform the DFT" mean? The signal whose DFT will be calculated is assumed to be $N$-periodic, and you get a discrete frequency representation of that signal (not the original non-windowed one). I believe that this is the reason why you state that FFT is not suitable on non-periodic signals. Nevertheless, you can still use it, we do it all the time. Maybe some answers in this question can help.

If you are interested in which frequencies are present at each instant of time, one can use several windows an overlap them, and take the FFTs of those signals. See STFT.

  Is there an easy way to remove drift/bias from the signal?

Do you mean DC value? In that case, you can just substract the mean value of the signal, or apply a high-pass filter with sufficiently low cutoff frequency.

  I would also welcome comments on any other aspects of this process if you feel it would make it more robust.

At some point in your script, you just set some FFT coefficients to zero. Note that this approach is generally not a good idea.


Meaning of tight frame in the context of texture synthesis

Since visoft requested i post the solution, here goes. It has been a while so not all the details are still fresh so feel free to ask anything if things are unclear!

Let's say we have an image $I$ and apply a filter, which is given by convolution by some $F$. The pixel values of $I$ convolved with $F$ are all the output of the inner product of $I$ with a translated (and less importantly: mirrored) version of $F$.

What i missed at the time was that each pixel should be seen as a coefficient of the new representation of the image. So instead of the 16 subbands i mentioned, the number of coefficients of our new representation is 16 times the number of pixels in the image.

Some notation: for pixel number $(a,b)$, let us write $F_{(a,b)}$ for $F$ translated $(a,b)$ steps (and mirrored) and $(I * F)(a,b)$ for the $(a,b)$-th pixel of $I$ convolved with $F$.

Also, let us call all the translated versions of $F$ our basisfunctions. One can check from the definition of convolution that if we convolve our filter output again with $F$, this results in multiplying every basisfunction $F_{(a,b)}$ by the coefficient $(I * F)(a,b)$ and adding them up. This results in
\begin{align*}
I * F * F &= \sum_{(a,b)} (I * F)(a,b) \cdot F_{(a,b)}\\
&= \sum_{(a,b)} \langle I, F_{(a,b)}\rangle \cdot F_{(a,b)}
\end{align*}
where the last equality follows by definition of convolution.

From this, one easily recognizes something that comes up using an orthogonal basis. Now the collection of $F_{(a,b)}$ for all $(a,b)$ is called a tight frame if the above mathematical expression is equal to $I$ itself! It is thus very similar to an orthogonal basis, except that one does not need the set $\{F_{(a,b)}\}_{(a,b)}$ to be linearly independent (if you are very sharp you will note that in this particular case it most likely will be linearly independent, but see the following).

In the paper, actually about 16 of these $F's$ are used. Let's number them $F_1$ up to $F_n$. This is called a tight frame if the following holds:
$$
I = \sum_n \sum_{(a,b)} \langle I, F_{n,(a,b)} \rangle \cdot F_{n,(a,b)}
$$

This expression is equal to (check this):
$$
I = \sum_n (I * F_n * F_n).
$$

Since convolution is equivalent to multiplication in the fourier domain, one can check this by going to the fourier domain, squaring all the $F_n$ pointwise and adding them up. That's simple to do. If that adds to an image that is equal to 1 in every pixel, we are dealing with a tight frame.

So concluding, i missed two things: convolving with one $F$ gives as much coefficients as there are pixels in the image, and convolving twice actually multiplies the basisfunctions with the coefficients and adds them up.

As I said, I wrote this quickly and it has been a while, so if anything is unclear feel free to ask again!


Here are three methods of finding duplicates in a database of audio files. If I've understood your question correctly, I think you want the cross-correlation approach.

Looking for identical files - md5 checksum (not a DSP solution I know)

Looking for identical audio samples but time-shifted - cross correlation

Looking for similar sounding audio but different samples - audio fingerprinting.


Is it possible to define Fourier or wavelet transforms on DNA sequences?

First, a code for paper A measure of DNA sequence similarity by Fourier Transform with applications on hierarchical clustering. Journal of theoretical biology is available at GenomeDFT (github).
Then, on the transforms. Discrete Fourier transforms can be defined on other domains than vector spaces. There are Fourier transforms on groups, rings or finite fields. On the latter, it is commonly called a number-theoretic transform (NTT). This possibility entails to making sense of group/ring/field axioms for the nucleotides. For instance, it should be possible to use a finite field of characteristic four. Apparently weird mathematical structures could be interesting when considering (C,G) or (A,T)  patterns, or the redundancy in the codon/amino-acid correspondence.
Some papers mention the Ramanujan-Fourier transform (A Novel Method for Comparative Analysis of DNA Sequences by Ramanujan-Fourier Transform, apparently with the same first author of the paper you mentioned). With less structured transformations, it is also possible to evaluate periodicity, as for instance in Categorical spectral analysis of periodicity in nucleosomal DNA.
There also were some works on wavelets, either continuous or discrete. On the first flavor, by A. Arneodo and colleagues, and more recently in 2019 with
Wavelet-Based Genomic Signal Processing for Centromere Identification and Hypothesis Generation. With discrete wavelet transform, there is for instance Identification of exonic regions in DNA sequences using cross-correlation and noise suppression by discrete wavelet transform, Application of discrete wavelet transform for analysis of genomic sequences of Mycobacterium tuberculosis, Wavelet analysis of DNA sequences.
Finally, on the mappings. There is for instance Digital signal processing methods for biosequence comparison, 1991 or Mapping Equivalence for Symbolic Sequences: Theory and Applications, 2009.
As asked by Marcus, your actual purpose might direct you to an appropriate combination.


Can a Kalman Filter Be Applied Using Measurement-Space Dependent Sensors?

Yes, sorta kinda.  Basically, if you can model your system accurately and the resulting model fits the Kalman filter paradigm, then yes.

However, I think it's a mistake to give the filter position and velocity -- were I doing this, I'd just give it the wheel position and let the Kalman filter determine the velocity -- which it will do well to the extent that your model matches the actual plant.

If it'll work at all, it'll be because you take the typical Kalman system model, $\mathbf{x}_k = \mathbf{F}_k \mathbf{x}_{k-1} + \mathbf{B u}_k + \mathbf{w}_k$, and you (probably, unless you have some prior idea of the acceleration being applied) let $\mathbf{B} = \mathbf{0}$, let $\mathbf{w}_k$ reflect your maximum expected acceleration, and let $k$ be the triggering instance.

The first fly in this ointment is that you have to calculate $\mathbf{F}_k$ and $\mathbf{w}_k$ according to the amount of time elapsed since the last trigger instant.  So for $\mathbf{x} = \left [ \begin{matrix} \theta & \omega \end{matrix} \right ]^T$, where $\theta$ is angular position and $\omega$ is angular velocity, you'd get something like $$\mathbf{F}_k = \left [
\begin{matrix} 1 & t \\ 0 & 1
\end{matrix}\right ]$$ and $$\mathbf{w}_k = \left [
\begin{matrix} t^2 & t \\ t & 1
\end{matrix}\right ] \sigma_\alpha^2$$ where $\sigma_\alpha$ is your expected variation in acceleration.

(Note that I'm not giving the derivation of $\mathbf{w}_k$, for which I only sorta-kinda-apologize -- it's tedious for me to translate my hen-scratches to properly formatted math, but it'll be good for you to derive it yourself.  Start with the definition of the Kalman-Bucy filter and grind through the math to find $\mathbf{P}(t)$ for a constant $\mathbf{w}(t)$ -- the above expression for the effective sampled-time $\mathbf{w}_k$ should drop out).

Note that the second fly in the ointment is that as the speed of the wheel drops the effect of the acceleration uncertainty ($\sigma_\alpha$, above) will get bigger and bigger compared to the position change from sample to sample.  At the point where the likelihood of overrunning the sample gets significant, you can't really use a "pure" Kalman filter anymore -- logically, this is because if you haven't triggered a pulse, then you know darned well that your position and velocity excursion are limited, but there's no way to effectively communicate that to a pure Kalman application.

If you need to track a really slow wheel, then you need to come up with some other, nonlinear, state estimator.

PART II

On reflection, I thought of a possibly-better way.  It's still not the best way, because you're jamming a nonlinear system with non-Gaussian uncertainty into a model of a linear system with Gaussian noise -- but it has the advantage both of being conceptually simpler, and of working down to zero speed.

The system model that you use to build a Kalman filter doesn't have to be time-invariant.  So model the system as $$\begin{align}
\mathbf{x}_k = \mathbf{F} \mathbf{x}_{k-1} + \mathbf{B u}_k + \mathbf{w} \\
\mathbf{y}_k = \mathbf{H} \mathbf{x}_k + \mathbf{v}_k
\end{align}$$

Note that $\mathbf{v}_k$ is modeled as time-varying.  In truth, the system is nonlinear (you get a pulse at specific values of $\theta$), but we'll just pretend.

Run the Kalman at a constant rate in time, with $\mathbf{F}$, $\mathbf{B}$, and $\mathbf{w}$ $\mathbf{H} = \left [ \begin{matrix} 1 & 0 \end{matrix} \right ]$ all constant.

But for each sample where a pulse is triggered, set $\mathbf{y}$ equal to the known position of the wheel, and $\mathbf{v}_k$ small, to reflect that you darned well know where the wheel is at that moment (size $\mathbf{v}$ to reflect any uncertainties in position due to mounting or irregular triggering or whatever).

For each sample where a pulse is not triggered, set $\mathbf{v}_k$ large, to reflect the fact that you don't know where between the triggers the wheel may be.  I'd start with $\mathbf{v}_k = \left( 3 \theta_0 \right)^2$, with $\theta_0$ equal to the spacing between trigger points -- and note that the multiply by three is coming directly out of my ear, so feel free to experiment with it.  Set $\mathbf{y}$ equal exactly halfway between the last known position of the wheel and the expected next position of the wheel (I'm assuming you know the wheel direction).

The net effect of this should be that if the wheel is turning fairly quickly the filter will have a darned good idea of the position and velocity, with "coasting" between readings.  If the wheel is turning slowly, it should coast OK if it's going at a constant rate, it should update at least fairly well if it speeds up, and if it slows down or stops, the fact that you're telling the Kalman filter that it's between positions "somewhere" should cause the Kalman's velocity estimate to slowly decay to zero (which is your best guess, intuitively), while the position estimate slowly settles on halfway between triggers (which again, is intuitively your best guess).


Blind source separation with one source rotated

Consider 2 points which are rotated version of each other, for these 2 point you could write 4 equation with 4 unknown and solve this system of linear equations (b=A*x, where b obtained from I1 and I2 at those 2 symmetric points and x are values of A and B at those points and the matrix A is obtained from Wij). Now you could solve this system of 4 linear equations for each pair of pixels, but by doing this you can't take into account the noise level unless you repeat your measurements!

To take into account the noise levels, you could use some measure of noise level like Total variation. You could build a cost function by adding the norm of (b-A*x) to a measure of noise level of your images then trying to find the unknowns by minimizing the cost function. In this case you cant use only a pair of points and you have to write the system of equations for all of points.


About the proof of an equality related to the DFT [sampling the DTFT to obtain the DFT]

The right hand side is the periodic repetition of the left, tells the last equation.

Starting from Fourier transform $$X(\omega)=\sum_{n=-\infty}^{\infty}x(n)e^{-j\omega n}$$ and $X(\omega)$ is periodic with period $2\pi$. So only the samples in the fundamental frequency range are sufficient.

Now we take $N$ equi-distance samples in the iterval $0\leq\omega\leq 2\pi$, so the spacing will be $\dfrac{2\pi}{N}$

Let $k$ be the sample index and $k=0,1,\dots,N-1$ and we evaluate $\omega$ at $\left(\dfrac{2\pi k}{N}\right)$

$$X\left(\dfrac{2\pi}{N}k\right)=\sum_{n=-\infty}^{\infty}x(n)e^{-j2\pi kn/N}$$

which can be written as

$$X\left(\dfrac{2\pi}{N}k\right)=\sum_{l=-\infty}^{\infty}\sum_{n=lN}^{lN+N-1}x(n)e^{-j2\pi kn/N}$$

Now change the index $n$ to $n-lN$ which results in

$$X\left(\dfrac{2\pi}{N}k\right)=\sum_{n=0}^{N-1}\left[\sum_{l=-\infty}^{\infty}x(n-lN)\right]e^{-j2\pi kn/N}$$

The sequence inside the square brackets is the periodic repetition of $x(n)$ of period $N$

$$x_p(n)=\sum_{l=-\infty}^{\infty}x(n-Nl)$$

Thus $$x_p(n)=\begin{cases}
               x(n),\;\;\;\;0\leq n \leq N-1\\
                     \\
               0,\;\;\;\;\;\;\;\;\;\;\;     \mathrm{elsewhere}\\
            \end{cases}$$


Sampling Theorem: How to know the value between two samples of a Signal

Interpolation using a Sinc function kernel.  This assumes that the sampled signal was perfectly bandlimited to below half the sampling rate.

Note that perfectly bandlimited signal are infinite in extent.  For finite-length "real world" signals, using a windowed Sinc interpolation kernel (thus a finite computation with a noise floor) is a common method to approximate points between samples (for upsampling, for instance).


How a diode-detector affects to Frequency Response?

But what does detector in the frequency domain?

Something really complicated.  We get trained to think about signal processing almost exclusively in the frequency domain, and we forget that the frequency domain stuff was only invented as a mathematical trick to make it easier to understand what's going on in the time domain.  If it's easier to do in the time domain, do it there.
In fact, envelope detection is really horrid from a frequency domain point of view, because not only are there nonlinearities in play, but there's nonlinearities and that perfectly linear capacitor that are interacting.  So generally if you want to analyze AM radio demodulation you try to find some equivalent operation that is linear, or at least less horridly nonlinear.
In this case, just considering the demodulation to be clipping isn't bad.  Figure that the clipping operation can be denoted as $y(t) = f\left( x(t) \right )$, where $f(x)$ is (thankfully) memoryless.  So you can assert* that it has a Taylor's expansion $$f(x) = \sum_{n=0}^\infty h_n x^n$$ where the $h_n$ are the terms of the expansion**
If we feed a sum of sinusoids through $f(x)$, the above expansion tells us that we're going to get just about every possible distortion product out of the thing.  We'll get stuff at the sum of all the sinusoid frequencies, at every possible combination of sum and difference frequencies, and at ever possible combination of the sum and difference of every possible harmonic of both sinusoids.
At this point, I'm just going to claim without proof that the element at the difference between the carrier and the sideband frequencies will be prevalent -- if the carrier and the modulating signal are widely spaced apart.  AM radio works because you're transmitting at 500kHz or more, and modulating the carrier with a signal that only occupies 0Hz to 3000Hz or so.
I'm not going to try to prove this in more detail because the math gets crazy in the frequency domain -- and because you can so clearly understand the operation in the time domain.

Why the example in audio editor doesn't work?

First, because you just multiplied $\cos \omega_c t$ with $\cos \omega_s t$  For AM modulation, you want something like $\left(\cos \omega_c t\right ) \frac{s(t) + 1}{2}$ -- you multiply the carrier by something that never falls below zero.
Second, because your carrier and your modulating tone are so closely spaced (although the signal may be recognizable).
Here's some plots.

Top, a picture of one 1kHz cycle's worth of your signal, as you've defined it.  Because you didn't shift your modulating tone, you don't have enough carrier to form an envelope.

My version with your frequencies, your modulating tone shifted by 1 and reduced by half.

My version with a carrier frequency of 20kHz, and the modulating tone shifted.

* You can assert this, but you'd be lying, because it has derivatives that go to infinity -- but you can make up something arbitrarily close to a clipping function that does have a Taylor's expansion -- just assume I've done that.
** And I use $h$ in honor of the horrid things I'm doing to the math...


How can I get the uncertainties for peaks on an image?

One way is to simply model each peak with a Gaussian, with mean $\mu_i$ and variance $\sigma_i$. In fact what you mean by uncertainty corresponds to the variance. You can iteratively fit Gaussians using e.g. EM-algorithm. In MATLAB you could easily do this with built-in fitting functions: https://www.mathworks.com/help/curvefit/gaussian.html

Alternatively, you could use the mean-shift algorithm to find the modes in the data, from which you could extract the probabilities:
http://de.mathworks.com/matlabcentral/fileexchange/39079-mean-shift-for-finding-modes

And finally, we might want to think of the uncertainty as the steepness of the maximum. In that regard, the magnitude of the first derivative would give us how strong the peak is. Then, a first order derivative, or a curvature filter can characterize the uncertainty. This is also doable in MATLAB with simple finite difference approximation schemes.


Applying a filter on an audio signal with Python

Apart from what Marcus suggests in the comments, ft[i] *= H(freq[i]) is probably wrong.  The overlap-add and overlap-save methods should deal with this correctly, but read on!

The FFT does not implement linear convolution. It implements circular convolution, which can approximate linear convolution provided certain conditions are met. One of those conditions is that $$N_{\tt fft} \ge N_{\tt filter} + N_{\tt signal} - 1.$$

In your case, you've effectively set $$N_{\tt fft} = N_{\tt filter} = N_{\tt signal}$$ which will introduce time-aliasing --- the problem that happens when the "certain conditions" are not met.  See here for an answer to a similar question.


Biquad Treble Shelf filter gain being weird

It's your scaling. Don't divide the coefficients by 2.336268076  (which is exactly the gain that's missing.

Instead use proper fixed point implementation techniques (scaling, shifting, "Q" notation, etc.) to deal with numbers larger than 1. Or much easier: use floating point, if you can.


Detecting moved and dirty areas between two snapshots/images

I love abstract art.  It speaks to something inside of me that is outside the library, it speaks in words but not words.  One of my favorite artists is Wassily Kandinsky.

An image of his that could meet the requirements of the problem, both with "text" and translation of a circle, is his "Circle within a circle".  The bouncing globe animation is pristine, not textured, and doesn't have the text.

Also, the real world is messy and complex, so outside of explicitly specifying "this must work only on computer-generated content", any viable method has to hand a few real-world curve-balls.

Here is a copy of it:

So using MS Paint, I moved one of the circles and hashed up the lines in a way that should be consistent with local text changes.

Here it is:

So lets open up the files, convert to grayscale, and subtract them.

Here is some code:

import cv2
import numpy as np
from matplotlib import pyplot as plt
from PIL import Image

#read in images
#name of input file
fname1 = './wassily-kandinsky(A).png'
fname2 = './wassily-kandinsky(B).png'

#read the files, convert to grayscale
img1 = Image.open(fname1).convert("L")
img2 = Image.open(fname2).convert("L")

#convert to numpy array for math
arr1 = np.asarray(img1)
arr2 = np.asarray(img2)

#smooth (to remove some of the manual-editing real-world images artifacts)
kernel_width = 11
arr1b = cv2.bilateralFilter(arr1, kernel_width,75,75)
arr2b = cv2.bilateralFilter(arr2, kernel_width,75,75)

#subtract images
arr_diff = np.abs(arr2b-arr1b)

#show them
plt.imshow(arr1, cmap='gray')
plt.savefig("figure1.png", dpi=300)
plt.show()

plt.imshow(arr2, cmap='gray')
plt.savefig("figure2.png", dpi=300)
plt.show()

plt.imshow(arr_diff, cmap='gray')
plt.savefig("figure_diff.png", dpi=300)
plt.show()

And this gives slightly different initial images, and allows a decent difference image.

New (smoothed, grayscale) "unchanged" image:

New (smoothed, grayscale) "changed" image:

The difference between the two:

Here is a zoom in of the "signal area":

Next steps (to do when I get time):

clean up circle signal, show how edge histograms clearly give locations
block location from difference image, show locations of "salt-and-pepper" which would come from alternate text in same location.


Solving Convex Optimization Problem Used for High Quality Denoising

Boyd has A Matlab Solver for Large-Scale ℓ1-Regularized Least Squares Problems. The problem formulation in there is slightly different, but the method can be applied for the problem.

Classical majorization-minimization approach also works well. This corresponds to iteratively perform soft-thresholding (for TV, clipping).

The solutions can be seen from the links. However, there are many methods to minimize these functionals by the extensive use of optimisation literature.

PS: As mentioned in other comments, FISTA will work well. Another 'very fast' algorithm family is primal-dual algorithms. You can see the interesting paper of Chambolle for an example, however there are plethora of research papers on primal-dual methods for linear inverse problem formulations.


Does whitening with PCA correspond to power spectrum flattening?

I believe the answer here is no.

Remember what 'whitening' is really doing: It is ensuring that the power/variance that correspond to each of your basis vectors are the same. Which basis vectors? The ones you chose.

In the DFT, you are choosing basis vectors comprised of complex exponential functions. In PCA, you are determining the basis vectors empirically. What is the relationship between complex exponential bases and empirical bases? You cannot say.

Thus, when you say you whitened in the spectral domain, you are really saying "I made sure the power/variance associated with all my complex exponential basis vectors representing my data are the same".

When you say you whitened in the PCA domain, you are really saying "I made sure the power/variance associated with of all my empirically derived basis vectors are the same".

Are your empirically derived basis vectors complex exponentials? There is no guarantee. Thus, whitening in either domain is not equivalent to whitening in the other.

  If so, are the eigenvalues of the data covariance (which I understand
  are also referred to as a spectrum) directly related to the mean power
  spectrum of my signals ? How ?

If I understand your question correctly, then the answer to this is a simple yes, for the simple reason of energy conservation. All the eigenvalues together will give you the total energy/variance in your data. The average of this energy, will have to correspond to your mean power in any domain, including the spectral.

Projecting your data onto either the DFT or PCA (empirical) bases dones not add/remove energy, but simply redistributes it. Thus, your average power has to always be the same.


Is there a way to separate 2 signals recorded on the same channel?

You should look into a generalization of ICA called Independent Subspace Analysis (ISA). Practically it is Principal Component Analysis (PCA) followed by ICA.

There are threads on dsp.stackexchange that discuss this as well as github repos (and here) that provide various implementations.


Confused in Bit rate, Symbol rate, Sample rate and Data rate

The $BT$ product is linked to the single-side $-3\rm dB$ down bandwidth and the bit rate $R_b$ with the relation below:
$$BT = \frac{f_{-\rm 3 dB}}{R_b}$$
Assuming you have a single-side $-3\rm dB$ bandwidth, you can extract your bit rate as:
$$R_b = \frac{f_{-\rm 3 dB}}{BT}$$
Also see the answer to this question. For more on GMSK see this paper $^1$ and this paper $^2$.

You have one bit per symbol, meaning for every bit you have one symbol. And their rate should be in that proportion as well. Or see it this way,
\begin{align}
\text{For} \ \ &n \ \ \text{bits you have} \ \ n \ \ \text{symbols}\\
\implies \text{for} \ \ &\Delta n \ \ \text{bits you have} \ \ \Delta n \ \ \text{symbols}\\
\implies \text{for} \ \ &\frac{\Delta n \ \rm bits }{\Delta t} \ \ \text{you have} \ \ \frac{\Delta n \ \rm symbols}{\Delta t}
\end{align}

See the sampling rate as the rate at which you are having your samples, and this value is in $\rm Hz$, or basically
$$
\frac{\Delta n \ \rm Samples}{\Delta t}
$$
The $\Delta t$s above are normally given as $1 \ \rm second$. The sampling rate is the number of samples per second, also known as sampling frequency (in $\rm Hz$).

Generally speaking, your bit rate is your data rate. The answer to this question details more on this.

You can use the maths and plug in your numbers.

$[1]$: A. Linz and A. Hendrickson, "Efficient implementation of an I-Q GMSK modulator," in IEEE Transactions on Circuits and Systems II: Analog and Digital Signal Processing, vol. 43, no. 1, pp. 14-23, Jan 1996.

$[2]$: K. Murota and K. Hirade, "GMSK Modulation for Digital Mobile Radio Telephony," in IEEE Transactions on Communications, vol. 29, no. 7, pp. 1044-1050, Jul 1981.


What features describe audio signals? (Besides frequency and amplitude)

Two remarks:

I am assuming you are plotting the real (or imaginary) part of the Fourier transform. It is much more common to work with the magnitude or squared magnitude (power spectrum).
The peak in the spectrum is a very poor measure of fundamental frequency (pitch). Take a piano note at 440 Hz, apply a notch filter to it to remove the 440 Hz component. Even if the peak in the spectrum is now at 880 Hz (roughly), it will still sound like a 440 Hz note - yes, at a frequency totally absent from the spectrum!

To answer your question, timbre is the property that distinguish sounds of same loudness and pitch. However, it is not an unidimensional characteristic with a well-defined unit. What you can extract are features such that things that "sound similar" have similar features.

Common features used for characterizing timbre include:

Spectral centroid, which indicate how "dark" or "bright" a sound will be perceived, spectral spread (a measure of bandwidth, tonality vs noisiness), and moments of order 3 and 4 (kurtosis and skewness).
Energy in each Bark band, or ratio of energies between consecutive Bark bands.
Any low-dimensionality descriptors of spectral envelope, including autoregressive coefficients, or Mel-Frequency cepstral coefficients.

See section 6 of this document.


What is the difference between a range image and a depth map?

Your understanding is correct, it's z-values vs distances.

A depth map contains the distances of points in the scene to a plane.

A range image contains the distances of points in the scene to another point. This point is very often the camera's position.

For a depth map the plane that the distances are computed with respect to can be any plane, but it is often parallel to the image plane. Examples of common planes:

The image plane
The focal plane
A plane closer to the object, such that the points closest to the camera have distance zero.

If you of the same scene have:

A range image with the distances to the camera
A depth map with  the distance to the image plane

these two images will look very similar, and it's not possible to conclude whether it's a range image or depth map, unless you have something you know is a plane in the image.


The FFT is faster than the naive DFT through matrix-vector products because it can reuse many intermediate results.
However, for inputs this sparse, there's really not much that can be re-used, even in the best case.
So, the most efficient way here is probably to work straight forward:

Allocate the space ($N=7408800$ values) for your output (initialize with zeros) $o$.
Go through your input $i$. Nested loops:

For every nonzero input index $n$:

For every output index $k$: add $i_ne^{j2\pi \frac{nk}N}$ to $o_k$.

Notes:

since with sparse input, your output is non-sparse, it's usually faster to traverse the input only once, but the output many times, instead of the other way around (that is because memory is much faster when addressed linearly)
to calculate the complex exponential, you can use SIMD-accelerated calculations, because the values are in consecutive addresses and have proportional phase. LibVOLK has a number of implementations that might be helpful for fast calculation here.


Here pretty good explanation.

I'll start by consider the 1D case of cubic interpolation (because its easier to explain) and then go onto the 2D case.

The basic idea of cubic interpolation is to estimate the values between any two points by as a cubic function, $f(x)$, with first derivative, $f'(x)$

$$f(x) = ax^3 + bx^2 + cx + d$$
$$f'(x) = 3ax^2 + 2bx + c$$

so that, when spliced to the adjacent cubic splines on the left and on the right, the interpolated function (the "zeroth derivative") and its first derivative are continuous.  That means, at the left boundary of some given spline, the zeroth and first derivatives are equal to the those of the right boundary of the adjacent spline immediately to the left.

Without loss of generality, we can assume that we are interpolating between $x=0$ and $x=1$.  On the left and right boundaries, we see that

$$f(0) = d$$
$$f(1) = a+b+c+d$$
$$f'(0) = c$$
$$f'(1) = 3a+2b+c$$

Solving for the polynomial coefficients

$$a = 2f(0) - 2f(1) + f'(0) + f'(1)$$
$$b = -3f(0) + 3f(1) - 2f'(0) - f'(1)$$
$$c = f'(0)$$
$$d = f(0)$$

Now consider the four points around the region being estimated, which are $(-1,y_{-1})$, $(0,y_0)$, $(1,y_1)$, $(2,y_2)$, and impose boundary conditions that

$$f(0) = y_0$$
$$f(1) = y_1$$

and define the two boundary derivatives as

$$f'(0) = \frac{y_1-y_{-1}}{2}$$
$$f'(1) = \frac{y_2-y_0}{2}$$

which is consistent to how they will be defined for the adjacent splines.

Then

$$a = -\frac{y_{-1}}{2} + \frac{3y_0}{2} - \frac{3y_1}{2} + \frac{y_2}{2}$$
$$b = y_{-1} - \frac{5y_0}{2} + 2y_1 - \frac{y_2}{2}$$
$$c = -\frac{y_{-1}}{2} + \frac{y_1}{2}$$
$$d = y_0$$

once the coefficients are known you can interpolate $f(x)$ for any point between the central 2 points, $(0,y_0)$ and $(1,y_1)$.

For bicubic interpolation the principle is pretty much the same but you estimate a surface using 16 points (4x4 grid) rather than just a curve. One simple way to do this is to first interpolate the columns and then interpolate the resulting rows.


Image thresholding while retaining "underneath" contents

This is by no means a complete answer but I think it will get you far enough to achieve what you want. Edit: maybe not.

find main colours using kmeans
assign main colours to different types (brown line, blue line, etc). Some appear in both
clean up with morphology
segment image

What I thought might work was to pick out the most common colours and classify each pixel by the nearest colour. Then classify each pixel by type (brown line, etc), where the colour class can belong to multiple type classes. This would hopefully solve the problem of gaps in the lines - i.e. the colour where blue crosses brown would belong to both type class and fill the gap.

Not sure if it really does a better job. The issue with kmeans is that it uses random starting points so the colour classes were changing every time the script was run. I think I've fixed that now.

I wonder if there is some code out there to connect the broken lines?

% Read file
I = imread('mapdsp.jpg');

% Reshape into 2D matrix of pixel colours
ir = reshape(I, [size(I,1) * size(I,2),3]);

% Choose exact same sample points everytime
rng(0);
ir_indx = randi(size(ir,1),20,1);
ir_startpoints = ir(ir_indx,:);

% Use k-means clustering to find most common colours
[ik,c,sk,sd] = kmeans(double(ir),20,'Start',double(ir_startpoints));

% Show colours
cim = reshape(uint8(c),[1,20,3]);
image(cim); title('Colours'); pause;

% Sort colours
chsv = rgb2hsv(c);
[~,cni] = sort(chsv(:,1),1);
cn = c(cni,:);
cnim = reshape(uint8(cn),[1,20,3]);
image(cnim); title('Sorted colours'); pause;

% Get index of index
[~,cni2] = sort(cni);

% Show pixels classified
in = reshape(cni2(ik,:), [size(I,1),size(I,2)]);
imagesc(in);  title('Pixel classification'); pause;

% Show each class
In = zeros(size(I));
for k = 1:20
    In(:,:,1) = I(:,:,1) .* uint8(in == k);
    In(:,:,2) = I(:,:,2) .* uint8(in == k);
    In(:,:,3) = I(:,:,3) .* uint8(in == k);
    image(uint8(In)); title (['Class ' num2str(k)]); pause;
    im_dist(:,:,k) = reshape(sd(:,k), [size(I,1),size(I,2)]);
    imagesc(-im_dist(:,:,k) .* (in == k));  title('Pixel distance'); pause;
end

% 12,19 and 20 blue lines
im_blue = (in == 13) | (in == 19) | (in == 20);
imagesc(im_blue); title('Blue lines'); colormap gray; pause;
im_blue = imopen(im_blue,strel('diamond',1));
imagesc(im_blue); title('Blue lines opened'); colormap gray; pause;
for k=1:3; In(:,:,k) = I(:,:,k) .* uint8(im_blue); end;
image(uint8(In)); pause;
imagesc(bwmorph(im_blue,'skel',Inf));

% brown lines
im_brown = (in == 1) | (in == 2) | (in == 1) | ...
    (in == 7) | (in == 8) | (in == 11) | (in == 12) | (in == 13);
imagesc(im_brown); title('Brown lines'); pause;
im_brown = imopen(im_brown,strel('diamond',1));
imagesc(im_brown); title('Brown lines opened'); pause;
for k=1:3; In(:,:,k) = I(:,:,k) .* uint8(im_brown); end;
image(uint8(In)); pause;
imagesc(bwmorph(im_brown,'skel',Inf));

% 9 - black line
im_black = (in == 10);
imagesc(im_black); title('Black lines'); pause;
im_black = imopen(im_black,strel('diamond',1));
imagesc(im_black); title('Black lines opened'); pause
for k=1:3; In(:,:,k) = I(:,:,k) .* uint8(im_black); end;
image(uint8(In)); pause;
imagesc(bwmorph(im_black,'skel',Inf));

% 16,17,18 green area
im_green = (in == 16) | (in == 17) | (in == 18);
imagesc(im_green); title('Grean area'); pause;
im_green = imopen(im_green,strel('diamond',1));
imagesc(im_green); title('Green area opened'); pause;
im_green = imclose(im_green,strel('disk',7));
imagesc(im_green); title('Green area opened and closed'); pause;
for k=1:3; In(:,:,k) = I(:,:,k) .* uint8(im_green); end;
image(uint8(In)); pause;


Costas loop (and other PLL-based synchronization algorithms) tend to lose lock when the frequency offset is higher than the pull-in range (i.e. maximum frequency offset the said algorithm can cope with). Thus, a coarse frequency recovery algorithm is typically included in practice. The general practice is to put the coarse frequency recovery before the Costas loop. There are many ways of doing coarse frequency recovery including FFT-based algorithms and the frequency-locked loop (FLL) (described in Harris paper. My favorite is FFT due to its generality. The FFT algorithm can be found on the GOES satellite user manual and is summarized in the figure below.

The algorithm works by:

Carrier detection by absolute-law (residual carrier), square-law (BPSK) or fourth-power law (QPSK). Square law detector involves squaring the IQ signal while fourth-power law uses the fourth-power of the IQ signal.
Finding the FFT
Averaging successive FFT sequences
Finding the peak by looking at 8 adjacent bins with the highest magnitude.

Below is a GNU Radio example for BPSK (square law). As you can see, an FFT is included at the very beginning to reduce the frequency offset (in my case it was Doppler shift/rate + carrier oscillator offset).

If you are interested, you could have a look at a paper I wrote a couple of years ago. The paper covers all aspects of synchronization including phase, frequency, and timing.
I hope this helps. Good luck.


SIS Particle Filter Python Implementation

Particle filter is a Bayesian filter. Any Bayesian filter requires process and measurement models so you also need to define them.
If you will use bootstrap particle filter, you just create initial samples with initial parameters from Gaussian distribution. Then, you propagate the particles using process model. This is your so called "prediction" step.
After obtaining your measurements for each time step, you calculate the weights. Using those weights you perform resampling.


The input to the system $x_1[n]=\delta[n-5]$ is zero everywhere and is only one when the argument of delta is zero (at $n=5$). The output $$y_1[n]=x_1[5n+2]=\delta[(5n+2)-5]$$ would be zero everywhere too, except when the argument of $\delta$ is zero: $$5n+2-5=0$$ which requires $n=\frac{3}{5}$. But since $n$ only has integer values, it can never be achieved. Hence, the output of the system is zero always.


How to calculate autocorrelation function of an image noise

If you only care about the image noise you can get negative values.
Besides, I don't understand why you want to stay between 0 and 1, white noise with STD of 1 cannot exist in this range.


Mathematics of Harris corner point detection

The meaning of that formula is really quite simple. Imagine you take two same-sized small areas of an image, the blue one and the red one:

The window function equals 0 outside the red rectangle (for simplicity, we can assume the window is simply constant within the red rectangle). So the window function selects which pixels you want to look at and assigns relative weights to each pixel. (Most common is the Gaussian window, because it's rotationally symmetric, efficient to calculate and emphasizes the pixels near the center of the window.) The blue rectangle is shifted by (u,v).

Next you calculate the sum of squared difference between the image parts marked red and blue, i.e. you subtract them pixel by pixel, square the difference and sum up the result (assuming, for simplicity that the window = 1 in the area we're looking at). This gives you one number for every possible (u,v) -> E(u,v).

Let's see what happens if we calculate that for different values of u/v:

First keep v=0:

This should be no surprise: The difference between the image parts is lowest when the offset (u,v) between them is 0. As you increase the distance between the two patches, the sum of squared differences also increases.

Keeping u=0:

The plot looks similar, but the sum of squared differences between the two image parts is a lot smaller when you shift the blue rectangle in the direction of the edge.

A full plot of E(u,v) looks like this:

The plot looks a bit like a "canyon": There's only a small difference if you shift the image in the direction of the canyon. That's because this image patch has a dominant (vertical) orientation.

We can do the same for a different image patch:

Here, the plot of E(u,v) looks different:

No matter which way you shift the patch, it always looks different.

So the shape of the function E(u,v) tells us something about the image patch

if E(u,v) is near 0 everywhere, there is no texture in the image patch you're looking at
if E(u,v) is "canyon-shaped", the patch has a dominant orientation (this could be an edge or a texture)
if E(u,v) is "cone-shaped", the patch has texture, but no dominant orientation. That's the kind of patch a corner-detector is looking for.

  Many references say it is the magnitude by which the window 'w' shifted...so how much is the window shifted?one pixel...two pixels?

Normally, you don't calculate E(u,v) at all. You're only interested in the shape of it in the neighborhood of (u,v)=(0,0). So you just want the Taylor expansion of E(u,v) near (0,0), which completely describes the "shape" of it.

  Is the summation over the pixel positions covered by the window?

Mathematically speaking, it's more elegant to let the summation range over all pixels. Practically speaking there's no point in summing pixels where the window is 0.


Making a low pass filter for irregular samples

I realize that this is not a filter, but I do not believe any filter will give you a reasonable result given the strange sample locations. This seems to be borne out by several references on the web. This one for example and this PDF one.

I suggest something something along these lines, using projections onto convex sets (POCS).

Algorithm Descroption

Assume that the known red pixel values $R_i$ are at $(x_i,y_i)$ so that your start red image is:
$$
I_{\tt red}(x,y) = \sum_{i=1}^N R_i \delta(x-x_i, y-y_i)
$$

Set $k=0$, and initialize your estimate of the fully sampled red image with $I_{\tt red}(x,y)$:
$$
\hat{I}^k_{\tt red}(x,y) = I_{\tt red}(x,y).
$$

Take the 2D FFT of $\hat{I}^k_{\tt red}$:
$$
F_{\tt red}(n,m) = {\tt FFT}\{\hat{I}^k_{\tt red}(x,y)\}
$$
Because the red channel is slowly varying, remove all frequencies above a cutoff:
$$
F^{lp}_{\tt red}(n,m) = F_{\tt red}(n,m) L(n,m)
$$
where $L$ is a low pass mask.

Form the inverse FFT to get $\hat{I}^{k+1}_{\tt red}(x,y)$
$$
\hat{I}^{k+1}_{\tt red}(x,y) = {\tt FFT}^{-1}\{F^{lp}_{\tt red}(n,m)\}
$$

Set the known pixels of $\hat{I}^{k+1}_{\tt red}(x,y)$ from $I_{\tt red}(x,y)$.

If the difference between $\hat{I}^{k+1}_{\tt red}(x,y)$ and $\hat{I}^{k}_{\tt red}(x,y)$ is small enough, stop, otherwise reiterate with $k=k+1$.

Results

An example run with just a $6 \times 6$ matrix is below. The left image shows the original pixels. The right image shows the POCS version.  The red dots show the original pixel locations in both images.

Different Initial Conditions

I tried setting the initial image to

Ired <- array(100000,c(ImageSize, ImageSize))

and

Ired <- array(128,c(ImageSize, ImageSize))

and the results are below, so I do not believe the initial values for the unknown pixels matter.

The plot scales the greyscale values, which is why the 100,000 example looks wrong in the left-hand plot.

R Code Below

# 32018
# Non-Bayer 2D interpolation

# Model pixels as:
# P(x,y) = (a_0 + a_1 * x + a_2 * x^2) * (b_0 + b_1 * y + b_2 * y^2)
#        = [a_0 a_1 a_2 b_0 b_1 b_2] * [1 x x^2 1 y y^2]^T
# And so we need to find the a_i and b_i from specific pixel and coordinate values.
#
# install.packages('ptw')

require(ptw)

pad <- function(x,N)
{
  p1 <- ptw::padzeros(x,N,side="right")
  p2 <- t(ptw::padzeros(t(p1),N,side="right"))

  return(p2)
}

Npoints <- 8
xy <- array(c(c(1,1), c(1,5), c(2,3), c(3,6), c(4,2), c(4,4), c(5,6), c(6,3)), c(2,Npoints))
xy_norm <- (xy - min(xy))/(max(xy) - min(xy))
Pxy <- array(runif(N, 0, 255), c(1,Npoints))

ImageSize <- max(xy)

Ired <- array(0,c(ImageSize, ImageSize))

for (idx in 1:Npoints)
{
  Ired[xy[1,idx], xy[2,idx]] = Pxy[idx]
}

Ihatred <-array(0,c(ImageSize,ImageSize))
Ihatred2 <- Ired
count <- 0

while (sum(abs(Ihatred - Ihatred2)) > 1)
{
  Ihatred <- Ihatred2

  Fred <- fft(pad(Ihatred, 2*ImageSize))

  Fredsize <- dim(Fred)[1]
  cutoff1 <- Fredsize / 8
  cutoff2 <- 7 * Fredsize / 8

  Fred[4:16,] <- 0
  Fred[,4:16] <- 0

  IFred <- fft(Fred, inverse=TRUE) / Fredsize / Fredsize

  Ihatred2 <- Re(IFred[1:ImageSize, 1:ImageSize])

  for (idx in 1:Npoints)
  {
    Ihatred2[xy[1,idx], xy[2,idx]] = Pxy[idx]
  }

  count <- count + 1
  if (count > 100)
  {
    break
  }
}

par(mfrow=c(1,2))
image(Ired, col= grey(seq(0, 1, length = 256)))
for (idx in 1:Npoints)
{
  points(xy_norm[1,idx] , xy_norm[2,idx] , col="red", lwd=5)
}
title('Original with only known red pixels')
image(Ihatred2, col= grey(seq(0, 1, length = 256)))
title('POCS solution')
for (idx in 1:Npoints)
{
  points(xy_norm[1,idx] , xy_norm[2,idx] , col="red", lwd=5)
}


DSP is a domain mixed between EE/CE and Math/CS. Hardware must be built, but also software written.

Some CS programs do touch DSP at the software level. An applied math program would touch even more the algorithmic side, since pretty much all DSP techniques are based on certain mathematical areas such as Fourier analysis, Complex analysis and Numerical analysis. An EE program might have a specialisation on DSP, covering both the hardware and software sides but at a more "applied", less "proof mathematical" level, than what mathematics programs do.

The field is that wide though that it might also be better to focus on only one thing and do it well, rather than think in a "jack of all trades" way, i.e. too broadly. That is, if your interest is in software and mathematics, then focus on those. If your interest is in hardware, then focus on that. Applied math is sort of a middle ground. It's both theoretical, but it will lack practice such as electronics design, even if it'd give you conceptual tools for understanding electromagnetism and digital logic.


The main issue is that the filters start with an initial position of (0,0). Then when you filter your actual data, the filter output has to ramp up to the actual position. In doing so, you can also observe some overshoot since your filter has a very short time response.

It would help if you were to start with a better estimate of the initial position by using for example the first position with:

dataOut_x = xcor_i(1)+filter(b,a,xcor_i-xcor_i(1));
dataOut_y = ycor_i(1)+filter(b,a,ycor_i-ycor_i(1));

The resulting plot would then look like:


Low pass filter output showing instability

I am going to explain Robert's answer below.

Apply Z-transform to both sides of the equation. You get the following:
\begin{align}
Y(z) [1 - 2z^{-1} + z^{-2}] = X(z)[1 -2z^{-6} + z^{-12}]
\end{align}

The transfer function is given as
\begin{align}
H(z) = \frac{Y(z)}{X(z)} = \frac{1 - 2z^{-1} + z^{}-2}{1 - 2z^{-6} + z^{-12}} = \frac{(1-z^{-6})^{2}}{(1-z^{-1})^{2}}
\end{align}

The denominator represents the poles of the filter $H(z)$. In this case, there exists two poles, both of which lie on the unit circle at $z=1$. For a filter to be stable, all the poles should lie inside the unit circle, i.e., the roots of the denominator polynomial of the transfer function $H(z)$ should lie inside the unit circle.

Perhaps, you can also try a different MATLAB implementation of the above filter:

a1 = [ 1 −2 1 ];
b1 = [ 1 0 0 0 0 0 −2 0 0 0 0 0 1 ];
lp_sig = filter( b1, a1, x );


How to get around the circular shift property of Discrete Fourier Transform?

Sorry No.

Circular shift property of the DFT (or actually the DFS, @robertbristow-johnson will love this!) is established by law; you cannot get away from it using other clever techniques...

May be you can introduce some redundancies (such as long set of samples but short windows on them, i.e., zero padded signals) you can do some tricks. But at a fundamental level, DFT will always perform a circular shift.


Designing an allpass IIR filter with arbitrary phase response

Depends a bit on what exactly you are trying to do. Generally, strict all-pass filters have reciprocal zeros and poles. This will always result in 180 phase shift per filter order so an arbitrary phase response that's strictly allpass cannot be designed. Example: let's say you want something that's in phase at DC,  sits at -90 degrees in the band and then goes slowly back to zero at Nyquist: this can't be implemented using strict all-passes.

Of course you can always approximate but the best method tends to depend on the specifics of your design problems and what types of errors you can tolerate.


Generating smoothed versions of square wave, triangular, etc

It seems the amplitude is not scaled properly. Rather than (2*A/pi) using (A/atan(1/delta)) seems more appropriate. In other words I propose:

y = (A/atan(1/delta))*atan(sin(2*pi*t*f)/delta);

Below is a figure illustrating the difference between the two scaling approaches. For low delta values the difference is not clear but for high delta values the amplitude departs from what is desired. The images can be reproduced using the MATLAB code below.

t=linspace(0,2*pi,500);
delta = 0.1;
A = 1;
f = 1/(2*pi);
y1 = A*sin(2*pi*t*f);
y2 = (2*A/pi)*atan(sin(2*pi*t*f)/delta);
y3 = (A/atan(1/delta))*atan(sin(2*pi*t*f)/delta);

figure; hold on;
h1=plot(t,y1,'r-','LineWidth',5);
h2=plot(t,y2,'g-','LineWidth',4);
h3=plot(t,y3,'b-','LineWidth',2);
legend([h1 h2 h3],{'Sine','Smooth square wave','Corrected smooth square wave'},'Location','SouthOutSide');
axis tight; axis square;
set(gca,'FontSize',15); grid on; box on;
drawnow;


Normalizing white noise to match pink noise's amplitude at a given frequency?

Okay, when discussing white noise or pink noise (or red noise or brown noise or flicker) or some other random process, there is this property called the power spectrum, in which white noise has a constant value for all frequencies.  But we integrate the power spectrum over all frequencies (from $-\infty$ to $+\infty$) to get power.  Integrating a constant, non-zero function over an infinite frequency range gives you infinite power.

So in reality, white noise does not exist.  (It puts out more power than the Sun.)

But if you bandlimit white noise to a finite bandwidth, you have a finite power.

Sampled signals (in discrete time) must have a finite bandwidth in order to be sampled properly.

If you were to synthesize bandlimited white noise by use of a random number generator that produces good random numbers with a uniform distribution between the values of $0$ and $\Delta$, you would have a mean of $\frac{\Delta}{2}$ and a variance of $\frac{\Delta^2}{12}$.  If you subtract from each uniform random number the mean of $\frac{\Delta}{2}$, then you would have the same variance, but a mean of zero.

Now if you were to add $N$ of those independent zero-mean uniform random numbers together, the variance would increase to $\frac{N \Delta^2}{12} \triangleq \sigma^2 $ and the distribution of the random numbers would closely approach that of a normal distribution (sometimes called Gaussian).

Now variance of a zero mean (or DC free) random process is the power, which is the area (or integral) of the power spectrum.  So from $-\frac{f_\mathrm{s}}{2}$ to $+\frac{f_\mathrm{s}}{2}$, you would have a flat power spectrum with an amplitude of $\frac{\sigma^2}{f_\mathrm{s}}$ so that the area or integral comes out to be $\sigma^2$.  (The amplitude in the power spectrum is zero for all $|f|>\frac{f_\mathrm{s}}{2}$.  This is what makes it bandlimited.)

So now you have "white" noise that is finite energy.  When you run that through a pinking filter (is that what you are doing to make pink noise?), then you have to multiply that constant amplitude of $\frac{\sigma^2}{f_\mathrm{s}}$ times the magnitude-squared of the pinking filter to get the frequency-dependent amplitude of the pink noise.  Then you have to integrate that from from $-\frac{f_\mathrm{s}}{2}$ to $+\frac{f_\mathrm{s}}{2}$ to get the power of the pink noise.  Whatever is the ratio of that power to the power of the bandlimited white noise going in, you compute the square-root of the reciprocal, and that is the scaling constant to multiply times your pink noise to make it have the same power as the bandlimited white noise going in.


The division should be by 100000, but I doubt that hardware is giving the result in actual milliwatts, maybe a normalization is required in the power being reported by the hardware as well followed by a subsequent normalization in bandwidth. This would require knowing how the hardware is generating this number,

considerations like, what is amount of duration over which hardware measures this spectrum, does it average out over those durations or does it provide cumulative value which needs to be normalized.


What is this effect/filter in given images?

Highlight blooming; the bright parts are bleeding into the shadows. It can be done with Photoshop's Diffuse Glow filter.


I have a noisy ECG Signal. want to use PSD to analyse the Signal to noise ratio

Why it is plotting against only 500hz and what information I can get on this??

Plotting against $500Hz$: Because that is the frequency of the sampled signal. We usually sample a signal at least two times of its max frequency in order to avoid aliasing. In your case $fs=1000Hz$ (sampling frequency), so max frequency is (no higher than) $500Hz$.
The information you can get: Which frequency corresponds to high powers. You can avoid power below $0$ frequencies. High frequency components are usually our signal while low power is noise.


Why Does Image Deconvolution Still Work with Image without Sharp Edges?

Think of it in the frequency domain.
Some of the data is highly damped.
Namely energy in High Frequency is damped with low values.

Yet, unless it is zero identically "There is hope...".

Basically, this is the idea behind deconvolution algorithms.
Restore the energy in frequencies before it was damped.

Practically, we are limited in DR.
Namely if it is damped really hard, even if it not zero, in order to restore it one would amplify it significantly which will also lead to noise.

Regarding the edges, well, assuming we have the right way to recover the enrgies edges will be restored.


Erosion compotition rule $\left ( f \ominus g \right ) \ominus h = f \ominus \left (g \oplus h \right ) $

Morphological Erosion is not associative as dilation. Much like addition is associative but not subtraction. Consider the following basic definitions first:$$(A \ominus B)^C = A^C \oplus \widetilde{B} \tag{1}$$
$$(A \oplus B)^C = A^C \ominus \widetilde{B} \tag{2}$$

where $\widetilde{B}$ is the reflection of kernel $B$ and $A^C$ is the compliment of image $A$, inverse of image $A$. Now let us evaluate the expression in the question:
$$(A \ominus B) \ominus C = ((A\ominus B)^C \oplus \widetilde{C})^C \tag{using eq. 1}$$

$$(A \ominus B) \ominus C = ((A^C\oplus \widetilde{B}) \oplus \widetilde{C})^C \tag{using eq. 1} $$

$$(A \ominus B) \ominus C = (A^C \oplus (\widetilde{B} \oplus \widetilde{C}))^C \tag{since dilation is associative}$$

$$(A \ominus B) \ominus C = A \ominus \widetilde{(\widetilde{B} \oplus \widetilde{C})} \tag{using eq. 2}$$

$$(A \ominus B) \ominus C = A \ominus (B \oplus C) \tag{distributing the reflection}$$

Hope that answers the question


What Is the Walsh Hadamard Transform? How Could One Use It in Image Processing?

NASA used to use the Hadamard transform as a basis for compressing photographs from interplanetary probes during the 1960's and early '70s.  Hadamard is a computationally simpler substitute for the Fourier transform, since it requires no multiplication or division operations (all factors are plus or minus one).  Multiply and divide operations were extremely time intensive on the small computers used on board those spacecraft, so avoiding them was beneficial both in terms of compute time and energy consumption.  But since the development of faster computers incorporating single-cycle multipliers, and perfection of newer algorithms such as the Fast Fourier Transform, as well as the development of JPEG, MPEG, and other image compression, I believe Hadamard has fallen out of use.  However, I understand it may be staging a comeback for use in quantum computing.
(NASA use is from an old article in NASA Tech Briefs; exact attribution unavailable.)


How do I find Gray Level Variance of the Image?

I would think that it is indeed the "normal" variance of the image.  You look at the distribution of the pixel values (i. e. gray levels), and compute its variance.


Testing whether a process is a Wiener process

That's not a good question to ask (you are ok and I'll vote up your question in a sec).
BM induces a measure on functions, called the wiener measure. Since this is a measure on functions, it is infinite dimensional. The thing with measures on infinite dimensional spaces is that they tend to be mutually singular: events that happen w.p. 1 wrt one measure may have zero probability wrt the other. For example, take BM $B_t$ and BM with some tiny constant added $C_t= B_t + \epsilon$.
Then $P(B_0 = 0)=1$ but $P(C_0=0)=0$. This situation is fairly generic.

What this means to you is that unless your process really is exactly BM (highly unlikely), as you take more and more data points your measure will converge to something that is mutually singular to BM and your hypothesis test will give a "no". This may happen even if your process is pretty "similar" to BM.

That's why this is not a good question to ask. Measures over function spaces are tough to work with.


You really didn't give sufficient context and background information to provide a complete answer, but I'll try anyway. You didn't say so, but I assume that $\sigma_W^2$ is the power of a watermark and that $\sigma_N^2$ is the power of additive noise. I furthermore assume that the "usable bandwidth" is the bandwidth for generating the watermark. If we additionally assume that the watermark signal is white with spectral density $W_0$, then its power inside the usable bandwidth $B$ is given by

$$\sigma_W^2=2W_0B$$

If $B_1=22050Hz$ and $B_2=4000Hz$ we get for the capacity ratio

$$\frac{\log_2(1+\frac{2W_0B_2}{\sigma_N^2})}{\log_2(1+\frac{2W_0B_1}{\sigma_N^2})}$$

which is not a fixed ratio, but it depends on the constants $W_0$ and $\sigma_N^2$. If they are not given, then there is no way to claim that this ratio is $1/5$.

However, if you can assume that $\sigma_W^2/\sigma_N^2\ll 1$ then the capacity can be approximated by

$$\frac{1}{2\log (2)}\frac{\sigma_W^2}{\sigma_N^2}=\frac{1}{2\log (2)}\frac{2W_0B}{\sigma_N^2}$$

and with the same assumptions as above the capacity ratio is equal to the ratio of bandwidths $B_2/B_1=0.18$, which is indeed close to $1/5$.

If any of my assumptions are not correct, please clarify your question.


MATLAB code to evaluate audio loundness over time

Any audio recording is generally assumed to have units that are proportional to Pascals (Pa).  This is the instantaneous pressure of the sound wave above (positive) or below (negative) the ambient pressure.

(To actually get Pascals from the WAV file, you need to know the calibration factor for your recording system.)

Loudness of sound (in air) is generally reported as Sound Pressure Level (SPL) in decibels. SPL (in dB) and pascals are related as

Pa_ref = 20e-6;  %reference pressure for SPL
SPL_dB = 10*log10( (Pa / Pa_ref).^2 );

If you use the instantaneous pressure values (ie, the values straight from the WAV file), this equation will give you the instantaneous SPL.  Based on the graph that you showed, you don't want the instantaneous SPL.  You want a time-smoothed SPL.

To get a time-smoothed SPL, the standard approach is to extract the envelope of the pressure signal.  There are lots of ways of doing this.  One way is to rectify and low-pass filter the signal.  In this case, your processing might look like this:

%load and calibrate the data
[wav,fs_Hz]=audioread(myfilename);  %load the WAV file
my_cal_factor = 1.0;  %the value for your system to convert the WAV into Pascals
wav_Pa = wav * my_cal_factor;

%extract the envelope
smooth_sec = 0.125;  %"FAST" SPL is 1/8th of second.  "SLOW" is 1 second;
smooth_Hz = 1/smooth_sec;
[b,a]=butter(1,smooth_Hz/(fs_Hz/2),'low');  %design a Low-pass filter
wav_env_Pa = sqrt(filter(b,a,wav_Pa.^2));  %rectify, by squaring, and low-pass filter

%compute SPL
Pa_ref = 20e-6;  %reference pressure for SPL in Air
SPL_dB = 10.0*log10( (wav_env_Pa ./ Pa_ref).^2 ); % 10*log10 because signal is squared

%plot results
figure;
subplot(2,1,1);
t_sec = ([1:size(wav_Pa)]-1)/fs_Hz;
plot(t_sec,wav_Pa);
xlabel('Time (sec)');
ylabel('Pressure (Pa)');

subplot(2,1,2)
plot(t_sec,SPL_dB);
xlabel('Time (sec)');
ylabel('SPL (dB)');
yl=ylim;ylim(yl(2)+[-80 0]);


How to circularly shift a signal by a fraction of a sample?

If you want the shifted output of the IFFT to be real, the phase twist/rotation in the frequency domain has to be conjugate symmetric, as well as the data.  This can be accomplished by adding an appropriate offset to your complex exp()'s exponent, for the given phase slope, so that the phase of the upper (or negative) half, modulo 2 Pi, mirrors the lower half in the FFT aperture.  The complex exponential shift function can also be made conjugate symmetric by indexing it from -N/2 to N/2 with a phase of zero at index 0.

It just so happens that the appropriate offset for phase twists or spirals, that complete an exact integer multiples of 2 Pi rotations in aperture, to be conjugate symmetric in aperture, is zero.

With a conjugate symmetric phase twist vector, the result should then end up as a circular Sinc interpolation for non-integer shifts.

Elaboration by OP:

Your choice of k = [0, 1, 2, 3, 4, 5, 6, 7, 8] is producing an asymmetrical complex exponential:

If you use k = [0, 1, 2, 3, 4, -4, -3, -2, -1] instead, you get a Hermite-symmetric complex exponential:

plot(fftshift(exp(-1j * 2*pi * 0.5/N * k)))

and now when you use the same exponential formula to shift by 0.5 or 3.5 samples, you get a real result:

plot ifft(fft(a)*exp(-1j * 2 * pi * 0.5/N *k))
plot ifft(fft(a)*exp(-1j * 2 * pi * 3.5/N *k))


response of algorithm in non-equidistant time

I would try to stick with a conventional discrete time control approach.  Trying to deal with missing samples in the control algorithm complicates things and without dealing with them the control output will not be very stable.  Can something like this work in your application?

Here the assumption is that you can detect missing samples and use an interpolation algorithm to fill in the gaps before applying the signal to your controller/compensator.  The "Control Signal" input is a reference that can be used to set the normal operating intensity of the Light Source.

I would characterize the transfer function with using normal control signals (impulse, step, swept sine, all with uniform sampling).  Then I would model the real signal (from the lamp transducer) with some impulse variations that match the characteristics of your missing samples.  By impulse variations I mean adding small error quatities to signal to represent error from your interpolator.  Test how well the transfer function regulates in the presense of this "impulse noise".


How to quantify delay in filtered output?

It's the transient and the group delay which are associated with the initial conditions and the LTI filter phase response $\phi(w)$ respectively.

The group delay (in samples) associated with the LTI filter is:
$$\tau = - \frac{d \phi(w) }{dw} $$

The output will be shifted by $\tau$ samples wrt to the input. You shall compansate this shift for synchronization.

This delay is constant for linear phase filters which are of symmetric FIR type filters. For IIR type filters, since their phase response is not linear therefore getting a fixed group delay is not possible.


How to explain this BER vs SNR behaviour?

In a BPSK system, the SNR at the output of the receiver's matched filter is $$\text{SNR}=\frac{E_b}{\sigma^2_n}.$$ In this formula, $E_b$ is the energy used to transmit one bit, and $\sigma_n^2$ is the power of the noise at the filter output.

Assuming an AWGN channel, the noise process at the filter input has constant PSD; this value is usually denoted $N_0/2$. The Wiener-Khinchine theorem relates the PSD and $\sigma_n^2$: $$\sigma_n^2=\frac{N_0}2,$$ so the SNR may be written as $$\text{SNR}=\frac{2E_b}{N_0}.$$

Note that the convention is to plot the bit error rate vs $E_b/N_0$ (note the missing $2$ in the numerator). The reason is that for quadrature systems, the total noise power in the receiver is $N_0$.

In a coded system, each information bit is transmitted using several coded bits. For example, a Hamming (7,4) code uses 7 coded bits to transmit 4 information bits; the code rate is $R=4/7$. In gneral, $R=k/n$, for a code that transmits $k$ information bits using $n$ coded bits.

To make a fair comparison between different codes, we usually keep $E_b$ constant. This means that the energy of each code bit depends on the code rate: $$E_c=RE_b,$$ or $E_b=E_c/R$.

In simulation, it is often convenient to set $E_c=1$. Then, we have $$\frac{E_b}{N_0}=\frac{E_b}{2\sigma_n^2}=\frac{E_c}{2R\sigma_n^2}=\frac{1} {2R\sigma_n^2}.$$


What are the uses of those three types of wavelet transformations?

The denomination you propose (Continuous wavelet transform, Discrete wavelet transform, Redundant wavelet transform) might be not precise enough. Yet here are the major differences.

Continuous wavelet transform: per se, they are a pure mathematical tool. You can use almost any wavelet you wish, with very little requirements. There exist discretizations, and they could be very redundant (roughly $10-100\times$ the original signal per dimension). You mostly use it on 1D signals due to the redundancy, but it is sometimes used in 2D and 3D with limited scales. It is used when one looks as fine details of the signals: fractal aspects, measures of regularity, precise onset or rupture detection in the signal's shape (continuity, changes in derivatives), adaptive filtering. One can just observe (analysis) or detect on 2D scalograms, or select an inverse transform (synthesis,), for instance for adaptive filtering.

Discrete wavelet transform: one only have a limited number of wavelets available. Often used in orthogonal or biorthogonal fashion, edp' when critical sampling or orthogonality is important: compression, statistical analysis, or when memory or computation is too expensive. It may induce wavelet-related artifacts with low-amplitude signals. Requires some technical skills to be used in detection.

Redundant wavelet transform: a mix between the two above. Limited choice in the wavelet, yet intermediate redundancy (roughly $1,5-10\times$ per dimension) and more important some invariance: translation, rotation, shear invariance, that is useful for detection, and often more robust to noises. Often used in restoration (deblurring, deconvolution, denoising),  feature analysis.


From the given transfer function you can see that the filter has a low pass character. The gain at DC is given by $H(1)=1$. This is also the maximum value of $|H(e^{j\omega})|$. So you can try to solve
$$\big|H(e^{j\omega_c})\big|^2=\frac{\max_{\omega}|H(e^{j\omega})\big|}{2}=\frac12\tag{1}$$
because this equation defines the frequency $\omega_c$ at which the filter gain equals $1/\sqrt{2}$, which is the definition of the $3$ dB cut-off frequency.
Your calculation of the magnitude of $H(e^{j\omega})$ is correct, so you can just continue by solving $(1)$.
As a final note, there are values of $a$ for which the $3$ dB cut-off frequency doesn't exist because the filter never reaches an attenuation of $3$ dB. You can try to find that limit of $a$ as a bonus exercise.


Fuse two sources of linear acceleration with a Kalman filter

That depends on how many of the accelerometer parameters (mostly drift and misalignment) you're trying to estimate.

If the IMU and the 'extra' accelerometer were in perfect alignment (and if their statistics are Gaussian), then the optimal combination of their outputs would be a simple weighted sum: $\vec {\hat a} = k_1 \vec a_1 + k_2 \vec a_2$ where (assuming equal scaling), $k_1 = \frac{\sigma_2^2}{\sigma_1^2 + \sigma_2^2}$, $k_2 = 1 - k_1$, and $\sigma_1$ and $\sigma_2$ are the noise variance of each accelerometer.

If you look at that first equation, you'll see that if one accelerometer is markedly quieter than the other, then unless it has other problems (misalignment, scaling, initial offset, etc.) your best bet is to just ignore the noisy accelerometer altogether.  Even at the point where one accelerometer has three times less variance than the other, the ultimate contribution of the "noisy" accelerometer ends up being $1/9^{th}$ of the total -- which is definitely tending into "why bother?" territory.

If you've got some oddball combination, i.e. you've got a really quiet accelerometer with terrible initial offset or misalignment vs. a really noisy one with low offset and misalignment, then life gets more complicated -- in fact, I'd have to ponder on exactly how you'd express that in a way that would make a sensible Kalman filter.


Difference between two Sine wave equations

The first equation:

x(t) = Sin(2*pi*F*t)

is for a continuous ("analogue") signal, where t may take any real value.

The second equation:

x(t) = Sin(2*pi*F*n*Ts)

(I assume that this is the correct version - it looks like you have a typo in the question ?)

is for a discrete ("sampled" or "digital") signal, where n is the sample number (integer) and Ts is the sample interval (inverse of sample rate). It has values only at the discrete sampling points where t = n * Ts.

In your MATLAB code you are effectively sampling a continuous function, so although you start off with the first (continuous) equation you end up with the second equation because you are only evaluating (sampling) the continuous function at a set of discrete points.

Note also that as Ts -> 0 the discrete version tends towards the continuous version.


How to calculate the total signal strength and phase shift of multiple signals in a reciever

This is simplified assuming there is not multipath occurring in the transmission (where the receiver would typically receive multiple copies of the same transmit signal at different delays resulting in multipath distortion which are typically addressed using channel estimation and equalization).

The OP's phase shift would all be relative since there is no common reference (clock) with the transmitter, so the comparison would be of the phase of the received signals to each other. A cross-correlation of the received signals in complex notation with the reference signal of what you know was transmitted would provide this.

This is accomplished by doing a multiply and accumulate of the received signal with the complex conjugate of the known ideal transmit signal, repeating the computation below for each possible offset $m$  (thus is the cross-correlation function specifically, correlating at each possible offset in time).

$$XCorr[m] = \sum r[n-m]t^*[n]$$

Where (*) indicates a complex conjugate.

The magnitude of the correlation peak would be proportional to the signal strength, and the shift $m$ is proportional to the delay. With two signals at the same shift $m$ (in the same bin), the phase of the correlation is proportional to the carrier phase, giving a fine delay indication by comparing the phase result from two different receiver signals if they were in the same bin. In this case the signals would need to not be correlated themselves, otherwise the result would be the vector addition of the two (appearing as a single received signal). For this reason the transmit signals from each transmitter would be uncorrelated by separating them in frequency, time or code; and you can identify each one separately in receiver using correlation.

Note the above computation for the cross correlation function can be done directly with FFTs given the FFT correlation property (which results in a circular cross-correlation):

$$XCorr[m] = \text{ifft}\{\text{fft}(r[n])\}\text{fft}^*\{(t[n])\}$$


Recorded an audio file using microphone via Windows Sound recorder, but has 2 channels that are slightly different. Why are they different?

I am not %100 sure, but it sounds like you are recording the file in stereo. If this is the case, then Windows Sound recorder is taking stereo input, which means it's accessing 2 AD converters which have the same input applied to them. The difference between the two signals in this case would be the noise and any offset between the two converters.

I would see if Windows Sound Recorder could make a mono recording, or else use another program to make the recording.


Looking for cycles of periods longer than the input signal length

If you plot frequency instead of period, you might be able to see that some of the increasing artifacts are due to the shrinking rectangular window on your signal.

An FFT decomposes the entire length of its input into sinusoids.  You have created a flat spot in your data input, and it takes a linear combination of a lot of sinusoids to represent those extended flat spots.  So the FFT of your latter shorter signal segments becomes more and more about the flat spot than your signal.

Added: There are other methods for estimating the parameters of a known small number of sinusoids within a known distribution of noise, such as MUSIC and ESPIRIT, that may be more suitable for estimation of sine waves not equal to or possibly shorter than the lowest DFT basis vectors.


Analytical way to generate an IIR that is the square root response of another IIR

The only time you'll be able to do this exactly is when you know that the poles and zeros of $H(z)$ are in complex conjugate pairs inside and outside of the unit circle.

Note that for systems with real-valued impulse responses, this means you need four of each: two inside and two outside the unit circle.

The keyword to search on is spectral factorization.

This means that if you can write $H(z)$ as:
$$
H(z) = G(z) G(z^{-1})
$$
where $G(z)$ is a minimum phase system, then $G$ just takes all the poles and zeros of $H$ that are inside the unit circle (and half those that are on the unit circle).

There exist approximate algorithms for finding a solution, but I am not well-versed in those.


Why do we need to construct Gaussian pyramid using SIFT detector

Actually, the purpose of all this is to approximate a Laplacian of Gaussian!
This computation is part of the corner detection of SIFT. You can find corners by examining extrema of the Laplacian of Gaussians (2nd order derivative). You use Gaussians for denoising, and a Laplacian to find inflection points. However, it is classical to not deal directly with computing this second order derivative, but to approximate it instead by difference of Gaussians for simplicity. This is sometimnes referred to as Gaussian blob detection.
You can find some illustration of this process here.
With this approach your difference of Gaussians (DoG) will have a fixed (small) size in pixels, so you're bounding the computational cost. But you still want to detect blobs of different sizes (because objects have varyng size,s, and the camera can be closer or further away from the scene), and you achieve this by "zooming out the image" with the various scales of the pyramid and keeping the DoG size fixed, instead of keeping the original image size and changing the size of your blob detector.


Help with denoising signal and periodogram analysis resources

Firstly, I am confused if I am supposed to filter my signals to get rid of any frequencies above the Nyquist frequency. My sampling frequency is 32Hz and my time series is somewhat noisy and has some artifacts. I am also unsure of which filter to select for this.

That ship has sailed.
Let $S=\left\{\left.\alpha e ^{i(\omega t+\varphi)}\right|\alpha > 0, 0\le \varphi <2\pi, \omega \in \mathbb R\right\}$, i.e. the set of all distinct complex sinusoid with an amplitude, frequency and phase.
Then $(S,\cdot)$ is a commutative semigroup (proof trivial).
Introducing the equivalence relation $\sim: a\sim b \iff a\left(\frac{n}{r}\right)= b\left(\frac{n}{r}\right) \forall n\in\mathbb Z$ ("two signals are identical after sampling with rate $r$"), we see that signals $s_l=\alpha e^{i(\omega_l t + \varphi)}, l=1,2,\ldots$ are $s_1\sim s_2$ if $\frac{\omega_1-\omega_2}{2\pi}=nr, n\in\mathbb Z$, i.e. we can't tell signals apart after sampling if their frequency differed by a multiple of the sampling rate.
Let's formalize this: $T=\left\{\alpha e^{i([ft+\phi]\mod 2\pi)}\right\}$ is a quotient semigroup of $S$, i.e. $T\preceq S$, and each of the elements is a leader of a $\sim$-equivalence class – and this homomorpism $S\mapsto T$ is in fact sampling, which, as we can see above, is not bijective.
Hence, all the original frequency components from $S$ got mapped to some component from $T$ with frequency normalized to the sampling frequency. That mapping is called aliasing.
What an anti-alias filter $h$ does is
$$h(s): S\mapsto S, h=\begin{cases}
s, & f< f_\text{nyquist}\\
0 & \text{else,}
\end{cases}
$$
and as you'll figure out when inserting $h(s)$ above is that this yields only the elements that are not aliased to a different frequency.
Thus, everything that "survives" $h$ will also "survive" aliasing without undergoing a change in frequency.
So, if you needed an anti-aliasing filter, it's now too late. Go and do your recording again.

Secondly, I wanted to look at the periodogram and the average power spectral density at each frequency within a group - but I am not sure if I understand the periodogram very well - if I have different time series lengths then my periodogram length will vary too, so I am not sure how this comparison can be made.

In that case, the periodogramm is not a useful mapping on its own – you'll need to add something like a truncation / padding operation to bring all signals to the same duration, for example. At which point the periodogram doesn't seem to be a sensible approach anymore.


Filter design - Rules of thumb regarding ratio of sampling and stop/pass frequency

In the most general case, there isn't really a formula for how long the filter should be. There are papers which show the empirical relationships between the filter transition width, passband ripple, and stopband attenuation. You typically have - for a fixed sampling frequency:

A smaller transition width (from the stopband to the passband), as a percentage of the sampling frequency, requires a longer filter. A 100 Hz transition width for fsamp=5kHz would require a smaller filter than a 100 Hz transition width with fsamp=15kHz - assuming other requirements are equal.
Smaller ripples in the passband, requires a longer filter
More attenuation in the stopband requires a longer filter

That said, for remez type designs there are a couple of published formulas that give you an estimated length of the filter. In Matlab you could use the remezord() function. These formulas work fairly well, but they do run into problems when the transition frequencies are near 0 Hz or fs/2.

There are some published papers that give equations for the length of filter using a Kaiser window for a windowed Sinc design - but you don't have as much control in the specification of the filter - if I recall correctly.


how to evaluate derivative of convolution integral?

Note that option (b) is not correct, and that it is also not equal to what you came up with. Option (b) is just the multiplication of $x(t)$ and $y'(t)$, not the convolution. Your solution and option (c) are both correct, assuming that all derivatives exist and that the convolution integrals converge, because with that assumption the following holds:

$$\frac{d}{dt}\left(x(t)\star y(t)\right)=x'(t)\star y(t)=x(t)\star y'(t)\tag{1}$$

where $x'(t)$ and $y'(t)$ denote the derivatives of $x(t)$ and $y(t)$, respectively.

An ideal differentiator is just an LTI system, so it can be written as a convolution, and since convolution is associative, $(1)$ must hold.

This can also most easily be seen in the frequency domain. Differentiation corresponds to multiplication with $j\omega$, so we have

$$j\omega\left(X(\omega)Y(\omega)\right)=\left(j\omega X(\omega)\right)Y(\omega)=X(\omega)\left(j\omega Y(\omega)\right)\tag{2}$$

which is equivalent to $(1)$.


How to express time delayed discrete signals

The problem with representing a discrete time just as a list like: $$x[n]=[6, 4, 0, 2]$$ is that you don't have a time reference.

The first index $n=0$ or the sample in $\textbf{bold}$ or $\underline{\mathrm{underline}}$ is usually taken as the start of the sequence.

For example, $$x[n]=[\mathbf{6}, 4, 0, 2]$$ or $$x[n]=[\underline{6}, 4, 0, 2]$$

So $x[n-2]$ is a delay of two samples, and $x[0]$ now happens two samples later: $$x[n]=[\underline{0}, 0, 6, 4, 0, 2]$$


Is possible reach the DFT if I have the DTFT?

The DFT is a sampled version of the DTFT only for finite length signals. Otherwise, there is no point in comparing the DTFT with the DFT because you can only compute the DFT for finite length (or periodic) signals. Your example $x[n]=\sin(\omega_0n)$ is an infinitely long sequence. For specific choices of $\omega_0$ it is periodic, but nevertheless, you can't obtain its DFT by sampling its DTFT. The operation of sampling a Dirac impulse is undefined because $\delta(\omega)$ cannot be evaluated at $\omega=0$ since the Dirac impulse is not an ordinary function but a distribution.


Regarding the harmonics: ASK has an infinite number of them. The ASK signal can be thought of as a square wave multiplied by a sinusoid. In frequency, this corresponds to the spectrum of the square wave shifted to the frequency of the sinusoid. Since the spectrum of the square wave is infinite, the spectrum of the ASK signal is infinite too.

In practice, the ASK signal is filtered to force it to fit in a specific bandwidth. This has the effect of distorting the signal, but if done carefully the distortion is not severe. In any case, ASK is only used in applications where bandwidth is not a constraint.

I recommend the following procedure to understand how to generate an ASK signal. First you need to map "bits" to "voltages" (also called "amplitudes"). Bits are abstract quantities that do not have physical substance, while voltages are physical. In ASK you typically map a bit 0 to 0 volts, and a bit 1 to $A$ volts. Now you can generate a sequence of voltages that corresponds to your bits:

$$ \text{bits}=\lbrace 0,1,1,0,1,1,0,1 \rbrace$$

$$ \text{amplitudes} = \lbrace 0, A, A, 0, A, A, 0, A \rbrace $$

Now you need to define a pulse $p(t)$ of duration $T_b$:

$$ p(t) = \begin{cases}1,\quad0\leq t\leq T_b,\\0,\quad\text{otherwise}\end{cases} $$

(You will eventually want to control the energy of the pulse by changing its amplitude).

The next step is to define the "baseband" ASK (that is, with no carrier) signal:

$$ s_{BB}(t) = \sum_{k=0}^7 a_k p(t-kTb) $$ where $a_k$ is the $k$-th amplitude (starting from 0), and $p(t-kT_b)$ is the pulse shifted $kT_b$ seconds "to the right". Take a moment to study this equation and you'll see that it is a train of pulses of amplitude 0 for a bit 0 and amplitude A for a bit 1.

Now you're ready to create the ASK signal, which is:

$$ s(t) = s_{BB}(t)\sin(2\pi f_c t) $$

To transmit more bits per pulse, all you need to is to increase the range of values of the amplitudes.

I know that this is different from the way most undergrad textbooks teach this subject, but I still think it makes more sense to think about it this way.


The $\sigma$ parameter is both. The Gaussian function can generate a scale-space where $\sigma$ is the scale parameter. It doesn't mean the image is scaled, instead it is the scale at which the features are being evaluated.

For example, with higher $\sigma$ the image is more blurred and therefore only larger image features contribute to the gradient histogram. SIFT is made scale-invariant by choosing the greatest feature response from a set of responses from different scales. The Gaussian function is used because with increasing scale one can guarantee no new local extrema (features) are introduced.

I'm not sure about the interpolated extrema, but I guess it is to obtain sub-pixel accuracy.


A speech communication channel as used in telephony typically has a frequency response of 300 Hz to 3 kHz. Although this rejects a lot of the energy in normal speech, intelligibility is still quite good - the main problem seems to be that certain plosive consonants, e.g. "p" and "t", can be a little hard to discriminate without the higher frequency components.

So you're probably looking for a "sweet spot" somewhere between using the full 20 Hz - 20 kHz bandwidth typically found in consumer audio and the most aggressive filtering used for voice comms (see above). I would suggest starting with a bandpass filter from say 50 Hz to 8 kHz. It will probably only improve SNR by a few dB at best, but it may help, particularly if you have a lot of high frequency background noise.


I have come across my ways.

$$\begin{align} P(e^{j\omega})&=1+2\cos(\omega) \\ &=\sum_{k=-\infty}^\infty \sum_{k=-1}^1 \delta[n-k]e^{-j\omega n} \end{align}$$

Knowing that non-zero only when $n=k$,

$$P(e^{j\omega})=\sum_{k=-1}^1 e^{-j\omega n}$$

Let $M=n+1$,

$$\require{cancel} \begin{align} P(e^{j\omega}) &=\sum_{k=0}^2 e^{-j\omega (M-1)} \\ &=e^{-j\omega (-1)} \sum_{k=0}^2 e^{-j\omega M} \\ &= e^{j\omega} \left(\frac {1-e^{-j\omega 3}} {1-e^{-j\omega}}\right )  \\&=\cancel{e^{j\omega}} \left( \frac {e^{j\omega 3/2}} {e^{-j \omega 3/2}} \right) \left( \frac {e^{j\omega 3/2} - e^{-j\omega 3/2}} {\cancel{e^{j\omega}} e^{j\omega 1/2} - e^{-j\omega 1/2} } \right) \\&= \frac{\sin(\frac{3\omega}{2})}{\sin(\frac{\omega}{2})} \end{align}  $$


How to "draw" the function/wave to send symbols using QAM?

The baseband QAM signal is complex, and the only way to draw it is by doing two drawings, one for the in-phase (real) component, and one for the quadrature (imaginary) component.

The passband QAM signal, though, is real, and it is a pulse-shaped carrier whose amplitude and phase depend on each symbol. Myself, I would "draw" it using Matlab or some other numerical program. I don't see the point in trying to draw it by hand.

For example, in Matlab, assuming QPSK and transmitting three symbols with a rectangular pulse shape, at a rate of one symbol per second and a carrier of 5 Hz:

t = 0:0.01:2.99;
p = ones(1,100);         % pulse shape
s = [1+1i, -1-1i, 1-1i]; % symbols to transmit
u(1:100:201) = s;
qam_bb = conv(u,p);      % pulse shaping to generate baseband signal
subplot(3,1,1);
plot(t,real(qam_bb)); title('Baseband In-Phase'); axis([0 3 -1.3 1.3]);
subplot(3,1,2);
plot(t,imag(qam_bb)); title('Baseband Quadrature'); axis([0 3 -1.3 1.3]);
qpsk = real(qam_bb.*exp(2*pi*5*t*1i));  % upconvert baseband signal
subplot(3,1,3);
plot(t,qpsk); title ('Passband QPSK');

produces this plot:


I don't think of the phase as providing information related to position or location.

A simple example might be illustrative. In 1-D, consider these two signals, x and y:

t = 0:0.001:0.1;
x = cos(2*pi*20*t) + cos(2*pi*30*t);
y = cos(2*pi*20*t + pi/4) + cos(2*pi*30*t - pi/6);

The magnitude spectrums of x and y are the same, but they're not the same signal; if you plot them, you'll see that they're different. In the frequency domain, you need both the magnitude and the phase to fully describe the signals. It doesn't make sense to refer to the harmonics as "being somewhere".


Is there a $z$-transform like for variable sampling rate signals?

Since the powers in the $z$-transform relate to delays or lags, in the philosophy of Laurent series, I would like to mention Puiseux series, a generalization of power series with negative and fractional exponents, or the Hahn (or Hahn–Mal'cev–Neumann) series, with even more flexibility in the  exponents. However, their use in signal processing remains quite shallow to me yet.

I suspect than rational exponents could be still tractable. Indeed, they are used on multiband ($M$-band) filterbanks, where you can spot weirdities like $X(z^{1/2})$ or $X(z^{1/M})$, but most of them are based only on rationals with $M^p$ denominators.
More generic powers (irrational) could be way more involved. The only mentions of Puiseux series I have found on IEEE Xlore are in the domain of control and automation (definitely not my cup of tea).


How to use function spectrum in R to perform spectral analysis

Hi: Your sampling frequency is $f_{s} = 10 $ but the frequency of the third component is $f_{3} = 6$.. Therefore, the nyquist criterion for  the maximum frequency is not being met because $2 * 6 = 12$ which is not less than $f_{s} = 10$. I'm new at this but hopefully this is the reason why. Someone will surely correct me if not.


Write a 440 Hz Sine Wave to .wav-File using Python and SciPy

After some experimenting i came up with the following:

#!/usr/bin/env python3

import numpy as np
from scipy.io import wavfile

sampleRate = 44100
frequency = 440
length = 5

t = np.linspace(0, length, sampleRate * length)  #  Produces a 5 second Audio-File
y = np.sin(frequency * 2 * np.pi * t)  #  Has frequency of 440Hz

wavfile.write('Sine.wav', sampleRate, y)

I pulled out one of my math-books and figured out how to write the sine-function:

Since sin(x) has a period of 2*pi, sin(2*pi*x) has a period of 1/Frequency of 1Hz. Thus, sin(440*2*pi*x) has a Frequency of 440Hz.

The last thing that was unclear to me was how to get the array-length and element-count right!

Opting for a File that has a sample rate of 5 seconds at 44.1 kHz Sample Rate, i need 5 * 44.1k samples, so that was already ok. But since sin(440*2*pi*x) performs 440 periods until t == 1, the array simply has to have a duration of 5 (5 seconds in my case to be precise). enlarging the array boundaries would have squished the array-values together, resulting in a higher pitched sound.

Sorry to anyone that has already scratched their head, i probably should have thought for a little bit longer!


Finding number of independent samples using autocorrelation

The autocorrelation you obtained shows that the signal becomes decorrelated with itself for specific time delays, which are the zero-crossings of the autocorrelation. The autocorrelation is not zero except at these time instants.

You can also define a threshold $\tau$ and say that, if the autocorrelation is less than $\tau$, then for your purposes the signal is decorrelated for any time shift $\Delta$ for which the autocorrelation is less than $\tau$.

When calculating the autocorrelation in Matlab (or similar), the lag is calculated in number of samples. For example, in your case it seems like the autocorrelation is zero for a lag of 200. Since your sampling interval is $T_s = 1/5000 = 200\,\mu$s, then your signal becomes decorrelated for a shift of 200 samples, or a time shift $\Delta = 200 T_s = 40\,$ms.


Compare feature vectors of an image in MATLAB

If you just need OCR to work, then please try the ocr function in the Computer Vision System Toolbox.

If you need to write your own, take a look at this example, showing how to train an SVM classifier using HOG features for recognizing hand-written digits.


Fourier transform before calculating the phase

The discrete Fourier transform (DFT) does by definition transform your image into the frequency domain, mapping it to a matrix of complex numbers.

Again, by definition, a complex number has a phase and a magnitude component.

Thus, you don't "need to perform FT before calculating the phase of an image", you get the phase of the image as a result of the FT.

It's a bit like saying "why do I need to perform an addition before calculating the sum of 2 number?".

I would suggest you to learn more about Fourier Transform


There doesn't seem to be much discussion of CRC in this forum but on the sister site math.SE, there are several postings. See, for example, CRC Computation or What is the theory behind seeding a CRC.


How to model state space for complex valued system correctly in SIMULINK (MATLAB)?

You can split up the real and imaginary part of the state into their own seperate states. Namely by defining $x_r=\mathrm{Re}(x)$, $x_i=\mathrm{Im}(x)$, $A_r=\mathrm{Re}(A)$, $A_i=\mathrm{Im}(A)$, $u_r=\mathrm{Re}(u)$, $u_i=\mathrm{Im}(u)$, $B_r=\mathrm{Re}(B)$ and $B_i=\mathrm{Im}(B)$ then the differential equation can also be written as

$$
\dot{x}_r+i\,\dot{x}_i=(A_r+i\,A_i)(x_r+i\,x_i)+(B_r+i\,B_i)(u_r+i\,u_i)
$$

which when split into their real and imaginary part gives

$$
\begin{bmatrix}
\dot{x}_r \\ \dot{x}_i
\end{bmatrix}=
\begin{bmatrix}
A_r & -A_i \\
A_i & A_r
\end{bmatrix}
\begin{bmatrix}
x_r \\ x_i
\end{bmatrix}+
\begin{bmatrix}
B_r & -B_i \\
B_i & B_r
\end{bmatrix}
\begin{bmatrix}
u_r \\ u_i
\end{bmatrix}.
$$

The output can be expressed using

$$
y=C\begin{bmatrix}
I & i\,I
\end{bmatrix}
\begin{bmatrix}
x_r \\ x_i
\end{bmatrix} + D\begin{bmatrix}
I & i\,I
\end{bmatrix}
\begin{bmatrix}
u_r \\ u_i
\end{bmatrix}.
$$


As illustrated by the simplified figure and the Matlab code below, such a graph can be obtained with a standard discrete wavelet decomposition, possibly with a periodic extension to preserve the size, over 3 to 5 levels of decomposition. Wavelet coefficients are simply concatenated in the standard low-frequency  /high-frequency  fashion, with scaling coefficients first, and followed by the wavelet coefficients.

load ecgData;
data = x(1:256,1);
waveletName = 'db4';
waveletLevel = 4;
dwtmode('ppd');
[C,L] = wavedec(data,waveletLevel,waveletName);
subplot(2,1,1);  plot(data);axis tight
subplot(2,1,2);  stem(C);axis tight


Set $\Delta=\frac{1}{99}\log_{10}\left(\frac{3000}{0.01}\right)$

Then the $n_{th}$ knob position should be

$p_n=10^{n\Delta-2}$

Where $n$ runs between 0 and 99 in integer step.


Suppose you have a signal $x[n]$, with $n \in {0, 1, ... N - 1} $. The same-size DFT is defined by:

$$
X[k] = \sum_{n = 0}^{N-1} x[n] e^{-j\frac{\Large{2 \ \pi n k}}{N}}
$$

The frequency resolution is going to be how many Hz each DFT bin represents. This is, as you have noted, given by $\frac{f_s}{N}$.

If on the other hand you had zero-padded your signal, such that $N_{zp}$ is greater than  $N$, then a more apt term of frequency granularity is given by $\frac{f_s}{N_{zp}}$

  Asking this because the spectrum is symmetric for real-valued inputs.

That is irrelevant. The frequency resolutions/granularities are given by the above.

  o, say I have fs=1000 Hz and N = 1024, where N is the number of FFT
  points. Now, is the frequency resolution 1000 Hz1024=0.9766 Hz or 1000
  Hz0.5∗1024=1.9531 Hz?

If your sampling frequency $f_s = 1000$ Hz, and you are taking an $N=1024$ (same-sized) FFT, then your frequency resolution is $\frac{1000}{1024}$, which is equal to 0.9766 Hz/bin. If your $N_{zp} = 1024$ (FFT length after zero-padding), then your frequency granularity is 0.9766 Hz/bin.