Spectral window when FFT is not a power of 2

500 is an exact integer multiple of the period of the cosine function used in your raised cosine window.  Therefore there are zeros in the spectrum at all multiples of the cosine's frequency (due to those basis vectors being exactly orthogonal).  512 is not an integer multiple of that cosine's period, so none of the correlations against the DFT basis vectors of a 512 length FFT will be zero (just small).  Thus the difference between the two spectrums.


When is a LTI IIR $H(z)$ system minimum phase

but what does the stability mean when talking about a system?

It means that all poles are INSIDE the unit circle.

and if so how does it turn out that when a system is causal and stable its also min phase.

Sorry, you got this wrong. Causal, stable and LTI does NOT imply minimum phase. A simple counter examples is a one-sample delay. It's causal, stable and LTI but definitely not minimum phase.
Minimum phase requires that all zeros are also INSIDE the unit circle as well.
Update
When you invert a system the poles become the zeros and the zeros become the poles. A minimum phase has both all zeros and poles inside the unit circle, that means that poles and zeros of the inverse will also be inside the unit circle and hence it's also minimum phase.
A non-minimum phase system has zeros outside the unit circle. If you invert the system, this zero turns into a pole that's outside the unit circle and hence the inverse is unstable (or non-causal).


How could one can say that frequency spectrum having Nyquist frequency 'x' is the most appropriate representation of a given signal?

The most important thing is to sample at more than twice your signal's bandwidth; the period is largely irrelevant. Say your signal with period 2 is a square wave -- its bandwidth is infinite and you'll never be able to sample it correctly (at least in theory). If your signal is a sine wave, then a period of 2 means its bandwidth is 0.5 and you'll easily sample it adequately.

Personally, whenever possible, I'd sample:

a bit above Nyquist (say, at $2.2B$, where $B$ is the signal bandwidth), to maximize the number of useful frequency bins, and
for as long as possible, to increase the spectral resolution.


What pattern maximises a JPEG file's size?

Pursuing on @OlliNiemitalo, JPEG allows some loss related to frequency quantization. So any file can be compressed to any size, hence you cannot maximize a JPEG size without setting some constraints, like quality. Let us assume a standard metric, like mean-squared error.

Now, let us check some of the JPEG features for image sparsification, and let us try to go the opposite direction:

JPEG assumes similarities in RGB color planes, correlated with a luminance/chrominance transform for instance. So you choose independent R, G, and B, you are good
8$\times 8$ patches are compressed almost independently, except for the average, predicted for the blocks before (raster order). If you have neighboring blocks with very different means, good
the more the high-frequency, the less a 8$\times$8 block as natural frequency, the better. They avoid long trails of zeros encoded by an EOB symbol.

One cannot easily fulfill all these properties together,  especially with RGB planes made of 8-bit integers.
Forgetting about the average-thing, generating  independent RGB images, with each pixel randomly $0$ or $255$ (ie high range). Quantization is likely to produced here a lot of false colors. This could give you very bad images for the JPEG in MSE. However, original and compressed ones could look as ugly.

You can find other fractal-based or fancy-patterned proposal in
What is the least JPG-compressible pattern? (camera shooting piece of cloth, scale/angle/lighting may vary)


References ISAR - Inverse Synthetic Aperture Radar

A couple of books I'd recommend that are similar to the Richard's book

Stimson's "Introduction to Airborne Radar" has a chapter on ISAR.
For a bit more detail you can look at
Wehner's "High Resolution Radar"

If you look at the Polar Format algorithm in SAR (I'd recommed Soumekh's - Synthetic Aperture Radar Signal Processing) it is very similar to the ISAR algorithms after you corrected the target's translational motion (this is not necessarily trivial). Soumekh's book doesn't really address the ideas of motion compensation, target tracking, and range cell migration - these are essential concepts if you are implementing an ISAR algorithm on Airborne data. If you are just trying to get an idea of ISAR processing then you just need a little bit about them.

A good book that covers a wide range of ISAR topics is:"Inverse Synthetic Aperture Radar Imaging: Principles, Algorithms and Applications" by Chen and Martorella.

One piece of advice, if you are reading papers - pay very close attention to the type of transmitted signal and the type of range compression that is used. Often authors just state the signal after "preliminary" processing is given by ... . If you do deramp on receive vs traditional matched filter, the nature of the signal is very different. Similarly, if a stepped frequency waveform is used then the resulting signal processing is different. This is one of the reason's I liked Soumekh's book.

One other book is:Inverse Synthetic Aperture Radar Imaging with MATLAB Algorithms by Caner Ozdemir. This book explains the concepts well, but I was unable to follow the mathematical development and the Matlab code is terrible. You can't easily modify it change the processing parameters.


Hough forests for object recognition are initially proposed by Jurgen Gal. The paper was named "Class Specific Hough Forests for Object Recognition".

Here is his homepage. Here are the source code and the details. I think you could make use of those.

Certainly better methods exist and one could well talk about them. However I am sticking to the Hough framework, as you mentioned that in the first place.


A sum of two sine waves can be expressed equivalently as the product of two sine waves.  This is a phenomenon commonly referred to in musical circles as a beat frequency, which can be used to tune instruments by ear.  For your example, your input can be expressed as a 103.5kHz frequency multiplied by a 1.5Hz frequency.

I wouldn’t use the term “squaring” for what you are doing as much as “clipping the ever living hell out of”.  This would probably be expressed best using concepts from intermodulation distortion.  However, in your application, the square nature of the output of the waveform can be analyzed more simply anecdotally.  I would expect to see a square wave at the 103.5kHz with glitches at the 1.5Hz zero crossings.  But notably, it would look a lot like a 103.5kHz waveform, with some side band nonsense caused by glitches.

I do not believe that the flip-flop can be expressed using any common signal processing techniques, but again it can be analyzed anecdotally.  As you’ve hypothesized, I would expect the time between edges to double, because you’re skipping every other edge.  However, this is a grosse generalization.  It would be exactly true for a sinusoid, less true for the input you’ve proposed, and untrue for a generic input.  So for you’re example, mostly energy at 51.75kHz.

At this point, you still have a square wave, roughly at 51.75kHz, with some side band garbage cause by the glitches.  One could apply an arbititary amount of filtering to extract the sine wave at this frequency, resulting in the output waveform you described.

Something else to keep in mind.  A frequency transform has practical limitations.  A proper transform requires all time domain points of the input signal.  Using less than infinite samples is equivalent to windowing, which causes side lobe patterns about dominant frequencies, which is a topic unto itself.  Additionally, a signal sampled at discrete time intervals (like on a computer) has a limited bandwidth.  When you “clip the ever loving hell” out of a signal, you are widening the bandwidth of the input signal significantly.  If care isn’t taken to ensure that the discrete time sample rate is high enough, a significant amount of aliasing distortion will occur, hindering proper analysis.

Sorry for being a bit of a rant, but your question touches on a lot of points.  Cheers.


What does the 'B' parameter represent in Wikipedia's page on window functions?

See paragraph from wiki. There is an explanation of B metric:

  Each figure label includes the corresponding noise equivalent bandwidth metric (B),[note 1] in units of DFT bins.

Note 1:

  Mathematically, the noise equivalent bandwidth of transfer function H is the bandwidth of an ideal rectangular filter with the same peak gain as H that would pass the same power with white noise input. In the units of frequency f (e.g. hertz), it is given by:
  $$
B_{noise}={\frac {1}{|H(f)|_{max}^{2}}}\int _{0}^{\infty }|H(f)|^{2}df.
$$


Resources on Solving Convex Optimization Problems in the Compress Sensing Field

There are few options:

Stephen Boyd, Lieven Vandenberghe - Convex Optimization.
This is the classic in this field. Very well written book.
Also have a look on other papers of Boyd on similar subjects such as the The Alternating Direction Method of Multipliers (ADMM).
They also have a great MOOC Course Stanford Online CVX 101 - Convex Optimization.
Amir Beck - Introduction to Nonlinear Optimization - Theory, Algorithms and Applications & Amir Beck - First-Order Methods in Optimization.
Amir is a great teacher and his book are both deep, practical and easy to read.
CMU Statistics 36-725 - Convex Optimization: Fall 2016 by Ryan Tibshirani.
Videos are available on YouTube.

I think once you skim through those you'll be able to handle most of the cases in Convex Optimization.


Evaluate frequency reponse of cascaded filter with bitwidth truncation in matlab

As bit truncation is not a linear operation, there's strictly speaking no frequency response, as the overall system isn't linear.
But: assuming somewhat benign input, I'd simply model the quantization noise as uncorrelated (uff, I know that's a big jump!)
linear system H1 -> + quantization noise -> linear system H2 -> + quantization noise
or equivalently
linear system H1 -> linear system H2 ->  + colored quantization noise -> + quantization noise
allowing you to combine H1 and H2 and then just add non-white noise of the appropriate variance.
Note that the variances of your two noise sources depend on the expected quantization error, which in itself depend on the bit-widening properties of H1, H2. You could math out and find that variance, or you could run these filters with enough test signal and simply measure the squared error between 12-bit truncated and untruncated output.
I'm not a great fixed-point filter designer myself, but I'm pretty convinced there's plenty of material out there helping you figure out optimal truncation.


If you have some numerical  methods background,  Richard Hammings, Digital Filters, is a good place to start.  The book is from 1977 but the math hasn’t changed.  The book has also been released by Dover so it has a very low price, unlike Oppenheim and Schaefer’s book and also unlike Oppenheim and Schaefer, is compact.
The math level isn’t advanced.

Hamming is one of the founders of Information Theory.  He essentially invented parity codes. He was also the first president of SIAM.


Band limited signal: Infinite time duration, time limited signal : infinite bandwidth

Here is an argument that doesn't require the uncertainty principle.

Assume $y(t)$ is a signal that is not identical to zero, with Fourier transform $Y(f)$. Let $s(t) = \text{rect}(t/T)$ a rectangular pulse of duration $T$ and Fourier transform $S(f) = |T| \text{sinc}(T f)$.

Then, $x(t) = y(t)s(t)$ is a time-limited signal with Fourier transform $X(f) = Y(f) \ast S(f)$. Since $S(f)$ has infinite support, so does $X(f)$.

Since any time-limited signal $x(t)$ can be defined in this way for a suitable choice of $y(t)$ and $s(t)$ (and a possible time-shift that doesn't affect the result), then we can say that all such signals have infinite bandwidth.

The same argument can be applied for the inverse Fourier transform: band-limiting a signal involves multiplication by a rectangular pulse in the frequency domain, or a convolution with a sinc in the time domain, resulting in an infinite-duration signal.


Is frequency resolution ("bin size" inverse) always directly proportional to sampling duration (window length)?

For orthogonality of basis vectors, the DFT (or FFT) basis vectors all have a different but exact integer number of full  cycles in the DFT length, e.g. if the FFT result bin for one basis vector has K full cycles, the next bin up will be for K+1 full cycles.

So if you want a DFT to separate 2 frequencies of sinusoids, the closer they are in frequency, the longer a time interval it takes for the two waveforms to become separated by exactly 1 full period or cycle, the longer a DFT (of FFT, or Goertzel, or wavelet, etc.) it takes to fit that time interval for measurement, and have at least a full result bin of separation.

Even if you don’t want to wait until two frequencies of sinusoids differ by one full cycle, say there is so little noise that you can tell apart 2 the sinusoids if they differ by only 0.1 degree (pi/1800.0) of phase change within some window, the rule still applies.  The closer the 2 frequencies are that you want to tell apart or resolve, the longer you will have to wait for them to differ by that 0.1 degree of phase shift or change that you can barely detect (regression fit within some statistical error bound or distribution, etc.)

Otherwise you can’t tell them apart, and thus claim a resolution or level of accuracy that can.

And the more noise, the bigger the difference you will have to wait for (maybe more than 1 degree, maybe even up to more than 2 full periods of difference within the measurement time window).

Note that even with analog filters (or IIRs), generally the narrower the bandwidth, the higher the Q, the longer the impulse response, thus the longer the measurement delay to respond to or detect that narrower frequency range.

So the “tightest” fit requires either more time, and/or low noise (including interference and measurement error), or a very high signal-to-noise ratio (or of the assumption of such, and hope that that assumption is actually correct.)


Bounds of the difference of a bounded band-limited function

As shown in the answers to this question, for continuous-time signals, the bound predicted by Bernstein's inequality is achieved with equality by a sinusoidal signal with a frequency equal to the upper frequency limit.

In analogy to this, in this answer I'll show a bound on $\big|x[n]-x[n-1]\big|$ for a discrete-time sinusoidal signal $x[n]$ at angular frequency $W$:

$$x[n]=A\sin(nW+\phi)\tag{1}$$

For the signal $x[n]$ given by $(1)$, the largest value of $\big|x[n]-x[n-1]\big|$ for all $n\in\mathbb{Z}$ occurs for two values $x[k]$ and $x[k-1]$ symmetrical to the point of the largest derivative of a (continuous) sinusoid, i.e., for $x[k]=A\sin(W/2)$ and $x[k-1]=A\sin(-W/2)$. Consequently, the bound for a sinusoid with frequency $W$ is given by

$$\big|x[n]-x[n-1]\big|\le 2A\sin\left(\frac{W}{2}\right)\tag{2}$$

This bound is achieved with equality for the signals

$$x[n]=A\sin\left(nW+\frac{(2l+1)W}{2}\right),\qquad l\in\mathbb{Z}\tag{3}$$

Note that the constant $A$ does not generally equal the maximum amplitude $B=\max |x[n]|$ of $x[n]$. Depending on the sampling phase, the maximum amplitude $B$ lies in the following interval:

$$A\cos\left(\frac{W}{2}\right)\le B\le A\tag{4}$$

Consequently, we have

$$A\le\frac{B}{\cos\left(\frac{W}{2}\right)},\qquad 0<W<\pi\tag{5}$$

Combining $(2)$ and $(5)$ we get

$$\big|x[n]-x[n-1]\big|\le 2B\tan\left(\frac{W}{2}\right),\qquad B=\max_n\big|x[n]\big|\tag{6}$$

Note, however, that $(6)$ is only useful for $\tan(W/2)<1$, i.e., for $W<\pi/2$ because for any $x[n]$ with $\max|x[n]|=B$ we must have

$$\big|x[n]-x[n-1]\big|\le 2B\tag{7}$$

I believe that the bound $(6)$ holds for all discrete-time band-limited signals with a maximum frequency $W$ (i.e., with $X(e^{j\omega})=0$ for $|\omega|\in (W,\pi]$), but I don't know how to show it.


Why does toning up a piano make it sound similar to that a string instrument

The human ear, in collaboration with the brain, recognises what it is that is being perceived (as well as where it is coming from) within the first few ms of the sound waveform.
This "what is perceived" (is it a piano? is it an organ? is it a poggo-stick?) is called Timbre.
This modification of the attack phase is sometimes used deliberately.
Here, it is not the delay that is effecting the modification. Notice his hand on the volume knob when the guitar now starts sounding like a violin. It is clearer from 01:43 onwards in the "Without echo" part of the video. Here, he is killing the guitar's distorted "twang" as the string is excited and "brings it in" just a little after that so that you get more of the tone and less of the "twang". It then sounds like a violin. You might also notice on a piano, if you simply play the waveform after the attack, it sounds more like a decaying pipe-organ rather than a piano.
By "speeding up" the sound, you are distorting those crucial first few ms that the ear/brain are using to recognise the instrument. For example, the instrument now has a sharper attack than a piano. The reason it begins sounding "stringy" is another byproduct of the speeding up as the combination of sinusoids (and their relative phase) might now be tending towards a sawtooth wave that is typical of stringed instruments.
If you are looking to pitch shift but keep the duration of the sample the same, you need to be looking at techniques that attempt to preserve both the amplitude and phase spectrum of a sound.
Librosa implements this with pitch_shift() but other implementations are also available.
Hope this helps.


Why does Matlab spectrogram look stripy when windowed input is zero-padded?

Since I don't have your original data, I'm going to give the answer from my intuition. Try Hann window function and 50% overlap and see what happens.
I have tried a human voice signal as input and whether zero-padding is applied doesn't change the results.
% signal = ...
winsize = 512;
win = hann(winsize, 'periodic'); % use hann window instead of rectangular window
shift = winsize/2; % 50% overlap
nfft = winsize;
% nfft=round(1.5*winsize);
srate = 48e3;
figure;
spectrogram(signal, win, winsize-shift, nfft, srate, 'yaxis');

Here are the results:

Forgive me for the Chinese xlabel and ylabel, they are respectively "Time (s)" and "Frequency (kHz)", and the colorbar label stands for "Power/Frequency (dB/Hz)". I don't see any visible difference between the two figures.
In addition, you can use a smaller FFT size. Generally, using such a long FFT size means that you need an extremely high frequency resolution.
$$
\Delta f = \frac{f_s}{\mathrm{NFFT}}
$$
where $\Delta f$ is the frequency resolution, which is 0.0168 Hz when $f_s$ = 78 Hz and NFFT = 4639. While a super high frequency resolution means a super low time resulution with the fact that $\Delta T = 1/ \Delta f$. In your case you only have 10 time sample results.

For the second question, 0 Hz means the direct current (DC) component. Maybe you should add some information about what kind of signal are you analysing.


There are many ways to estimate the instant of glottal closure (GCI).  Here (PDF link) is a recent journal article that reviews many of them.  The freely-available VOICEBOX toolbox for Matlab includes a couple GCI estimation techniques.

For most approaches, the fundamental idea is the same: the speech signal is can be viewed as an excitation signal resonating an all-pole filter, and the location of the excitation instants is the instant of glottal closure.

A simple GCI estimation algorithm might work like this.  The all-pole filter (encapsulating the contributions of the vocal tract and glottal source) can be estimated using linear prediction.  Inverse filtering the speech signal with this (time-varying) filter will give the linear predictive residual signal, which is an estimate of the excitation signal.  For normal voices, this signal will look like a noise impulse train.  Peak-picking this signal will give good results for these voices during sustained vowels.  Voice offsets and onsets will be more troublesome, as well as voices that are very breathy or irregular.

(A quick note about your figure: the signal marked as the residual seems to include glottal source signal, and so it not a linear prediction residual in the same sense as what I refer to in the above.)


The networks generally learn to discriminate without the information of a generative model. This is why, I would say, they tend to learn not to be invariant. This is also the reason why they perform very well at fine-grained classification tasks. Their level of invariance is largely influenced by the data, which is provided. So if the data is augmented so as to cover all rotations, then the network starts to be invariant against rotations. This is also similar for other types of non-trivial transformations. Such a property is desirable, as we do not want to be invariant against the transformations that are not shown to us.

Having said that, they do, intrinsically develop certain invariance. This is mostly due to the pooling mechanism. So if you were to design your pooling method (corresponding to certain norm), some more invariance can be achieved.

Also, scattering transform of Mallat construct translation-invariant operators on $L^2(R^d)$, which are Lipschitz-continuous to the action of diffeomorphisms. They can be made to be theoretically invariant to actions of compact
Lie groups.

Please check Goodfellow et. al. for the degree of invariance, the networks could develop.


How to determine, calculate and plot a probability distribution for a given set of numbers?

A kernel density estimation based on your data:

Depend on your programming language, you may find the corresponding fitting functions with different distributions.

EDIT

I tried to use the allfitdist function with statistic toolbox in Matlab in your data fitting.

x = [24.988 25.224 25.212 25.066 25.310 24.963 24.826 24.944 25.490 25.176 24.740 24.988   ...
24.994 24.722 24.510 24.863 25.249 24.931 25.218 24.907 25.459 25.292 25.163 25.447 ...
25.072 25.360 24.976 25.103 25.237 24.845 25.237 25.206 25.109 24.969 24.976 25.194 ...
25.218 24.975 24.777 24.832 24.564 24.637 24.919 24.931 24.944 24.994 25.378 25.042...
25.182 25.151 25.097 25.103 25.030 25.018 24.783 24.722 24.667 24.963    ];

[D PD] = allfitdist(x,'PDF');

Unfortunately the fitted probability density is not good enough:

The curves for each possible distribution are all overlapped. Like I mentioned yesterday, you may need more sample number to obtain a better estimation.


When you express the cosine as $\cos(2\pi f_kt)=\frac{1}{2}(\exp(j2\pi f_c t)+\exp(-j2\pi f_ct)$ and insert it into the convolution integral you get:

$$
h(t)*\cos(2\pi f_ct) = \frac{1}{2}(\exp(j2\pi f_ct)H(f_c)+\exp(-j2\pi f_ct)H(-f_c))
$$

In case your $h(t)$ is real-valued, you have $H(f_c)=H(-f_c)^*=|H(f_c)|\exp(j\phi)$, and you get
$$h(t)*\cos(2\pi f_ct)=\cos(2\pi f_ct+\phi)|H(f_c)|.$$

In case you have a complex impulse response, this does not hold.


Finding the impulse response of a system

One can see that the given expression can be decomposed as
$$
\begin{align}
H(j\omega) &= \frac{1+0.5 e^{-j\omega}}{1-1.8 \cos(\frac{\pi}{16}) e^{-j\omega}+0.81 e^{-j2\omega}} \\
&= \frac{A }{1-0.9 e^{j\frac{\pi}{16}} e^{-j\omega}} + \frac{B}{1-0.9 e^{-j\frac{\pi}{16}} e^{-j\omega}}
\end{align}$$

where $A = B^*  = 0.5 - j 3.9375 = 8/16 - j 63/16 $.

Then the impulse respons will be:

$$h[n] = A (0.9 e^{j\frac{\pi}{16}})^n u[n] + A^* (0.9 e^{-j\frac{\pi}{16}})^n u[n]$$

which can be simplified as:

$$h[n] = 2 \cdot 0.9^n |A| \cos( \frac{\pi}{16}n + \angle{A}) u[n]$$

where $\angle{A}$ is the phase angle of $A$. Following is the resulting sequence plotted, from $n=0$ to $n=35$, using MATLAB/OCTAVE.


Step response of a differentiating system

If your system is an ideal differentiator with input-output relation

$$y(t)=\frac{dx(t)}{dt}\tag{1}$$

then its transfer function is

$$H(s)=\frac{Y(s)}{X(s)}=s\tag{2}$$

From the equation in your question you obtain for its step response

$$A(s)=1\tag{3}$$

which in the time domain corresponds to a Dirac delta impulse:

$$a(t)=\delta(t)\tag{4}$$

This is also intuitively clear, because a step input has a derivative of zero everywhere except at $t=0$, where it is discontinuous.


Gaussian weighted distance between pixel amplitudes: motivations and sources?

Maybe not an exact answer, but I'll to give a direction.

What you are using is essentially an RBF-Kernel. First, it has a ready interpretation as a similarity measure (satisfies Mercer's conditions). In essence, other such kernels can be used. For instance, an inner product (linear kernel) would define an angular distance between pixels/features. For a deeper understanding, I would dive in the literature of kernels.

Having said that, the simplest form of such relations are called inverse distance weighting (Grady could have used this as well), where any identity like $1/d(g_i, g_j)^p$ can be used. These methods also cover a wide spectrum.


What are the values on the x-axis when I plot the impulse response of a lowpass filter?

For digital filters the x-axis is samples.


Autocorrelation of a stochastic process which is a sum

Your proof is incorrect.  the step indicated below is not true, i.e.:
$$
E[X_n X_m] \not = E[X_n] E[X_m]
$$
in general.


How to interpret these filter coefficients?

Sorry, but that does not make any sense

A 3rd order BW lowpass can be described most simply with a single number: the cutoff off frequency
It's the not numerator a denominator since this would have 4 numbers and the numerator would have the shape of x*[1 3 3 1] and the denominator should start with 4.
It's not the poles and zeros since all the zeros are at -1 and one pole is conjugate complex
It sure as heck isn't state space either
The transfer function is a boiler plate text book copy and also mostly wrong for a third order BW (which has one biqaud and one single real pole section).

You are not missing anything obvious and having you guessing blindly doesn't serve any useful purpose. You have to get back to the client and tell them that you need proper specification or you can't get the work done.


Check whether a given equation is linear


Simple and efficient algorithm to detect frequency and phase of a sine signal

Note: I originally posted this answer for the Stack Overflow copy of this question, before realizing that it had also been asked here.  It somewhat duplicates pichenettes' answer, but I felt it still worth (re)posting here, since it includes some extra details.  (Whether those details are useful or not, I'll leave for you and the OP to judge.)

If you know your signal is a pure sine wave, you can just use zero crossing detection.  Each cycle of the sine wave will have two zero crossings: one from negative to positive, and one from positive to negative.  This will be a lot simpler and more efficient than trying to do something fancy like Fourier transforms.

There are a few details to keep in mind, though:

It's OK for the signal to be slightly biased, but if the bias might exceed the amplitude of the signal, you'll need to correct it somehow.

You can either do this before sampling with an analog high-pass filter, or you can track a moving average of your sampled signal and use it as the "zero level" to compare to.  Or you can, instead, look at zero crossings of the difference between successive samples (which correspond to the maximum and minimum of each cycle in the signal) to avoid any bias issues.
If your input (or ADC) is noisy, your samples might randomly fluctuate around the zero level when the signal is close to it.  In such cases, naïvely comparing successive samples might detect multiple zero crossings where there's really just one.

One way to fix this issue is to smooth your signal before processing it, but it may be easier and more efficient to implement hysteresis in your zero-crossing detector.  That is, you'd only detect a positive-to-negative crossing when the signal dips below some pre-set threshold level −ε, and a negative-to-positive crossing only when it rises above +ε, like this pseudo-C code:

boolean isPositive = (firstSample > 0);
while (running) {
    int signalLevel = getNextSample();
    if (isPositive && signalLevel < -threshold) {
        isPositive = false;
        handleZeroCrossing();
    }
    else if (!isPositive && signalLevel > +threshold) {
        isPositive = true;
        handleZeroCrossing();
    }
}

In fact, you may not even need an algorithm for this at all, since, as noted in the comments below, you can implement zero-crossing detection with hysteresis in hardware with a simple Schmitt trigger, which basically converts the sine wave input into a square wave signal with the same frequency and (almost) the same phase, which you can then read as a simple digital input.  You might even be able to use the output of the Schmitt trigger to drive the MCU interrupt pin directly.


An ideal pilot signal for channel estimation (without any apriori further knowledge of the channel) would transmit power uniformly across the channel and be easy to replicate in the receiver for channel estimation, as channel estimation cannot be determined where no signal exists for comparison, as well as carrier and timing recovery for a synchronous digital modulation data link. PN sequences are good choices due to their characteristics that mimic white noise. Further and due to the same properties resulting in an auto-correlation that appears as a single impulse as expected with white noise, they are also useful for synchronization such as estimating the start of a data packet. The PN sequence also provides fast convergence for timing and carrier recovery in a digital data link, but you may also consider adding a certain duration of a 10101010 sequence after the PN sequence (if the overhead allows) which will provide further robustness for timing recovery.
A chirp could also be used providing both properties of having a white spectrum (on average across the frequency band of interest for which the pulse is applied) and a sharp autocorrelation function (as demonstrated by the matched filter output being a “compressed pulse) although the use of PN codes could more easily offer a greater set of uncorrelated patterns that could be used in the same pilot time interval which could be used for other features related to packet identification (if desired). A chirp does not readily provide for carrier recovery since it is difficult to distinguish carrier offset from time offset, nor does it provide information on the data rate that the toggling of the PN sequence would provide.
The primary use and purpose of a pulse shaping filter is for applications where the spectrum must be contained. (With no pulse shaping filter, and the transmission of rectangular pulses, the ideal matched filter becomes and integrate and dump filter (or equivalent moving average process) and there is no penalty in doing this other than the much wider transmit spectrum. Distortion to pulse shaping is a non-linear distortion mechanism such as  driving amplifiers too hard into compression rather than linear distortions such as pass-band ripple, reflections etc that can be corrected for with channel equalization (as these distortions end up being similar to multipath that the channel introduces).   The resulting spectrum after a pulse shaping filter is also flatter across a narrower spectrum (just over B wide for a symbol rate of B symbols/sec) while a rectangular pulse will have about a 4 dB roll-off over that same bandwidth. For this reason the use of a pulse shaping filter for channel estimation patterns will result in more efficient channel estimation applying more equal power across the channel. When timing is of primary concern (such as with GPS), pulse shaping will degrade time resolution since for that the widest bandwidth is desired (for the sharpest autocorrelation peak). To create a matched filter, it is typical to decompose the pulse shaping into the cascade of two identical filters (thus square root filters), with one in the transmitter and the other in the receiver.
Here is a link that describes the trade space with pulse shaping further, but again first consider what your spectrum containment needs to be and then design from that; the tighter the pulse shaping the more complicated the receiver will be both with the duration of the matched filter and in the timing recovery algorithm/loop.
Pulse shaping with RRC : Number of taps


Separate key-related frequencies from audio signal

Attenuating the harmonics of your instruments isn't really a viable option, because you'll end up annihilating the signal itself. The harmonics here are what are sampling the behaviour of your instrument.

I'm not sure if I'd call it simple, but a possible solution to separating the contribution of pitch in such a classification system might be cepstral liftering.

This is, calculate a cepstral representation of the recording of each instrument (the Fourier transform of the logarithm of the frequency response) and separate the periodic components (pitch) from the aperiodic (perhaps timbre?). Such separation is known as liftering - the analog of filtering in the cepstral domain.

The separated signals should roughly correspond to a periodic signal and an impulse response better characterising the instrument.

This is an approach used in to separate the harmonics of speech (carrying pitch information) from the transfer function of the vocal tract (which generates a filter forming vowels amongst other characteristics).

I'm not sure if this will work in your case, but I think it would be worth exploring.


Trying to apply Parks-McClellan optimal FIR filter design

There are a few problems with your designs. First, the desired bandwidths and transition bands (between passbands and stopbands) are very narrow. Second, some of the center frequencies are very low, especially the ones at 20/100000 and 300/100000. And finally, your chosen filter length is very large. I do not know the implementation of the Remez exchange algorithm that you use, but I've seen several implementations that run into serious numerical problems for such filter lengths. You can see in your designs that around Nyquist there are some small peaks in the filter response. These are due to numerical problems.

So you should try to make the passband and the transition bands a bit wider and reduce the filter length to see if you can get a reasonable design. From there you could try to make the specifications a bit more stringent and increase the filter order and see how far you can go.

In Octave, the following specification gives a reasonable result (center frequency: 0.1, filter length: 300, '1' corresponds to Nyquist):

h=remez(300,[0,.07,.085,.115,.13,1],[0,0,1,1,0,0]);


In pulse code modulation what is pulse word represent?

Think of a DAC as a thing that is emitting a train of pulses, each new one hard up against the last one, and each one exactly as wide as a sampling interval (sampling interval = $T_s = \frac{1}{f_s}$).
The code in PCM is the height of each pulse.
And if the name seems odd now, just remember that it made perfect sense in the 1925 to 1948 time period when the basis for this stuff was being laid down.


Why is the Hamming and Hanning window the most widely used for speech spectograms?

any window (almost any) is better than a boxcar window, unless your spectrum is flat.  speech signals are not flat.

I suspect that Hamming and Hann are commonly used because you can avoid scalloping loss for 50% overlap. One gets a sufficiently favorable tradeoff between time resolution and frequency domain dynamic (side lobe level) range. There is also the fact that these windows are cited in papers that go back to the 70s.

You can also apply the widows with a 3 point convolution (3 coefficients instead of N) in the frequency domain, which might have had some advantages in hardware in the past when memory was a constraint.

Classical fixed windows are chosen according to characteristics like sidelobe level. They aren’t like adaptive windows that use the data itself.  It’s relatively pedestrian processing.

in summary, boxcar windows are for flat spectra.


Regarding the nature of white noise in images

A signal is just a varying quantity over an axis.

Think of a greyscale image.

Think of a single row of that image. Can you see how the intensity evolves over the distance on the horizontal axis?

In an image, you get two axes (x and y), whereas in sound, you just get one.

In (mono) sound, your instantaneous intensity is the only quantity varying. In an image, that quantity might be brightness alone, or brightness, hue and saturation, or red, blue and green intensities, or whatever channels your color model defines.

Anyway, the point about white noise is it's autocorrelation function; and that's easy to define on the two axes of an image, exactly the same as for sound:

Rows of an image: If you shift the row(s) by any amount (other than 0), then the expected sum-of-products between the original and shifted row is 0.
Columns of an image: If you shift the columns(s) by any amount (other than 0), then the expected sum-of-products between the original and shifted column is 0.

If both properties are given (and in fact, the correlation in each direction, even if you "sliced" your image arbitrarily in diagonals), then your 2D noise is white.


Is there a way to double all frequencies?

That's typically referred to as "pitch shifting" and it's typically a combination of re-sampling and time stretch/compression. Re-sampling is straight forward but time stretch/compression is a bit more tricky and the best methods depends on the application requirements, the characteristics of the signal and the amount of stretch/compression you'll need. See http://blogs.zynaptiq.com/bernsee/time-pitch-overview/ and https://en.wikipedia.org/wiki/Audio_time-scale/pitch_modification to get started


You can alternatively create a DFT matrix in matlab using this code:

exp(-1j*2*pi* ((0:N-1)/N).' * (0:N-1))

And the IDFT matrix thus:

1/N * exp(1j*2*pi* ((0:N-1)/N).' * (0:N-1))

As the only difference betweenm DFT and IDFT is the sign and a scaling factor.

You could alternatively just do:

ifft(eye(N))

But this doesn't get around needing the full DFT matrix, and honestly I don't see an obvious way of perfoming a DFT by matrix multiplication without actually having the DFT matrix. Perhaps if you were writing your own code you could leverage the symmetry of the DFT matrix so you only needed half of it but I don't think this idea will get you anywhere in Matlab.


Applying Convolution (In 2D) Theorem Swaps Quadrants

It's probably due to how you padded the filter before DFT'ing it.  I bet you had it in a corner instead of in the middle.  Since convolution via DFT is circular convolution, it's not so much of a quadrant swap as it is a wrap around.  I believe that if you put the non-zero taps in the middle that it will fix the problem.


Simple $\sin(2\pi 1000t)$ Fourier transform in PSpice not behaving as expected

PSPICE is showing the "Double Sided" Spectrum with the assumption (as it is) of a real signal. This means the spectrum is positive frequency axis only and the result will be double of what you would expect to see with a "Single Sided" Spectrum, basically matching the amplitude of the sinusoidal component represented.


How to Calculate Total Variation (TV) of an Image?

The Total Variation of an image $I$ can be computed using two formulas:

$TV(I) = \sum_{x} \| \nabla I (x) \|_1 $ (anisotropic TV);
$TV(I) = \sum_{x} \sqrt{ \| \nabla I (x) \|_2^2 }$ (isotropic TV).

In practice, both formulas yield almost the same result.

Using matlab, what you propose implements the first formula, while the second one can be obtained by sqrt(gradI_x.^2 + gradI_y.^2).


proving that a log gabor filter has 0 DC offset

Log-gabor are filters defined similarly as gabor filters in the sense that their envelope consist in a Gaussian in Fourier space. This is advantageous because this makes them optimal with respect to the compromise between localization (in space) and detection (of the mean frequency).

The difference is that log-gabor (as their name implies) are defined in the log-space frequency domain:

This makes sense as this relevant feature (frequency) may for some applications better optimized when the precision in frequency is proportional to the mean frequency.
this has the drawback that they have no simple analytical formulation in the space domain as simple Gabor filters,
In perception, this is very advantageous: for instance the human ear is in a large range sensitive to the relative increments of frequency (this is called the Fechner-Weber rule). This is also used in vision, see for instance the filters obtained in this this computational neuroscience paper mimicking the receptive field of simple cells in the primary visual cortex:

To specifically answer your question, log-Gabor filters have indeed the property of having a zero DC offset:

To show this consider that this offset is the value of the spectrogram at the origin (for null frequencies). In log-frequencies, the DC corresponds to the limit of the log function at zero, that is minus infinite. Necessarily, this implies that this value is zero for a finite bandwidth of the filter. In (simplified) mathematical terms, the envelope is :
$$f \propto \exp( - \log(f)^2 / B^2 /2$$
so that for $B>0$, we have  $f(0) \propto \lim_{f \rightarrow \infty} \exp( - \log(f)^2 / B^2 /2) = 0$.
This is relevant in perception. Indeed, remember that the DC component represent the average luminance value in the image, that is some global configuration (daylight vs nightlight, some specific exposure setting in your camera/eyeball) more relevant to specific tasks such as keeping the circadian rhythm (= syncing with the day rhythm). When doing a perceptually relevant task, such as detecting the shape of an object or tracking the motion of an object, these features are independent of the global configuration: e. g. a cat at night has the same shape as a cat during day!


Problem with a code of FIR Lowpass filter

The problem is, as you've already pointed out, the choice of the summation limits. Obviously, the array index shouldn't become negative. Look at the definition of convolution:

$$y(n)=\sum_{k=-\infty}^{\infty}h(k)x(n-k)\tag{1}$$

Now assume, as in your case, that the filter length (i.e. the number of coefficients) equals $M+1$, and that the input signal is a finite length sequence in the range $0\le n<N$. Then you have two conditions on the summation index $k$:

$$0\le k \le M\quad\text{and}\quad 0\le n-k \le N-1$$

With these conditions the sum in (1) becomes

$$y(n)=\sum_{k=\max\{0,n-N+1\}}^{\min\{M,n\}}h(k)x(n-k)\tag{2}$$

It is the sum in (2) that you need to implement.

Also note that if you want to implement convolution, then the output signal will be longer than the input signal, so you would need to define two different variables len. But if you only want to compute the output signal as long as there is a non-zero input signal present, then your implementation is fine.

Another improvement would be to have a routine that just takes a scalar input value (the current input sample) and computes the corresponding output sample. In this way you don't need to know the length of the input signal and you don't need to store the complete input signal in an array because you just process each incoming sample directly. Of course you need an internal buffer storing the last $M$ input values. Your version is appropriate for offline-processing of not too long signals.


What Information Is Lost by Histogram Equalization?

The expectation of information is called entropy. The loss of information can hence be understood as difference in entropy between source and processed image, assuming no random effect was added.

Together with my answer on why contrast is not an appropriate measure of entropy, this gives us the simple answer:

The loss in information is simply the number of bits necessary to count the different images that have the same processed image (i.e. after histogramm normalization).

Now, the math here gets pretty different whether you assume that pixels have real-valued intensities or discrete ones.

In the continuous case, you'd want to realize that histogram equalization applies a function to the pdf of the intensity, thereby adding negative differential entropy.

In the discrete case, multiple intensities of the source image will be mapped to a single intensity in the output: You can make statisitics of these ("empirical CDF") and from that derive how much information you lose.


You could represent the grid with 5 parameters:

(x0, y0) to represent the offset from your image origin to the grid
origin
(xs, ys) to represent a multiplication factor from your image
plane to your grid plane. If your grid is spaced apart every 100
pixels on your image plane, these values might be 100,100
theta to
represent the rotation of the grid

You could fit the grid to point in image by varying these 5 parameters and optimize these parameters by calculating a cost function at each step.

An example cost function would be the total Euclidean distance between known point and grid points. The grid parameters that best overlap your given point should represent an optimal solution in most cases, especially in cases where you have several "correct points" and few of the "outlier points".

For efficiency you could ensure that your x0 value is [0,xs] and your y0 is [0,ys]. Since it sounds like you have a guess for the size of your grid, you could input a decent guess for xs and ys and also constrain these values from changing too much.


Does using the Goertzel algorithm actually give better frequency resolution?

Your understanding is correct.

Goertzel algorithms give almost exactly the same result as 1 bin of a DFT or FFT of the same length or number of samples (and where the FFT twiddle factors are generated by a trig recursion), when used for frequencies that are exactly integer periodic in the Goertzel length.  But many forms of the Goertzel algorithm provide only the magnitude, not the complex result or the phase of the FFT 1 bin result.  Also, computationally, a generic Goertzel may be numerically slightly less stable than a generic FFT.  For non-integer-periodic-in-aperture frequencies, the result is equivalent to a sinc interpolation between bins of a DFT or FFT of the same length (which can be slightly more accurate of an interpolation than the more typical parabolic interpolation of FFT results).

The interpolation might be said to increase resolution in a graphical sense (more plot points), or making it easier to visually spot the maxima, but not in the information theoretic sense, nor for better separating 2 close spaced spectral lines as 2 separate peaks.


Speech dereverbaration via maximum kurtosis adaptive filtering

As was advised to me  by A_A in the comment section I will close this thread with a "self-answer".
The code present in the question is a code for a speech signal de-reverberation adaptive filter based on the Kurtosis signal of the signal's Lp residual. The original idea isn't mine (the reference is in the question). But the problem I faced was mainly due to the fact that my sensor did not convert to the desired result. I drastically increased the number of iterations the adaptive filter does and the result converged.(from n=500 to n=2000) There are others parameters that can be tweaked too (like mu) but I haven't had the time recently to thoroughly evaluate their effects.

Here's a link to download the code hope it helps. (expires by April 18th)

Code download link

Here's the link for the RIR code I used, I did not write the code though. I downloaded the file from matlab file exchange. The original author is Stephen G. McGovern, he also wrote a paper with the theoretical explanations for the RIR filter.

RIR code


Is it possible to find precise peak frequency using FFT phase shifts?

For a set of exact equations that describe how a pure real tone works in a DFT check out equations (23)-(25) in my blog article:

DFT Bin Value Formulas for Pure Real Tones

The phenomenon you have found can be seen clearer if you consider the complex definition of a real valued sinusiodal.

$$ \cos( \theta ) = \frac{e^{i\theta}+e^{-i\theta}}{2} $$

So, when you are near a peak bin, one of the numerator values is dominant, and you have an approximation of a complex tone.  Using equation (24) in this article shows how a the distance to the bin $\delta$ rotates the bin value.

DFT Bin Value Formulas for Pure Complex Tones

Hope this helps.

To answer your last question: Yes, I've done it different ways with more coming.

Check out:

Exact Frequency Formula for a Pure Real Tone in a DFT
Two Bin Exact Frequency Formulas for a Pure Real Tone in a DFT
Improved Three Bin Exact Frequency Formula for a Pure Real Tone in a DFT
A Two Bin Exact Frequency Formula for a Pure Complex Tone in a DFT
Three Bin Exact Frequency Formulas for a Pure Complex Tone in a DFT
[Jacobsen's Frequency Formula as an Approximation]
[Candan's Tweak of Jacobsen's Frequency Formula]
[Tweaking Macleod's Frequency Formula]

Your $k$ is the fraction of a bin distance from a bin on the frequency scale of cycles per frame.

Thus, from (20):

$$ k = -\frac{N}{2\pi} \delta $$

The negative sign is unfortunate, I should have done it differently in the article.

$$ \delta = -\frac{2\pi}{N} k $$

Your $\phi$ is the angle of the nearest complex bin value, so I will use $\theta$ in the formula for the phase angle from the signal definition in (24):

$$ \phi = \arg(X[]) = -\delta (N-1)/2 + \theta $$

$$ \phi = \frac{2\pi}{N} k (N-1)/2 + \theta $$

$$ \phi = \frac{N-1}{N} k \pi + \theta $$

So, for N being large, and $\theta$ being zero, $\phi$ looks like $k \pi$.

Remember, this is based on using a complex tone as an approximation for the real value tone, so it won't be exact.  I am not sure why your ratios are slightly larger than one when they should be slightly smaller.


How to find f given that f(some_image) = processed_image

Fit 3 polynoms (or any other model) using the individual channels of the output image: r' as a function of (r, g, b), g' as a function of (r, g, b), b' as a function of (r, g, b). I would also suggest not including in the data set all pixels with a 0 or 255 value (if your image is 8-bit) - these values might have been clipped.


Orthogonal signals in frequency and time domains

Preamble: a scalar product is a positive semi-definite bilinear form, and we will suppose here we want the standard definitions.

Using their Fourier domain is (one of) the simplest way to construct orthogonal signals. It suffices to choose signals with disjoint frequency supports. Note that is not a necessary condition: some signals can be orthogonal with the same frequency supports, like $\cos$ and $\sin$, and more generally Hilbert transform pairs.

This is due to the Parseval-Plancherel theorem (inner product and norm preservation under Fourier transform): with disjoint frequency support, the scalar product is zero, and hence the vectors are orthogonal.

So in the following results, we take a random spectrum for a real signal. We then segment it into three subsets of points, mutually exclusive, and whose union give the full frequency axis. And then do an inverse FFT to obtain three subsignals whose inner products are:

S1.S2 = -4.7705e-18
S2.S3 = 1.7564e-17
S3.S1 = -1.2645e-17

The corresponding signals and spectra are displayed below, and the code follows (in Matlab):

function seDsp51117(nSample)

if nargin < 1
    nSample = 512;
end
nSampleHalfSpectrum = nSample/2+1;
halfSpectrum = randn(nSampleHalfSpectrum,1)+1i*randn(nSampleHalfSpectrum,1);
halfSpectrum(1) = real(halfSpectrum(1));
halfSpectrum(end) = real(halfSpectrum(end));
nSampleSpectrum1 = ceil(nSampleHalfSpectrum*rand(1));
nSampleSpectrum2 = ceil((nSampleHalfSpectrum-nSampleSpectrum1)*rand(1));
nSampleSpectrum3 = nSampleHalfSpectrum - nSampleSpectrum1 - nSampleSpectrum2;

idxSpectrumRand   = randperm(nSampleHalfSpectrum);
bndSpectrumRand = round(linspace(1,nSampleHalfSpectrum+1,4));

idxSpectrum1 = idxSpectrumRand(bndSpectrumRand(1):bndSpectrumRand(2)-1);
idxSpectrum2 = idxSpectrumRand(bndSpectrumRand(2):bndSpectrumRand(3)-1);
idxSpectrum3 = idxSpectrumRand(bndSpectrumRand(3):bndSpectrumRand(4)-1);

[signal,halfSpectrumSelect] = IFFTReal(halfSpectrum);
[signal1,halfSpectrumSelect1] = IFFTReal(halfSpectrum,idxSpectrum1);
[signal2,halfSpectrumSelect2] = IFFTReal(halfSpectrum,idxSpectrum2);
[signal3,halfSpectrumSelect3] = IFFTReal(halfSpectrum,idxSpectrum3);
figure(1);clf
subplot(4,2,1)
plot((signal));axis tight;grid on
xlabel('Original signal')
subplot(4,2,3)
plot((signal1));axis tight;grid on
xlabel('Subsignal S1')
subplot(4,2,5)
plot((signal2));axis tight;grid on
xlabel('Subsignal S2')
subplot(4,2,7)
plot((signal3));axis tight;grid on
xlabel('Subsignal S3')

subplot(4,2,2);hold on
plot(abs([halfSpectrumSelect]),'o');
plot(abs([halfSpectrumSelect]),'x');
xlabel('Original spectrum')
axis tight;grid on
subplot(4,2,4);hold on
plot(abs([halfSpectrumSelect]),'o');
plot(abs([halfSpectrumSelect1]),'x');
axis tight;grid on
xlabel('Subsignal spectrum for S1')

subplot(4,2,6);hold on
plot(abs([halfSpectrumSelect]),'o');
plot(abs([halfSpectrumSelect2]),'x');
axis tight;grid on
xlabel('Subsignal spectrum for S2')
subplot(4,2,8);hold on
plot(abs([halfSpectrumSelect]),'o');
plot(abs([halfSpectrumSelect3]),'x');
axis tight;grid on
xlabel('Subsignal spectrum for S3')

sp12 = dot(signal1,signal2);
sp23 = dot(signal2,signal3);
sp31 = dot(signal3,signal1);

disp(['S1.S2 = ',num2str(sp12)]);
disp(['S2.S3 = ',num2str(sp23)]);
disp(['S3.S1 = ',num2str(sp31)]);

function [signal,halfSpectrumSelect] = IFFTReal(halfSpectrum,idxSpectrum)
if nargin < 2
    idxSpectrum = 1:length(halfSpectrum);
end
halfSpectrumSelect = zeros(size(halfSpectrum));
halfSpectrumSelect(idxSpectrum) = halfSpectrum(idxSpectrum);

%% Hermitian glue
wholeSpectrum = [halfSpectrumSelect;conj(flipud(halfSpectrumSelect(2:end-1)))];
ratioImagReal = 10^(-6);
tmp = ifft(wholeSpectrum);
normRatio =  norm(imag(tmp),1)/norm(real(tmp),1);
    signal = real(tmp);
if normRatio > ratioImagReal
warning('Check the imaginery value');
    end


1) Normalize your image to range $[0,255]$.

2) Select a threshold and threshold the image. For your image, what worked is: $\tau=[140-150]$.

3) Compute a Euclidean distance transform.

4) Apply watersheds segmentation.

If I apply this procedure, here is what I get:

Not perfect, but maybe a good start. The result looks similar to performing a Voronoi diagram on a smoothed/noise filtered version of the image.

Now to each region, one could fit a circle, using any algorithm (edge extraction and fitting or Hough transform etc.) For now, I used a simple circle fitting to the edges:

Again, not the best, but could be a sufficient starting point.

I also think that you would find the following useful:

  Segmentation Methods for Light Microscope Image

  Detecting Spheric Particles


How to predict the cramped frequency of a digital filter based on an analogue frequency?

I've come up with a solution.

Let $ f_a = $ analogue frequency, $ f_d $ = digital frequency, $ f_s = $ sampling rate and $ f_c = $ corner frequency with $ \omega_c = \frac{\pi f_c}{f_s} $ and $ \omega_a = \frac{\pi f_a}{f_s} $ then:

$$ c = \frac{\omega_c}{\tan(\omega_c)} $$

$$ \omega_d = \arctan\left(\frac{\omega_a}{c}\right) $$

Giving:

$$ f_d = \frac{\omega_d f_s}{\pi} $$


That $\angle$ is "angle". As the text right below it indicates, it's the phase of a complex number, in this case the phase of the sum of squared $p_k$.

Your second question: I really don't think pseudocode makes this any clearer, but sure:

factor = 1
sequentially for each timestep n:
    s = r[n] * factor
    t = integrate_over_last_N_and_decimate_by_N(s)
    u = remove_bits(t)
    v = non_coherently_integrate_over_last_K(u)
    w = non_linear_discriminator(v)
    x = loop_filter(w)
    y = nco(x)
    factor = generate_carrier(y)

I mean, that's a pretty straightforward block diagram.


Where to get transtion matrix for Kalman filter?

Simple: if you do not have a model, you cannot apply the Kalman filter. Or you could and make up a model, but you cannot expect any of the optimality properties of the filter to hold. Based on that paper, it seems they assume that $\mathbf{x}_k$ is known or at least directly measureable...


How to determine the original size of the FFT used in unknown OFDM symbol

If you use an FFT size to generate a PSD estimate far longer than the IFFT used to generate the OFDM signal, then you'll see the "valleys" between the sinc maxima in the spectrum, and then it's easy to just go in there and count subcarriers, and their spacing.

This will, however, require a high degree of frequency stability in your system, and won't work well with OFDM systems where small number of consecutive symbols are bursted out.

You can make that a little more robust and work on longer autocorrelation functions, like here.

You probably want one of the more advanced methods of OFDM parameter estimation, like those implemented in gr-inspector (hint: that's free&open source; don't reinvent the wheel).

The relevant blog post of the author of that software  would be https://grinspector.wordpress.com/2016/07/08/week-7-ofdm-prototype/ .

You basically want to read

  Bouzegzi, A., Ciblat, P., & Jallon, P. (2010). New algorithms for blind recognition of OFDM based systems. Signal Processing, 90(3), 900–913. doi:10.1016/j.sigpro.2009.09.017

They propose clever (as in: outperforming above roughly sketched autocorrelation-based methods) ways to synchronize and estimate OFDM parameters.


Transfer Function to go from a Step Input to a Linear Ramp between the two step values

A simple moving average with a rectangular window will turn a step into a ramp, but no linear filter can give you a constant slope, because linearity implies that scaling the input (step) will scale the output (slope). Instead, what you want is a non-linear slew-rate limiter that steps towards the target at constant rate (eg. take the minimum of the "distance to target" and your allowable rate, then accumulate).
If you never want to move "too slowly" then you likely need a state-machine. If you are in steady-state and the target has moved "far enough" from your current value, then transition to moving state immediately. When you reach the target, transition back to steady-state. If you are in steady-state but the target only moved slightly, then transition into a timer-state where you still check if the target is far enough to start moving immediately, but also start moving when the timer expires. The downside is that small steps (or shallow ramps) will cause you to wait for the timer, but this is unavoidable unless you want to constantly (in the worst case every time-step) switch between moving and not moving (at which point you are basically doing PWM).
If you want to smooth the acceleration, then this can be done by taking the output of the slew-rate limiter and filtering it further. One (or two if you also want to avoid jerk) passes of moving average with short rectangular window often works well for this and will give you symmetric acceleration vs. deceleration and will still reach the target in finite time. Exponential moving averages (aka. one-pole low-pass filters) can also be used, but in that case you will only reach your target asymptotically.


What do high and low order have a meaning in FIR filter?

1- an FIR filter can be described by the following input/output relation:
$$ y[n] = \sum_{k=0}^{N} b[k]x[n-k]$$

where $N$ is the order of the FIR filter (its length-1) and those coefficients $b[k]$ of length $N+1$ are the FIR filters impulse response, or the filter coefficients in practice.

As one can see high order FIR filter means $N$ to be large and low order means $N$ to be small. Note that most of those coefficients can be set $0$ but the last one at $x[n-N]$ would still make a high order filter (specifically a filter with a pure delay).

The delay associated with FIR filters (of symmetric impulse response types) is directly related to the order $N$ of the filter.

2-
An FIR low-pass filter is a mathematical algorithm and its cutoff frequency is normalized to be between $0$ and $\pi$ radians mathematically and between $0$ and $1$ for software such as MATLAB. When you want to find the physical cutoff frequency of such an FIR filter , you should consider the physical sampling rate and relate it through $$f_c = \frac{\omega_c F_s}{ 2\pi} $$
where $\omega_c$ is the cutoff frequency (in radians per sample) of the mathematical filter, and $f_c$ is the physical cutoff frequency in Hertz.

3- Impulse response of the FIR filter is directly available in the coefficients $b[k]$ as I've explained. When you are dealing with an IIR filter or an LTI filter which does not posses a LCCDE expression, then practical computation of the impulse response that can be done in a number of ways. The simplest could be to create a signal that defines the unit impulse as delta = [1, zeros(1,L-1)] then the $L$ points of the impulse response can be computed directly by h = filter(b,a,delta) . Note that if impulse response $h[n]$ does not decay fast enough then you should increase the length $L$.


If by "performance" you mean the quality of the processed (speech) signal, then this question cannot be answered. The reason is that the quality of the output of a spectral subtraction algorithm depends on so many factors, and the precise subtraction rule is just one detail. What is much more important is the noise estimation algorithm, certain smoothing operations, etc.

What you can say is that the dynamic range of algorithms working with magnitudes is much smaller than when powers are used. This can be a great advantage when such systems need to be implemented in fixed-point arithmetic.


How to implement compressed sensing reconstruction?

I have had encountered this problem before. In using l1-magic toolbox, the first thing you need to do is to modify the file l1eq_pd.m.

To do so find all linsolve instructions and remove the third parameter, for example on line 141 change the following line from

[dv,hcond] = linsolve(H11p, w1p,opts);

to,

[dv,hcond] = linsolve(H11p, w1p);

  This is because the structure opts specify the property of the matrix
  A in solving Ax=b linear system of equations and omitting the
  structure opts in the linsolve function will allow Matlab to choose
  the appropriate solver automatically [Reference]

If the errors persists, your own code should be scrutinized.

I suggest to have a look on this examples codes.


fundamental frequency contour of a syllable

Syllables might require a certain duration to be recognizable by humans.  Over part of that time, voiced syllables might have a pitch, or a clear periodicity in waveform.  Pitch can be described as having a fundamental frequency or F0 (the reciprocal of the perceived repetition period). Note that there a lots of other answers ont this Q&A site on various ways to measure pitch, and the potential problems related to such.

Pitch can change over time.  If you measure pitch at several points in time (say 9 points for example) over the duration where a syllable is voiced, you might find that the measured pitch in fact has changed (or not).  With those pitch measurements over time one could interpolate a curve (spline or polynomial regression, etc.) thru those measured points.  You might end up making a bunch of different looking curves from different snippets of audio.  One could also publish nice looking graphs of those interpolated curves and call the shapes of those curves contours.

Feed a lot of those nice looking pitch/F0 frequency contours (plus a large amount of other hopefully related data) to a machine learning algorithm, and it might be able to infer some possibly interesting correlations between things.

Or not.


Let me try to give you a short introduction to the concept of bi-orthogonal bases from a mathematical point of view. You are already using the concept, but your notation and nomenclature is obfuscating what is really going on. I will make this answer independent of the wavelet context you have, so that it is more general. You will certainly be able to apply it to your situation however.

You have a non-orthogonal basis, i.e. a set of linearly independent vectors spanning your $L^2(K)$ Hilbert space $\mathcal{H}$. That vector space induces another vector space, that of all linear forms $F:\mathcal{H}\to K$. We can make this space a dual space of $\mathcal{H}$ if we take the Hilbert space inner product $\langle \cdot, \cdot\rangle : \mathcal{H} \times \mathcal{H} \to K $ to generate a linear form $\langle \phi, \cdot \rangle : \mathcal{H} \to K$ for every $\phi \in \mathcal{H}$.

In an effort to simplify the notation a bit and make things more intuitive, we can introduce Dirac's Bra-Ket notation for vectors and their duals. Instead of $\phi \in \mathcal{H}$, we write the Ket $\left| \phi \right\rangle$ and instead of $\langle \phi , \cdot \rangle \; \phi \in \mathcal{H}$ we write the Bra $\left\langle \phi \right|$. Applying the linear form $\left\langle \psi\right|$ to $\left| \phi \right\rangle$ is then simply the notational fusion to a braket standing for the inner product $\left\langle \psi | \phi \right\rangle$.

In the bi-orthogonal basis business, this distinction is often ignored, making everything more confusing that it should be. Consider a basis that is not orthogonal with respect to the canonical inner product. We can of course use the usual knowledge of linear algebra to still construct everything of interest and expand vectors with respect to this basis, but it tends to get messy quickly.

That's where a neat trick comes in. We simply invent a new inner product, i.e. a non-degenerate, non-negative bilinear form $B: \mathcal{H} \times \mathcal{H} \to  K$, so that $B( \phi_n, \phi_m ) = \delta_{n,m}$ for your basis $\{ \phi_n, n\in\mathbb{N}\}$ and the Kronecker symbol $\delta$. Because the basis is complete, $B$ is entirely determined by this condition. We can also write this bilinear form using a linear operator $B:\mathcal{H}\to\mathcal{H}$ with the same name. Then we get $$\left\langle \phi_n | B | \phi_m \right\rangle := \left\langle \phi_n | B \phi_m \right\rangle = B( \phi_n, \phi_m ) = \delta_{n,m}$$.

At this point, it makes sense to define the adjoint $A^\dagger$ of a linear operator $A$ like follows: $$\left\langle A^\dagger \phi_n | \phi_m \right\rangle = \left\langle \phi_n | A \phi_m \right\rangle$$

Of course it follows from the definitions above: $$\left\langle B^\dagger \phi_n \right| = B(\phi_n,\cdot)$$
and therefore $$\left\langle B^\dagger \phi_n | \phi_m \right\rangle=\delta_{n,m}$$
which means $\left\langle B^\dagger \phi_n \right|$ is the dual of $\left|\phi_n\right\rangle$ under $B$, whereas the dual under the canonical inner product is still $\left\langle \phi_n \right|$.

Now what you can often see, and I cannot stress enough how much I think this is a bad idea, is that the dual under $B$ is mapped back to $\mathcal{H}$ as a dual under the canonical inner product. That means you construct a second basis $\left| \tilde{\phi}_n\right\rangle =\left| B \phi_n \right\rangle$ which is bi-orthogonal to the first basis with respect to the canonical inner product:
$$ \left\langle \tilde{\phi}_n | \phi_m \right\rangle = \delta_{n,m}$$
This obfuscates the real meaning of the second basis and puts both on the same footing and suggests that both bases live in the same space and transform in the same way, which they don't.

So, what is the alternative? Not taking the $B$-dual back to $\mathcal{H}$ makes things much more intuitive. You want to expand $\left|\psi\right\rangle$ in the basis $\phi_n$? Just use the $B$-dual: $\left|\psi\right\rangle=\sum_n \left| \phi_n \right\rangle \left\langle B^\dagger \phi_n | \psi \right\rangle $. Or you want to project onto a vector $\left| \kappa \right\rangle$? Just use the $B$-dual: $\left| \kappa \right\rangle \left\langle B^\dagger \kappa | \psi \right\rangle$. So whatever you do, just use $B$ as your inner product and you can work like with an orthogonal basis.

And to come back to your original question. You say that $V_0$ is not orthogonal to $\oplus_n W_n$. Yes, but only with respect to the canonical inner product. With respect to $B$, they are orthogonal. And so your argument for the orthogonal basis case works.

I hope all this does not confuse you. I've tried to be as clear as possible, but it's not a trivial topic, so don't be frustrated if you have to read the answer several times until it makes sense. If I can clarify anything, please let me know.


Why does raised cosine pulse through Butterworth filter introduce (significant) ISI?

That's because Matlab's rcosdesign creates a root-raised-cos filter, not a raised cos filter. If you want no ISI (ie meet Nyquist's criterion), you need to do the match filtering of the reception part :

%%
alpha = 1;
R = 10;
hcos = rcosdesign(0.25, 30, R);
x = upfirdn(sign(randn(1000,1)), hcos, R);
y = upfirdn(x, hcos, 1, R);
eyediagram(x, 2*R);
hold on
eyediagram(y, 2*R);

eye of x

eye of y


How to verify the obtained SNR of the noisy signal after adding noise?

Here is a common/straight-forward way using the discrete time domain samples. If you have the noise free signal, $x[n]$, and you created a noise signal, $w[n]$, then you can calculate the SNR by using the formula:
$$\hat{\text{SNR}}=\frac{\sum_n |x[n]|^2}{\sum_n |w[n]|^2}$$
Edit
An unstated assumption above is that $\mathbb{E}\big[x[n]\big]=\mathbb{E}\big[w[n] \big]=0$. In general, you'd use the variance not the sum of squared values, as to take care the possible non-zero means:
$$\hat{\text{SNR}}=\frac{\text{var}(x[n])}{\text{var}(w[n])}$$


In order to shift in MATLAB you need to play with the indices of the signal vector.
for example, given vector x - vX of length L, to shift it in time such that the first sample is the fifth you'll do this:

vXShifted = x(5:L);

Pay attention that the signal is shorter by 4 samples.

Now, sometimes, for periodic signals, a cyclic shifting is required which can be done using circshift() function in MATLAB.


You can notice that:

the top left looks like the original image, albeit smaller
the top right shows mainly horizontal features
the  bottom left shows vertical trends
the  bottom right is much less clear

What you show "is not wavelets" per se. It is another image, normally with the size of the original one, with four quadrants. However, it illustrates some wavelet features. The top left is a coarse approximation of the image, resulting from filtering and downsampling,  obtained from a scaling function. The three others are details, with at least one high pass component. In the classical version, this image combines low- and high-pass on rows, and the same on columns, each followed by downsampling on both directions, which explains the above items. If you iterate the process on the low-pass image, under some conditions, you get a separable two-dimensional wavelet transformation. The key point is that in your image, all the original information remains.


Is a high-passed signal the same as a signal minus a low-passed signal?

In general you can't simply subtract a low-pass filtered version of a signal from the original one to obtain a high-pass filtered signal. The reason is as follows. What you're actually doing is implement a system with frequency response

$$H(\omega)=1-H_{LP}(\omega)\tag{1}$$

where $H_{LP}(\omega)$ is the frequency response of the low-pass filter. Note that $H_{LP}(\omega)$ is a complex function. What you probably want is

$$|H(\omega)|=\big|1-|H_{LP}(\omega)|\big|\tag{2}$$

but this is generally not the case when (1) is satisfied.

Now write $H_{LP}(\omega)$ as

$$H_{LP}(\omega)=|H_{LP}(\omega)|e^{j\phi(\omega)}$$

where $\phi(\omega)$ is the phase response of the low pass filter. If you defined a filter as

$$H_{HP}(\omega)=e^{j\phi(\omega)}-H_{LP}(\omega)=e^{j\phi(\omega)}\left(1-|H_{LP}(\omega)|\right)\tag{3}$$

then this new filter would indeed satisfy the magnitude relation (2). So if you do the subtraction as in (3), you can realize a high-pass filter with a magnitude response that is complementary to the response of the low pass filter. This means that you need to filter the signal by an all-pass filter with frequency response $e^{j\phi(\omega)}$ before subtracting from it the low-pass filtered version of the signal.

In practice this is very simple if the low-pass filter has a linear phase response, because then the phase term is given by

$$e^{j\phi(\omega)}=e^{-j\omega\tau}\tag{4}$$

and corresponds to a simple delay. FIR filters can have an exactly linear phase response, and the ones designed by the Parks-McClellan algorithm do have a linear phase. It is advantageous to make the delay $\tau$ in (4) an integer, which is achieved by choosing an even filter order $n$: $\tau=n/2$

So what you have to do is the following:

design a linear phase FIR low-pass filter with an even order
filter the signal using filter(); let's call the result $x_{LP}[n]$
delay the original signal by $\tau=n/2$ samples; let's call the delayed signal $x_d[n]$
generate the high-pass filtered signal by subtraction: $x_{HP}[n]=x_d[n]-x_{LP}[n]$

Here is a very simple illustration in Matlab/Octave

h_lp = fir1(100,.3);                        % low-pass design
h_hp = [zeros(50,1);1;zeros(50,1)] - h_lp;  % high-pass design by subtraction
[H_lp,w] = freqz(h_lp,1,1024);
[H_hp,w] = freqz(h_hp,1,1024);
plot(w/2/pi,20*log10(abs(H_lp)),w/2/pi,20*log10(abs(H_hp)))
grid, axis([0,.5,-100,5])

EDIT:

Since you used the filtfilt command, the phase is artificially eliminated, and Equations (1) and (2) above become equivalent, because all frequency responses are in fact squared magnitudes of the designed responses. So, apart from small differences due to the filter design process, numerical errors, and small differences caused by the filtfilt function (automatically chosen initial conditions minimizing start and end transients), the result of subtracting filtered data from the original data should closely resemble direct filtering with a complementary filter. Since this is not the case in your example, I suspect that the filter design routine gives you trouble due the extremely high filter order. If you do the same with simpler filters (I chose $n=100$), you get what you would expect. In the figure below you see a section of the data in blue, the output of the low pass filter in green, and the result of the subtraction of the output of the high pass filter from the original data in red. The green and red curves are virtually identical.

x = load('data.txt');           % data to be filtered
h_lp = fir1(100,.3);            % LP impulse response
h_hp = fir1(100,.3,'high');     % HP impulse response
y = filtfilt(h_lp,1,x);         % apply low pass filter
yh = filtfilt(h_hp,1,x);        % apply high pass filter
yd = x - yh;                    % low pass by difference with high pass filter
n = 1:length(x);
plot(n,x,n,y,'g.',n,yd,'r')
axis([3500,4000,140,150])


This is somehow a broad question. An interesting book is Digital Processing of speech signals by Rabiner & Schafer, in which you can find some good explanations to get started. It probably is out of date when looking into more advanced techniques, but I think it's good to understand the main characteristics of speech.

In short, every sound we produce is determined by the disposition of the vocal tract. In particular, if the vocal folds are tense and don't vibrate, the sound is said to be unvoiced. However, if they vibrate, the sound is voiced. You can sense the difference by placing your hand in your neck while you say "Ahhhh" or "Ssssshh". In Spanish, all vowels are voiced and all consonants are unvoiced (this is not totally true, please refer to the comments). I'm not sure about other languages though. You can think about a generator (glottis) which is connected to a filter (vocal tract) in a way that sounds take their characteristic shape. This is in fact the basis for speech production models (such as the simplest).

If you look closely to a speech signal, you can easily find differences between the different phonemes. Voiced phonemes are said to be pseudo-periodic. This means that, if you look at a fragment of the signal which is short enough (about tens of milliseconds long, e.g. 40ms), you can find the superposition of harmonic sine signals. This means that:

The signal is periodic with a fundamental frequency $f_0$. This is also known as the pitch.
Harmonics will appear with frequencies $2f_0, 3f_0, 4f_0$... which energy will depend on the exact phoneme.

If you plot the spectrum of this speech signal, you can find something like the following (source).

This resembles, for example, the harmonic structure of a brass instrument sound.

In contrast, unvoiced signals don't have such a clear harmonic structure. In fact, they have no fundamental frequency, and can be conceived as a noisy signal. You can sing by saying "Ahh", but you can't by saying "Ddd". Unvoiced sounds can be conceived as coloured noise. They usually have less energy than voiced.

Taking these differences into account, we can think about a way to detect them automatically.

  Voiced sounds are periodic while unvoiced aren't.

To detect periodicity, you can use methods such as autocorrelation (which has been widely treated in this site).

  Voiced sounds have more energy than unvoiced.

You can divide the signal in 20-40ms windows, compute the energy for each one and decide according to a threshold.

  Unvoiced sounds are noisy and tend to be prominent in high frequencies, while voiced aren't so.

Therefore, unvoiced will cross with zero more often than voiced. This is called Zero Crossing Rate.

While these simple methods are not too accurate by themselves, are a good place to start from. Furthermore, I am sure you can get good results by combining them. In any case, this post does not substitute a good book, as the one I mentioned before. It's just some kind of introduction to show you some possibilities.


Filter, window function and audio signal: how to link them to each other?

Your sample code looks like it is trying to create a windowed Sinc FIR filter.  A Sinc function (the transform of a "brick-wall" rectangle filter) has infinite support, thus it must be windowed to create a finite length FIR filter kernel.  There is pseudo-code to generate such a filter, using a Hamming window, on my DSP web page.

You can apply such a FIR filter to audio by convolution with your audio samples (either directly convolution, or using FFT/IFFT overlap-add/save fast convolution on blocks of data).

If you are using iOS frameworks such as Core Audio, note that you can also use the vDSP_biquadm API in the iOS Accelerate framework for IIR filtering.


Averaging magnitude squared coherence across multiple time series

I believe I might have a solution... Editing A Donda's work above, I think a simple rearrangement of the last two lines is necessary:

sxya = (abs(sxy1) + abs(sxy2)) / 2;
ca = sxya.^ 2 ./ sxa ./ sya;

So here, the average of the magnitudes is taken, rather than the average of the complex CPSDs. Then, this average can be squared.


Wrong values calculating FFT with EEG Bands using Numpy

There are a few things going on:
1.) Your data is actually only length 5 (made up of 5 longer components). This means when you call fft_freq = np.fft.rfftfreq(len(data), 1.0/fs) your array containing the frequencies of the FFT is only length 3! I think the main issue I'm having is that I'm not sure what your data actually represents: multiple electrodes?
2.) I think you're using np.where wrong. Basically it's not returning a slice which appears to be what you're trying to get.
3.) Even if a single array in your data represents an FFT of a single channel, the length is not nearly long enough to get anything useful. If we take the frequencies of data[0] we get [  0.   12.5  25.   37.5  50.   62.5  75.   87.5 100. ]. Not even one bin per band so you need to increase the length of your data. This is probably why freq_ix returns an empty array sometimes.
I was basically doing the same thing here in lines 45-68. A little different but more or less the same idea.
Also highly recommend using the MNE package. :)


Can you recommend an efficient encoding algorithm for sound with long periods of silence?

Thresholding should reduce the file size to the limits of the encoder being used. AAC/H.264 are using the state of the art general purpose audio compression codecs (which could only be beaten by application specific optimized codecs, such as speech vocoders etc).

May be your files are already minimal? Put some numerical data such as the sampling rate, average bits per sample, recording duration and file sizes...


Histogram Equalization - mapping a certain portion of pixels to pure black and white

Compute the image histogram, then its running integral (i.e. an array T with as many entries as there are gray levels, such as T(k) = sum_{i=0, 1..k} (histogram(i)).

The last element of the integral is, of course, the area of the image.
For the lower threshold, start from the left and find the level at which T exceeds your desired fraction. Analogously, for the high threshold, start from the right and find the level at which T goes below your high fraction.


You forgot to mix or heterodyn the complex FFT output up to the carrier frequency.  Then the negative index FFT bins will be below the carrier in frequency.  FingerIO appears to use a very high audio frequency carrier.


How to find the time of beginning of the slope of the signal?

I observe the following about your signal:

of interest are only frequencies below 20 kHz
there's but a single, very dominant tone in that – you've got excellent SNR
The development of frequency over time is either constant, or a linear function of time

So, from that, I'd propose the following steps:

Low pass filter appropriately to restrict your bandwidth to below 15 to 20 kHz
employing a parametric spectrum estimator to find the frequency of your single tone
post-process that frequency estimate to decide whether you're still in the linear or in the constant frequency regime.

Regarding these steps:

might be easiest with just a simple second-order PLL. If you don't know where your tone is going to be at the start of your observation, start with a very high loop filter bandwidth, and as soon as you have a lock, reduce the bandwidth drastically. That second-order control loop actually directly gives you an error signal that you can use for 3.
that's a problem for a bit of low-pass filtering and estimating the frequency/time slope from that, and then checking whether the hypothesis still fits, or the slope has suddenly changed.


You are correct. However, the coefficients $X_k$ could be samples taken from an analog signal. These could be used in the receiver to recover the analog signal using interpolation.


Your bipolar sequence at one sample per symbol is equivalent to a perfectly synchronized sequence of symbols with no intersymbol interference (ISI). Intersymbol interference in this context refers to the tendency for symbols that precede and follow the current symbol to interfere at the receiver. That is, at the correct sampling instant for symbol $N$, the receiver also sees some of symbols $N-1$ and $N+1$ (and possibly more in each direction) mixed in. This causes an increase in symbol error rate.

In general, ISI is caused by some non-ideal channel impulse response between your transmitter and receiver. This can be caused by several phenomena, including:

A multipath communications channel (e.g. a wireless channel). In this case, the receiver observes a signal that can be modeled as your originally transmitted signal convolved with some impulse response.
The use of a non-ideal pulse shape. On the AWGN channel, the optimum receiver consists of a filter whose impulse response is matched to the transmitted pulse shape. The receiver therefore observes a stream of pulses where each has the shape of the originally transmitted pulse shape convolved with itself. The choice of the pulse shape is important to ensure that the receiver doesn't encounter any ISI after applying the matched filter.

This is why the root-raised cosine filter is often used. When convolved with itself, the result is a raised-cosine response, which has zero ISI (it is exactly zero for all multiples of the symbol period $T$, apart from $t=0$).

With that said, there are some practical considerations when designing your pulse shape. When you actually transmit it, you'll need to oversample it in order to meet the sampling theorem's reconstruction requirements. A typical oversampling rate might be 4.

The rolloff factor controls how much bandwidth the RRC-shaped signal will occupy. In general, the more excess bandwidth, the better you can expect your synchronization algorithm to perform (i.e. you would usually observe smaller error from your timing recovery subsystem with a larger rolloff factor). However, this comes at a cost of more bandwidth required to carry the signal.


I've used the first technique before. The idea is that you have some chunk of time domain samples where you know only noise is present and use that to get an estimate of the noise power, $P_n$. Then you determine the total received power, which will have contributions from both signal and noise, $P_{s+n}$. When you take the ratio of these, you get:
$$\frac{P_{s+n}}{P_n}=\frac{P_s}{P_n}+1$$
For high SNR, the extra $1$ doesn't matter too much but to be exact the SNR estimate should be:
$$\hat{\text{SNR}}=\frac{P_{s+n}}{P_n}-1$$


Projecting a 3D point into new camera coordinates

Your formula is correct. Here are some precisions:

the matrix $P$ takes a scene point to an image reference coordinate frame;
the matrix $K$ takes the coordinates of the point in the reference image frame and makes the necessary frame changes to fit with the actual camera (origin of the frame, pixel size...);
$X$ is a 3D point in projective coordinates, i.e., $X = (x, y, z, 1)^T$, otherwise the dimensionality of the point is not correct with the projection matrix;
$x$ will be actually a column of 3 rows: $x = (u,v,w)^T$, and the final image coordinates are given by diving $u$ and $v$ by $w$.

A few answers to your comments:

the equation assumes an ideal pinhole camera model. Thus, the coordinates are undistorted. You can actually find all these equations by yourself using basic triangular geometry and Thales's theorem;
you don't add the projective coordinate ("the 1") twice! The product $P \times (X,1)^T$ will create as output a projected (2D) point that has 3 coordinates, i.e., a 2D point in projective coordinates. What you need only is to "remove" the projective coordinate by dividing by $w$ at the end;
to get the order always right, the best way is to think of these equations as I wrote then in plain English. You take a 3D world and project it onto the world plane given by the camera position ($P$, extrinsic parameters). Then, you consider what happens for this particular camera model/sensor (intrinsic parameters). The extrinsic parameter matrix describes the camera with respect to the world, while the intrinsic parameters are always teh same wherever you put your camera.


Entropy of a dynamical system and source code length

I can try to answer your first question (but I'm not sure I understand exactly what you mean).

A source produces messages with average information $H$. You seem to have a binary source, with two messages $m_1$ and $m_2$.

Now you want to encode the messages to a sequence of 1s and 0s (subsequently called "bits"). This is called "source coding". Even though we know that, in principle, we should be able to encode $N$ messages using $NH$ bits (for large $N$), it is not always obvious how to do so. There are algorithms such as Huffman and Lempel-Ziv that can find good (even optimal) encodings. Note that all good coding algorithms encode messages in groups, not individually, and they approach $NH$ only when $N$ is large.

If your question is "what's the smallest number of bits that I need to encode $N$ messages from a source with entropy $H$", the answer is $NH$ (in theory), and somewhat more in practice, depending on the optimality of your encoding algorithm and the value of $N$.

Example 1: A source produces 8-bit messages $b_7b_6b_5b_4b_3b_2b_2b_1b_0$ uniformly at random, but $b_7$ is always 0. The entropy is 7 bits/message. Let's say we want to transmit 100 messages. Even though the source produced 800 bits, we only need 700 to transmit the messages.

Example 2: You flip a biased coin with $p(\text{heads})=0.3$ 1000 times. The entropy of the experiment is 0.88 bits/flip. To convey the result of the 1000 flips, you'll need on average 880 bits. Some particular results will need more bits, some will need less. To figure out how many bits you need for a particular sequence of flips you need to specify your encoder algorithm.

Example 3: You have a biased 3-sided die with probability distribution $p(1)=0.5$, $p(2)=0.25$, $p(3)=0.25$. The entropy is 1.5 bits/throw. It's pretty easy to come up with an encoder:

result  encoding
1    -> 0
2    -> 10
3    -> 11

You should be able to convince yourself that, on average, this encoding results in 1.5 bits per throw. However, let's say that you're unlucky and you get 1000 '3' in a row. That particular sequence will need 2,000 bits to transmit. That's why the entropy can only predict the average sequence length, not the length you'll get for a particular set of messages.


What does frequency weighting an acceleration signal mean? And How do I frequency weight the raw acceleration data collected?

...I do not know where this comes from...

It looks like that's not a scientific article, so it doesn't have references at the end.  Read it carefully, and see if it references military or industry standards for vibration dosage.  Then read those standards and see if they refer back to any studies.
It'll be "fun", like black-water plumbing.  Think of it as the intellectual equivalent of cleaning a blocked sewer line.

I do not understand what the plot means and how to get weightings from it. But when I think about it, I feel like it has something to do with the acceleration data collected over 0.5 - 80 Hz and get the weightings from there. The thing is the weightings are not constant in the plot for the mentioned frequency region.

Probably you get the frequency spectrum of your accelerations in $x$, $y$ and $z$, then you multiply each of those by the individual weightings shown, then you calculate your $VDV$ from that.
I was going to say you could use Parseval's Theorem and just do the computation in the frequency domain -- but there's that $a^4$ term in there, so you may just need to weight it in the frequency domain, convert to time domain, and do the math.
But this smells strongly of being an established field, which means it has established methodologies; you should not invent your own method.   The established methodologies are often different from what an expert in signal processing would do, in detail if not in principle, and sometimes that different can be significant.  Even if your own special way is better in every regard, it'll still be "wrong" because it'll be different from how everyone else does it.
So find the source material via a literature search: either dig into the scientific literature (starting with the references in your second paper), or if you're trying to meet some regulation into the relevant industry or government standards.  Somewhere in all the muck the method will be spelled out.  If it's a regulation, it'll be written in techno-legal, and chances are you'll need to find an explanation of it someplace.  If you can't understand the actual defining document, then you can search on the title or designation of that document, i.e. "how do I implement SAE-234566-3?".

Also, how do I use the weightings with the raw-data collected to get frequency-weighted acceleration data?

Explained above -- convert your measurements into a spectrum, multiply said spectrum by the weightings in your graph.

Do I multiply the discrete raw-data with the the corresponding weights for each axis to get frequency weighted acceleration?

Yes.  Then, presumably, you take the weighted acceleration vector at each frequency point and you get it's Euclidian length, and it's that Euclidian length that you sum up.


Zero padding in frequency domain, any tricks to speedup?

It looks like your choice is based on whether (approx) 5*NLogN (IFFTs+FFTs) is greater or lesser than N^2 (linear convolution).

You could use coarser interpolation methods (cubic, spline, etc.) to upsample the spectra but the error would likely be audible as a clear loss in quality.


How to conclude LTI, causality and BIBO stability of a system represented by a differential equation?

I take your equation:
$\frac{d^2y}{dt^2}-3\frac{dy}{dt}+2y(t)=x(t)$

Laplace transform will be:
$ s^2Y(S) - 3sY(S) +2Y(s) = X(s) $

Now I can find transfer function: $H(s) \triangleq \frac{Y(s)}{X(s)}$

$H(s)=\frac{1}{s^2-3s+2}$

Poles are the solution of denominator = 0 and zeros are the solution of numerator = 0. In this case there are not zeros because numerator is always equal to 1, denominator = 0 is a second order equation that can be easily solved.

$ Poles[H(s)] = \{s=1 ; s=2\} \\ Zeros[H(s)] = \{\} $

There are two real poles.

To determine if the system is BIBO stable you can see if all poles have real part $ < 0 $. In this is example there is at least one pole that has real part  $ \geq 0 $ (both in reality) so this system is UNSTALBE.

To be causal the number of zeros must be $\leq$ of the number of poles in this case we have $ 0 \leq 2 $ that it's true and means that your system is CAUSAL.

Time invariant system is a system that doesn't depends on absolute time but only on relative time from time origin. In a transfer function generally there are not information to determine if it's a time invariant system or not beacuse it should be an external specification of your problem. Btw in DSP you consider systems time invariant almost ever.

Finally this is only a quick overview of zero-pole analysis you should study the argument as well and the prove of the concepts I  used.


Training custom classifiers using trainCascadeObjectDetector

Check out this tutorial on how to use trainCascadeObjectDetector. It describes common pitfalls, and includes examples.

There is also a graphical tool for labeling training images on MATLAB Central, which makes preparing your training data much less painful.

By the way, the link you posted is actually the documentation for vision.CascadeObjectDetector, which uses the classifier produced by trainCascadeObjectDetector. The doc page for trainCascadeObjectDetector is here.


Are Chebyshev Type 1 filters very unstable at very high orders?

With infinite-precision arithmetic they will be stable. However, since you don't have infinite-precision arithmetic, you will have quantization issues even if you use 64-bit precision. These quantization issues can make your filter unstable. Even if your filter is stable, perhaps you will not get the frequency response wanted because of these quantization issues.

First solution : split your order-30 IIR filter into 15 order-2 IIR filters and cascade them. This should fix your problem.

https://www.dsprelated.com/freebooks/filters/Series_Second_Order_Sections.html
https://www.mathworks.com/help/signal/ref/tf2sos.html

Alternatively : Instead of cascading, you can parallelize the order-2 IIR filters.

https://www.dsprelated.com/freebooks/filters/Parallel_First_and_or_Second_Order.html

Finally, do you really need an order-30 chebyshev filter? I find it hard to believe to you really more than an order-8 or order-10 IIR filter.

Edit : You edited your post to mention that you already use biquads. Have you tried an simpler filter, like an order-4 IIR? Make sure it works for an order-4 IIR filter (i.e. 2 cascaded biquads) and then try increasing the order. Also, are you sure that your factorization algorithm is correct? Maybe, there's a simple bug hidden somewhere and it will be easier to spot it with a reduced-order filter.


How does transformation matrix and camera matrix change when indexing pixels from -1 to 1 instead of 0 to N-1

Before attempting to answer this, let's first see what camera transformation matrix does. Consider the 3d point
$$ P_1 = [X_1, Y_1, Z_1]^T $$
Pre-multiplying by camera matrix gives
$$ K \cdot P_1 =
\begin{bmatrix}
\frac{W/2}{tan(hfov/2)} & 0 & W/2 \\
0 & \frac{H/2}{tan(vfov/2)} & H/2 \\
0 & 0 & 1 \\
\end{bmatrix}
\begin{bmatrix}
X_1 \\
Y_1 \\
Z_1
\end{bmatrix}
$$
$$ =
\begin{bmatrix}
\frac{W/2}{tan(hfov/2)} X_1 + W/2 \\
\frac{H/2}{tan(vfov/2)} Y_1 + H/2 \\
Z_1
\end{bmatrix} $$
As we can see, $W/2$ at $K[1,1]$ scales $x$ value and the $W/2$ at $K[1,3]$ adds an offset. Thus, if $x \in [-1,1]$ then by this scaling and offsetting, we get $x \in [0,W]$.
Conversely, when $x \in [0,W]$, we should use the camera matrix as defined earlier. Hence when $x \in [-1,1]$, we shouldn't add the scaling or offsetting. Thus, only the camera matrix needs to be changed to
$$K = \begin{bmatrix}
\frac{1}{tan(hfov/2)} & 0 & 0 \\
0 & \frac{1}{tan(vfov/2)} & 0 \\
0 & 0 & 1 \\
\end{bmatrix}$$


Signal correlation vs Pearson correlation between signals

Thank robert bristow-johnson for the right question.
In my case signal coefficient shows a some level of similarity between noise-like signals because of appreciable mean values. That is as you compare two electrical noises with some direct voltage: of course they are as more similar as high equal direct voltage they have.
Since I need to compare specifically shapes, I should use Pearson method only (which "normalizes" data itself in terms of mean values).


Nonlinear tricks for blocking wide band noise

A very common non-linear filtering technique for removing noise while preserving edges is median filtering. It basically replaces the samples inside a window by their median. Total variation denoising is another non-linear technique that works especially well for denoising piece-wise constant signals. You can find some more information about this method and Matlab software here.