Slow Down Music Playing While Maintaining Frequency

Yeah some of us can do it, you can speed up or slow down without affect the pitch, some guys call this applications of Time Stretch, there different ways to do it, you can do in frequency domain or time domain, you will need choose what is best for you, you will find some advantages and disadvantages of each.

Time Domain:

In Time Domain you can try some techniques like:

TDHS(Time domain Harmonic Sacaling)
SOLA (Synchronous Overlap Add)
PSOLA(Pitch Synchronous Overlap Add)
WSOLA(Waveform Similarity Overlap Add)

Pros: Is fast, some algorithms are easy to understand, good quality in monophonic sounds.

Cons: Generally you'll need a very nice pitch track to splice in the right position, it is hard to do :-(, so if your pitch track fails or not work in Poliphonic sounds this algorithms will give a lot of glitches/artifacts in the output sound.

Frequency Domain:

All time stretch that i know in frequency domain are based in phase vocoder techniques.

Pros: Will work in polyphonic or monophonic sounds.

Cons: Can be painful understand all the math, implementation is a bit hard, is not so fast like time domain codes, for voice i prefer time domain results, some tricks to improve the result of the standard phase vocoder are not shared.

I can say that the window and the hop size are one of the key for the phase vocoder quality, generally we choose 4x overlap to resynthesis, one hann window of size 4096 is enough for my ears (of course if u have processing power for this sizes), the standard phase vocoder can add some reverberant effetcts, to try avoid this kind of problems you maybe need lock the phase.

For datailed take a look in the Miller Puckette and Portnoff paper


Why is histogram equalization called as "equalization"?

To answer your question histogram equalization is called like this because its function is to produce an equalized histogram (that is an uniform probability density function).

There are different algorithms that may approach this function, and obviously there is a problem in the example that is shown:

In fact, the algorithm used there will always have trouble producing a flat histogram when mapping the function to a discrete set of 256 (except when the input histogram is already flat!)
As a consequence, the cumulative histogram shows discrete jumps which produce this non-regular output histogram.

Why is this simple algorithm still used? Because,

the output image will often good enough for the initial purpose,
that is to maximize the range of intensities used in the image - though the histogram may look ugly...
it is simple!

There are simple alternatives. The simplest is to use a linear interpolation of the cumulative function.

Theory: In general terms, the method may be easily derived if we know the probability distribution function $dP_{i}$ of variable $a_i$ representing the luminance image at pixel $i$. Indeed,  by choosing the non-linearity as the cumulative distribution function transforming any observed variable $\bar{a}_i$ into:
$$
z_i(\bar{a}_i)= P_{i}(a_i \leq \bar{a}_i) = \int_{-\infty}^{\bar{a}_i} dP_{i}(a_i)%
$$
This is equivalent to the change of variables which transforms the luminance image to a variable with uniform probability distribution function in $[0, 1]^M$.
Evaluation: In physiology (retina for instance), this distribution is assumed to be sampled for each "pixel" (photoreceptor in mammals) over time (for instance, something like "the probability distribution over last minute"). In image processing, when handling one single image, you can assume stationarity over space, that is compute the histogram over all pixel values.
Algorithm: now that we have all in place, some example algorithm in python:

Sample code:

import numpy as np
from scipy.misc import lena
I = 1. * lena()
I = (I - I.min())/(I.max()-I.min())
import matplotlib.pyplot as plt
# the cumulative histogram
n_quant = 40
(hist, bins) = np.histogram(I.ravel(), bins=np.linspace(0., 1., n_quant+1), normed=True)
P = np.cumsum(hist/np.sum(hist))
# equalization
I_z = np.interp(I, np.linspace(0., 1., n_quant), P)
(hist_z, bins) = np.histogram(I_z.ravel(), bins=np.linspace(0., 1., n_quant+1), normed=True)
plt.bar(np.linspace(0., 1., n_quant), hist, color='r', alpha=.3, width=1./n_quant)
plt.bar(np.linspace(0., 1., n_quant), hist_z, color='b', alpha=.3, width=1./n_quant)
plt.title('Histogram of Lena')
plt.xlabel('luminance')
plt.ylabel('probability')
plt.xlim((0,1.))
plt.savefig('equalization.png')

The output (red is raw histogram, blue is equalized). As you can see, there are variations, but it looks better :-). Notably, you have less "jumps" around most popular values. The parameter n_quant is important to tune smoothness versus precision.


2-D Wavelet output changing when using 'figure' command in Matlab

The reasoning is simple. MATLAB defaults to a 'jetblue' colormap when you use imshow(im).

Typically you need to specify the colormap being used when you use imshow.

For this, add another command like:

$$colormap(gray(256));$$

to your current script.


Why do we care about the multiplicity of poles and zeros in rational Z-transforms?

Multiple poles have amplitude and phase profile different to that of single poles.


How can I change the formula to detect the barcode in every angle? - MATLAB

If the filter you have works well enough for you, you could simply use imrotate to create several rotated versions of the source image, and use the filter you have on each of them (make sure to use bilinear to get sensible derivatives)


For spectral features like roll-off, flux, centroid, etc, the process is as follows:

Before delving into the extraction process, a pre-processing phase is sometimes necessary, this includes:
* Converting stereo signals to mono by summing both channels.
* Resampling the signal to a target sample rate (if a collection of signal of heterogenous sample rates are used).

Now, we can assume that we are dealing with a mono signal of known sample rate.

At each time $t = 0, \Delta, 2\Delta...$ at which the feature has to be calculated (for example, $\Delta = 0.04s$ if 25 values of the feature are to be extracted per second of audio):

Extract the block of signal in the range $[t, t+l]$. $l$ is the analysis window duration and must be greater or equal than $\Delta$. There is a compromise to be found between using short $\Delta$ and $l$ (good temporal resolution) and using a long $l$ (good frequency resolution).
Multiply this signal chunk by a windowing function such as the Hann window. This prevents spectral leakage which might cause unwanted results (for example, report abnormally high values of spectral spread for a purely sinusoidal signal).
Eventually, pad the resulting signal with zeros. The purpose is twofold: make sure the number of samples processed at the next step is a power of 2 (in case the FFT library you use is rustic and does not handle other sizes); and in some applications, artificially boost the frequency resolution through spectral interpolation.
Compute the FFT of the windowed, padded, block of signal.
Discard the upper half of the coefficients returned by the FFT function (they correspond to "negative frequencies" and are of no interest).
Compute either the magnitude or square magnitude of the remaining coefficients. This gives the spectrum on which spectral features are computed.
From there, use whatever formula describes the feature to be extracted. The formulas usually involve $A(k)$ (squared magnitude or magnitude in the $k^{th}$ FFT bin) and $f(k) = \frac{k . \text{sample rate}}{\text{FFT size}}$ (Center frequency of the $k^{th}$ FFT bin).

There is little consensus in the audio signal processing community as to whether the magnitude or power spectrum should be used for spectral features - I have routinely seen both used in research code.


Fitting fixed number of complex exponentials to complex signal

Complex A and B are a bit messy. What about using the signal model

$$ f(t) = \sum_{k=1}^{2} A_k\, e^{j(w_k t + \phi_k)} $$

where $w$ is the frequency and $\phi$ is phase.

Since sinusoids are orthogonal, probably the easiest thing to do is to Fourier transform the signal and pick out the two largest peaks in the positive part of the spectrum. These will correspond to the two main sinusoids.


Since your shift is going to be in the 0-5 pixel range, I am assuming you are going to want a method that gives you sub-sample accuracy.

I have a suggestion for you to try.  Read my article Exponential Smoothing with a Wrinkle and apply the technique to your 1D signal.  Look for some matching zero crossings in the difference case near the center of your data frame.  You should be able to determine the zero crossing by linear interpolation of the two nearest points, or do linear regression on the nearest four or five.  The linear interpolation should be sufficient.  You should be able to measure the shift at several zero crossings. If it is a true shift, all the values you calculate should be the same.  Average them for noise mitigation.

This method is easy to code and takes way fewer calculations than the method you are using.

BTW, that is a very good presentation of your question.

Ced


Frequency response is $1/(1−ce^{iw})$, if you want to get amplitude information you just need to get the magnitude of this complex function and its phase will be the inverse of the phase of $1/(1−ce^{iw})$ or

$$-\arctan (-c\sin w/(1-c\cos w))$$

You can compute the magnitude of the transfer function by:

$$ 1 / sqrt( (-c\sin w)^2 + (1-c\cos w)^2 )$$


This is up to how you model noise and channel.

In the classical, most-commonly-used models, the channel is actually noise-free and the noise is only added at the receiver.

If you have a different channel, that's kind of specific to how your special problem is being modeled, and we don't know that.

But since you say, you're new to communications systems:
The cases you're most likely to encounter are actually that the noise only happens at the receiver.

If that's not the case, look for the context in which you're working. For example, if you're in the context of diversity combining, and are considering a channel that has interferers, e.g. because  you're trying to implement optimum combining, then your overall noise is usually described by a covariance matrix – in which you don't really care about where the noise energy actually comes from (there's the assumption that receiver noise is uncorrelated and hence manifests as diagonal entries in the matrix).


Calculate ifft using only REAL forward FFT

EDIT: This answer is basically a more detailed version of Hilmar's answer, which I noticed too late. Below there are more details about the derivation and the rearrangement of the outputs of the real FFT.

I assume you use the definitions of the FFT and IFFT as commonly used in DSP:

$$\textrm{FFT: }\quad X[k]=\sum_{n=0}^{N-1}x[n]e^{-j2\pi nk/N}\tag{1}$$

$$\textrm{IFFT: }\quad x[n]=\frac{1}{N}\sum_{k=0}^{N-1}X[k]e^{j2\pi nk/N}\tag{2}$$

For a real FFT we assume that $x[n]$ is real-valued. We also assume that the FFT length $N$ is even (otherwise you don't get an FFT bin at Nyquist).

From $(1)$ and $(2)$ we see that

$$\text{IFFT}\big\{X[k]\big\}=\frac{1}{N}\left(\textrm{FFT}\big\{X^*[k]\big\}\right)^*\tag{3}$$

Since your FFT routine only accepts real-valued inputs, we need to write $X[k]$ in terms of its real and imaginary parts:

$$X[k]=X_R[k]+jX_I[k]\tag{4}$$

Furthermore, if we assume that the inverse FFT of $X[k]$ is real-valued, we can rewrite $(3)$ as

$$\begin{align}\text{IFFT}\big\{X[k]\big\}&=\frac{1}{N}\textrm{Re}\Big(\textrm{FFT}\big\{X_R[k]\big\}-j\,\textrm{FFT}\big\{X_I[k]\big\}\Big)\\&=\frac{1}{N}\Big[\textrm{Re}\big(\textrm{FFT}\big\{X_R[k]\big\}\big)+\textrm{Im}\big(\textrm{FFT}\big\{X_I[k]\big\}\big)\Big]\tag{5}\end{align}$$

A little extra inconvenience lies in the fact that a real FFT routine only computes $N/2+1$ values whereas Eq. $(5)$ assumes FFT results of length $N$. This Matlab/Octave script shows how to correctly arrange the results of the individual FFTs:

N = 64;     % must be even!
N2 = N/2;

x = randn(N,1);

X = realfft(x);
XR = real(X); XR = [XR;XR(N2:-1:2)];
XI = imag(X); XI = [XI;-XI(N2:-1:2)];

x2a = real(realfft(XR)) / N; x2b = imag(realfft(XI)) / N;
x2 = [x2a + x2b; x2a(N2:-1:2) - x2b(N2:-1:2)];

In the above script we assume that there's a function realfft.m defined by

function X = realfft(x)

N = length(x);
if rem(N,2), error('N must be even.'); end

X = fft(x,N);
X = X(1:N/2+1);

end


Inverse Z transform of a left sided signal

Note that the order of your numerator equals the denominator, first perform a long division to simplify, before performing partial fraction expansion:

$$ H(z) = \frac{ 1 + 1/6 z^{-1}}{1 - 0.25 z^{-1}} = -\frac{2}{3} + \frac{5/3}{1 - 0.25 z^{-1}} $$

(Note: long divison already simplified the expression that it does not require a further partial fraction expansion to be applied, hence skip step-2 now)

Now apply the inzerve Z-transform and note that for left sided inverse sequence wil be $-a^n u[-n-1]$ for a term like $1/(1-az^{-1})$.


DSSS direct sequence spread spectrum - encoding each raw data bit (question)

In your title and your question overall, you confuse bits and symbols. Not a big deal – we can just work with the idea that we use a constellation where one symbol is equivalent to one bit (e.g. BPSK). But remember that there's no restriction that says you can't DSSS e.g. a 64-QAM modulated symbol stream. (In fact, that's what IEEE802.11b does, iirc).

I'll hence use "symbol" in place of "bit" where I cite you below, whenever it doesn't make a difference. When it does, I'll try to explain.

  Will this mean that we have to find hardware that generates a clock rate (for the PN generator) having 15 times the frequency of the raw symbol rate?

Yes! We call that higher rate chip rate, by the way.

  Initially, I was thinking that a high frequency clock could be used as a reference. However, conventional flip flop frequency dividing methods divide frequencies by even amounts, like divide by 2, divide by 4 etc. There is no option for 'divide by 15' (or divide by an odd number).

Well, you can certainly build a digital counter that counts to 15 and resets – your 1/15 clock divider. So, I think you might be imagining problems where there aren't problems. It's not as trivial as a $2^N$ counter, but honestly, a $2^N-1$ counter is just either a $2^N$ counter where you reset all bits to 0, but the lowest bit to 1 (giving you a single count of headstart), or really just a counter with a $N$-input NOR gate. Or really, just a counter that can compare to an arbitrary value. All these aren't hard to build. This isn't the 1960s!

  So if I want to use LFSR PN sequence patterns for encoding raw symbols, then what kind of method should be used? Is this the usual problem faced when using maximal LFSR generators?

The thing is even easier: You've got a device that repeats itself every 15 states, and changes at the high clock rate. Why not simply use that exact device to build your frequency divider? You could just compare all the state of SR to the one state on which you output a pulse. Done!

But again, integer frequency division really isn't that hard to implement at the rates you're probably worried about.

As another matter of fact: you'd probably not even use the higher frequency as reference, but a lower frequency, from which you synthesize the symbol and the chip rate, and quite possibly also the RF local oscillator. Since the relationship between chip and symbol rate is a rational one, a relatively simple integer-N synthesizer might be used, but frankly, aside from spur performance, there's little reason not to go to fractional-N synthesizers – which, by the way, are especially made for such systems.

  Due to the odd-valued number of bits, is there generally going to be an issue with coming up with a system for generating the clock rate that is 7 times higher than the raw data rate? I can understand that generating 8 times data rate would be relatively easy to do, because 8 is a 2^N value.

There's no hard tasks here. Count to 7, emit a pulse, reset your counter. Counters are relatively easy to build capable of high clock rates.

  Would fractional DDS type frequency generation methods be used instead?

Not quite sure what fractional DDS is, but a DDS method involves jumping through a table at high rates, to produce lower-frequency tones. I don't even think that's easier than a odd-number-counter to build a frequency divider. DDS is generally used to produce tones, not digital clocks.


It isn't. The convolution sum is

$$y[n]=\sum_{k=-\infty}^{\infty}x[k]h[n-k]=\sum_{k=-\infty}^{\infty}\alpha^ku[k]u[n-k]\tag{1}$$

The step $u[k]$ in $(1)$ can be removed if the index $k$ starts at zero, and the step $u[n-k]$ can be removed by replacing the upper summation limit by $n$. For clarity, one should multiply the result by $u[n]$ because if $n<0$ the result of the convolution is zero:

$$y[n]=u[n]\sum_{k=0}^{n}\alpha^k\tag{2}$$


Lighting for computer vision automisation

What I learned so far that my lighting is not good. It's coming from a fluorescent lamp from that's right above the camera. Apparently I would have to diffuse the light with some sort of material. The reflections of the light are not doing us well.

A "Ring Light" type of light will cast uniform light and minimal shadows to your object. It is used regularly in photography. The link to the product is for demonstration, not to suggest that specific product.

Of course, the other thing with the ring light is that the camera shoots through the ring and it's not casting shadows.

Using diffusers is not efficient for relatively large distances (a few meters), or rather, using a diffuser means that you have to use a much stronger light behind the diffuser to flood a large area.

Now, once the light is flooding the area uniformly, place a sample tape rectangle and extract the Red, Green, Blue (RGB) values of the pixels occupied by the tape. To determine your threshold do a histogram on the RGB channels and try to characterise the modes on the R and the B channel.

A more involved way of doing this is to scan the image once and convert the RGB values to Hue Saturation Value (HSV). You can then do one histogram on the H channel which will have two modes, one on the angle that corresponds to Red and one on the angle that corresponds to Blue. Consequently, you now have one variable (H) with two thresholds and two deviations around the thresholds.

Hope this helps.


How this 36*1 came and how we calculated it?

HOG is an algorithm which:

works on a portion of the image, called "detection window";
divides the "detection window" in a certain number of cells;
associates an histogram (of oriented gradients) to each cell.
Each histogram has N orientation bins. WLOG, here, we'll consider N=9);
"Normalization step" (to make histogram values not affected by lighting
variations): it creates blocks (where each block is made up of 4 cells
(2x2)), and it overlaps them;
creates a vector for each block (so, 4 cells x 9 orientation = 36 components for block)
creates a final vector, associated to the overall "detection window";
the final vector is passed to an SVM classifier.

At each iterations, the "detection window" scans the whole original image, and the HOG algorithm is applied at each scan.

  and is it compulsory that we always need 9 bin vector? Is it a fixed size for HOG?

No, the usage of 9 orientation bins for each histogram is not mandatory; it's only a trade-off to which the authors came through tests, as stated by MimSaad.

In this paper, it is explained that increasing the number of orientation bins improves performance significantly up to about 9 bins; beyond this value, perfomance doesn't improve considerably.
Figure 4b shows this result.


Why does Hilbert filter distorts the shape of the signal?

A Hilbert transform is dispersive, as in any constant phase shift over frequency, which will distort a signal that has bandwidth. This is because a fixed DELAY, which would have no distortion, has linear phase versus frequency. If you fix the phase versus frequency to be constant (as we do in a Hilbert transform), the delay at each frequency will be different and thus will cause significant distortion to the signal.

This is demonstrated with my slide below that is used for explaining why linear phase filters are important specifically, but also helps demonstrate what a constant phase shift versus frequency would do. Note that delay specifically is the derivative of phase versus frequency:

$$delay = \frac{d\phi}{d\omega}$$

In the slide below, the first three Fourier components of a square wave signal are shown, along with their addition when they all add up in time (which would happen with a linear phase filter as the delay in time given the formula above would be the same for each frequency). The "Non-linear Phase" case shows what happens if the Fourier components are summed after going through different time delays, which occurs in non-linear phase filters, AND occurs if the composite signal goes through a system where the phase is constant for all frequencies (as the time delay for each frequency must be different in order for the phase to be the same).


This is tricky. Loudness is a very complicated perceptual phenomenon and playback of sound in a room is a very complicated physical phenomenon. To do this right, you would need calibrated sound pressure signal measurements at the listener's ears and then run these through a suitable perceptual model (here is a good overview article: http://www.tcelectronic.com/media/1014016/skovenborg_2004_loudness_m.pdf)

For sine waves you could use the ISO 226 curves but with calibrated signal they are only correct for 1kHz. Let's say compare to sine waves at 200 Hz and one has ten times the amplitude: If playing back the first one 30 dB SPL the second one will be about 4 times louder. If you do this at 80 dB SPL the second will only be 3 times as loud at best. This gets even more complicated when multiple frequencies are involved.

Things get a lot more messy if you include the playback system as well. For pure sine waves the perceived loudness (and energy at the ears) can vary by as much as 30 dB (factor of 8 or so) just by moving a few centimeters. That's a basic effect of room acoustics. Furthermore most playback systems are fairly non flat (either by design or by sloppiness). Most cheap set of computer speakers will not provide a lot of output at 10 kHz.

Your best shot is probably to assume the signal are "normal" (music, speech, something broad band), apply an A-weighting, assume some average playback level and calculate a loudness difference from there. This is a gross simplification but depending on the application this may or may not be good enough.


Typically to compare segmentation algorithms, people use an NPR index. NPR stands for Normalized Probabilistic Rand.

An over segmented or an under-segmented image will typically have negative values of NPR whereas a properly segmented image will have positive values of NPR index.

This type of comparison is typically carried out for segmentation algorithms of the split-and-merge type.

Following this paper and the further publications of the group would be advantageous.


How does integration of a signal change the spectrum?

It is the basic property of Fourier Transforms. Given that $g(t)$ and $G(f)$ are Fourier transform pairs i.e.
$g(t)  \rightleftharpoons G(f)$
then
$\int_{-\infty}^tg(\tau)d\tau \rightleftharpoons \frac{1}{j2\pi f}G(f)$

This assumes that $G(0)=0$. If $G(0)$ is not zeros then the integral of $g(t)$ has a Fourier transform that includes a Dirac delta function.


Plot the frequency response of MTI delay line canceller using Octave

Magnitude |H(z)| is always defined as an absolute value of your transfer function H(z) values which are complex. What you are plotting is the vector that contains complex numbers. You need to take their absolute value and that will give you magnitude.

If you call the angle function then you are going to get phase. Transfer function H(z) is always complex. Try to print values of real(H_z) and imag(H_z), then you will see real and imaginary values of $H(z)$. Now magnitude |H(z)| is sqrt(real * real + imag * imag) for every sample...


Splitting of signal into different components for calculation of PSD

Hint:

$$f_i=8+i\cdot 2\;\text{Hz},\quad i=0,1,\ldots,11$$

Note that $f_0=8\,\text{Hz}$ and $f_{11}=30\,\text{Hz}$. So you have 12 frequencies in the range $[8,30]\,\text{Hz}$, including both edge frequencies.


Estimating the bandwidth of various signals

Another way would be fitting each of your signals' PSD with a Gaussian probability density function whose standard deviation $\sigma$ and mean $\mu$ you estimate, and use the full width at half maximum (FWHM) of this approximation as your half-power bandwidth. This is known accurately with the following closed form:
$$
\mathrm{FWHM}=2\sqrt{2\ln(2)} \sigma\approx 2.355\sigma
$$
This approach is often used in practice to measure things like bandwidth expansion.


Standard deviation of the spectrum of white noise

true white noise has infinite power (because the power spectrum with height $ \frac{\eta}{2}$  has infinite area) which means infinite variance $\sigma^2$ which means infinite standard deviation.  white noise is a conceptual instrument.

to do any analysis of a system with white noise in it, you must somehow determine the (single-sided) bandwidth of the system, $B$, and then the power or variance is $\eta B$.


What channel model to choose while simulating a system in passband with frequencies over GHz range?

However, when I add an AWGN channel, I get pretty low BER of 10^-5 in passband as compared to around 10^-4 BER in baseband under same AWGN channel of SNR 10 dB. Why is it happening?

You've got a bug in either passband channel or baseband channel model, or in how you transform your noise or your signal – end of story!

  I thought that there would be more attenuation in passband as attenuation increases with respect to frequency

If your passband model does that, then your baseband model has to do the same – they are mathematically equivalent.

If your passband model is frequency-selective, then your baseband model has to do the same.


Butterworth High Pass Filter Amplitude Response

If $H_{LP}(j\omega)$ is the frequency response of a Butterworth low-pass filter, then its inverse would go to infinity for $\omega\rightarrow\infty$ (because $H_{LP}(j\omega)$ goes to zero for high frequencies). So that's not what you want.

The correct high-pass frequency response is obtained by replacing $\omega/\omega_c$ by $-\omega_c/\omega$, resulting in

$$|H_{HP}(j\omega)|=\frac{\left|\frac{\omega}{\omega_c}\right|^n}{\sqrt{1+\left(\frac{\omega}{\omega_c}\right)^{2n}}}$$

which approaches unity with increasing frequency, as should be the case for a high-pass filter.


I would suggest to apply a median filter like

y=filter(1,0.05*ones(1,20),x)

And then apply Matlab function "findpeaks"

As an alternative to second step you perform a segmentation based on different from zero areas like in

ne0 = find(A~=0);                                   % Nonzero Elements
ix0 = unique([ne0(1) ne0(diff([0 ne0])>1)]);        % Segment Start Indices
ix1 = ne0([find(diff([0 ne0])>1)-1 length(ne0)]);   % Segment End Indices
for k1 = 1:length(ix0)
    section{k1} = A(ix0(k1):ix1(k1));
end
celldisp(section)                                   % Display Results

And then perform calculation on each segment. For example sum all samples and divide by number of samples.


Convolution theorem with a kernel smaller than the image

If you want to do the convolution in the frequency domain, you have to pad the kernel with zeros to make it the same size as the image.  There is really no way around it.  The upside is that you only need to take the FFT of the filter once, and then you can re-use it for multiple images.

If you are concerned with speed, you would have to do some analysis here. For small separable kernels doing the convolution in the spatial domain may be faster.


Conversion between DCT to DFT frequency domains

Both are linear transforms. Both can be expressed as a matrix. These matrices  are invertible. So you can have the product of one matrix with the inverse of the other. In practice, you could not actually calculate inverse matrices, so you'd have to stick with solvers. However, such products are a little op-art. Coursey GNU octave.


This is more likely a question of what the data reporting application is trying to achieve, rather than a DSP question.

Yes, there is a different term for Gate or Window. In fact, if you "subsample" and low pass filter, as per your first intuition, then you are effectively doing decimation. For more information please see: http://en.wikipedia.org/wiki/Decimation_%28signal_processing%29

The rest of your suggestions would be equivalent to subsampling followed by non-linear filtering such as median, min, max which would have the corresponding effect on the reported data.

For example, a sliding window median filter would "absorb" any spikes in the reported data and the min and max filters would tend to report the highest or lowest values. For more information please see: http://en.wikipedia.org/wiki/Nonlinear_filter and http://en.wikipedia.org/wiki/Median_filter

So, it is really a matter of what does that GUI element try to convey to its user.

Hope this helps.

Supplementary Information:
Yes decimation does imply a stable sampling frequency but most DSP operations can be defined over irregular domains as well.

Irrespectively of that however, I think that we are indirectly looking at the real problem here which is one of visualisation.

The real problem, or at least my understanding of it so far, is that you are trying to convey a comprehensive view of a window of the history of a signal at a rate lower than its actual rate of change. Furthermore, this signal is not used in any further computation. If that is mostly accurate and assuming that your signal is coming from a stochastic process then no single value will be able to represent this history reliably.

Therefore, I would suggest two alternatives:

To use a single boxplot that would represent your last signal's frame
history in a comprehensive way including its min, max, median and
quartile information. For more information please see:
http://en.wikipedia.org/wiki/Box_plot
To convey the signal's history to your user via a Sparkline which
could furthermore be annotated with min,max,median lines. Please note that Sparklines can also be created via the use of special fonts, which makes them really easy to quickly insert to a GUI. For more
information please see: http://en.wikipedia.org/wiki/Sparkline and https://www.fontfont.com/how-to-use-ff-chartwell or http://nsfmc.github.io/chartjunk/

Hope this helps


Finding impulse response of an RC circuit from its step response

The answer (How to calculate the impulse response of an RC circuit using time-domain method)  provides a direct time-domain solution of an RC circuit for the impulse reponse $h(t)$. Now this new answer modifies it to solve for the step-response $s(t)$ instead and then computes the impulse response according to :

$$h(t) = s(t)' $$

The differential equation of the first order circuit was derived as
$$ \boxed{ y'(t) +\frac{1}{RC} y = x'(t) } \tag{1} $$

The step reponse $s(t)$ is defined as the output $y(t)$ of Eq.(1) when the input $x(t)$ is a unit-step function $$ x(t) = u(t) \implies y(t) = s(t)  $$

Let's apply a one stage direct solution to obtain $s(t)$.

The homogeneous solution is found from
$$ y'(t) +\frac{1}{RC} y = 0  \tag{2} $$

The characteristic equation : $s + \frac{1}{RC} = 0 \implies s = - \frac{1}{RC}$.

The (causal) homogeneous solution is :

$$y_h(t) = K e^{-t/RC} u(t) \tag{3}$$

Then, the particular solution $y_p(t)$ will be from the method of undetermined coefficients as follows:

For the particular input $x(t) = u(t)$ we may assume a particular solution as $y_p(t) = A u(t) + B \delta(t)$, then plug this assumed solution into Eq.(1) to find out the coefficients $A$ and $B$. This yields $A=0$ and $B=0$, hence the particular solution is found to be identically zero.

Then since the total solution is:
$$s(t) = y_h(t) + y_p(t) = y_h(t)$$

to find the value of unknown $K$ in Eq.(4) we need an initial condition on the output $y(0^+) =s(0^+) $ ; found from the circuit physics: apply KVL around the loop of voltage source, capacitor and resistor, with the fact that the capacitor voltage is fixed at zero at time $t=0$. Then we find $y(0^+) = s(0^+) = 1$ which yields $K=1$

Hence the step-response is:
$$ \boxed{ s(t) = e^{-t/RC} u(t) } \tag{4}$$

Then the impulse response of Eq.(1) is found to be
$$ h(t) =  s(t)' = \left( e^{-t/RC} u(t) \right)' $$

$$ \boxed{ h(t) = \delta(t) - \frac{1}{RC} e^{-t/RC} u(t) } \tag{5}$$


Use MATLAB to Restore a Signal from a Given Degraded Signal Using Tikhonov Regularization

The idea is to represent all operation sing Matrices.
Once it is done, it is easy to solve the problems as a Least Squares problems.

The way to represent Convolution Operation using a Matrix is by Toeplitz Matrix.
For 1D it is pretty straight forward to do (Just pay attention to boundary).

So let's take the simple model in the comment:

$$ g = f \circ h + n \rightarrow g = H f + n $$

Now, the LS solution is easy:

$$ \arg \min_{f} {\left \| H f - g \right \|}^{2}_{2} = { \left( {H}^{T} H \right) }^{-1} {H}^{T} g $$

The matrix form for the regularized model is:

$$ \arg \min_{f} {\left \| H f - g \right \|}^{2}_{2} + \lambda {\left \| \Sigma f \right \|}^{2}_{2} $$

Where $ H $ is the matrix Toeplitz Matrix generated from the Blur Kernel and $ \Sigma $ is the Toeplitz matrix generated form of the derivative operator (In this case, keeps the signal smooth to regulate the High Pass property of the first term).

This is actually Tikhonov Regularization and it is easily solved.

$$ \arg \min_{f} {\left \| H f - g \right \|}^{2}_{2} + \lambda {\left \| \Sigma f \right \|}^{2}_{2}  = { \left( {H}^{T} H + \lambda {\Sigma}^{T} \Sigma \right) }^{-1} {H}^{T} g $$

All the above is pretty straight forward to calculate in MATLAB once the the $ H $ and $ \Sigma $ matrices are formed (Given the Convolution Kernel, the function convmtx can assist you generate the matrix).

The way to create those matrices is pretty simple.
Given a signal of the length $ N $ (Samples) implies those matrices are $ N \times N $ (Assuming the output has the same number of samples).
Where the $ i $ row of the matrix applies the linear combination for the $ i $ output sample.

For instance, given the Gradient Kernel and a signal of length 5 the matrix becomes:

$$ \begin{bmatrix}
1 & -1 & 0 & 0 & 0 \\
0 & 1 & -1 & 0 & 0 \\
0 & 0 & 1 & -1 & 0 \\
0 & 0 & 0 & 1 & -1 \\
0 & 0 & 0 & 0 & 1
\end{bmatrix} $$

One must pay attention to the rows on the beginning and the end to apply the required boundary conditions.
In the case above I chose zero boundary conditions (Last row).

One should note, that in practice, the convolution matrix is not feasible due to memory limitations.
Hence it must be represented by a Sparse Matrix (Since it is mostly zeros) or the problem should be solved using a solver which doesn't require the matrix explicitly (Like Conjugate Gradient).

Another note would be that if Circular Convolution is after (Imitating DFT Based Convolution) the matrices should not only be Toeplitz but Circulant Matrices which means they enforce periodic boundary conditions (See the circular property in MATLAB padarray function).

Here is a MATLAB Code for the exercise.
I only used 800 samples due to memory constraints.

% Audio Deconvolution
% See http://dsp.stackexchange.com/questions/26433

%% General Parameters and Initialization

% The follwoing code must be executed before any of the following cells.
% The cells must be executed one after one in the order they are written.
clear();
close('all');

set(0, 'DefaultFigureWindowStyle', 'docked');
defaultLooseInset = get(0, 'DefaultAxesLooseInset');
set(0, 'DefaultAxesLooseInset', [0.05, 0.05, 0.05, 0.05]);

titleFontSize   = 14;
axisFotnSize    = 12;
stringFontSize  = 12;

thinLineWidth   = 2;
normalLineWidth = 3;
thickLineWidth  = 4;

smallSizeData   = 36;
mediumSizeData  = 48;
bigSizeData     = 60;

randomNumberStream = RandStream('mlfg6331_64', 'NormalTransform', 'Ziggurat');
subStreamNumber = 57162;
set(randomNumberStream, 'Substream', subStreamNumber);
RandStream.setGlobalStream(randomNumberStream);

addpath(genpath('RawData'));
addpath(genpath('AuxiliaryFunctions'));

%% Loading Data

load('input.mat');
load('handel.mat');

vInputSignal    = input;
vRefSignal      = handel;

vBlurKernel     = [0.0545; 0.2442; 0.4026; 0.2442; 0.0545];
vGradientKernel = [-1; 1]; % Forward Finite Differences
% vGradientKernel = [1; -1]; % Backward Finite Differences

%% Parameters
% Limiting the Number of Samples due to Memory Constraints
numSamplesToProcess = 800;

numSamples = size(vInputSignal, 1);

numSamples = min([numSamplesToProcess, numSamples]);

vInputSignal    = vInputSignal(1:numSamples);
vRefSignal      = vRefSignal(1:numSamples);

% Creating the Operators in Matrix Form
mBlurKernel     = convmtx(vBlurKernel, numSamples);
mGradientKernel = convmtx(vGradientKernel, numSamples);

% Enforcing size to match Convolution using 'same' property
mBlurKernel     = mBlurKernel(1:numSamples, :);
mGradientKernel = mGradientKernel(1:numSamples, :);

% Enforcing Circulant Matrix (Like DFT Based Convolution)
mBlurKernel(1, (end - 3):end) = vBlurKernel(1:4);
mBlurKernel(2, (end - 2):end) = vBlurKernel(1:3);
mBlurKernel(3, (end - 1):end) = vBlurKernel(1:2);
mBlurKernel(4, (end - 0):end) = vBlurKernel(1:1);

mGradientKernel(1, end:end) = vGradientKernel(1);

regFactor = 0.02; % Regularization Factor - Lambda

%% Least Squares Solution

vLsRestoredSignal = ((mBlurKernel.' * mBlurKernel) \ mBlurKernel.') * vInputSignal;

%% Regularized (Tikhonov) Least Squares Solution

vRegLsRestoredSignal = (((mBlurKernel.' * mBlurKernel) + (regFactor * (mGradientKernel.' * mGradientKernel)))\ mBlurKernel.') * vInputSignal;

%% Display Results

hFigure = figure();
hAxes   = axes();
hLineSeries = plot(1:numSamples, [vRefSignal, vInputSignal]);
set(hLineSeries, 'LineWidth', normalLineWidth);
set(get(hAxes, 'Title'), 'String', ['Reference and Input Signal'], ...
    'FontSize', titleFontSize);
set(get(hAxes, 'XLabel'), 'String', 'Samples Index', ...
    'FontSize', axisFotnSize);
set(get(hAxes, 'YLabel'), 'String', 'Samples Value', ...
    'FontSize', axisFotnSize);
hLegend = legend({['Reference Signal'], ['Input Signal']});

hFigure = figure();
hAxes   = axes();
hLineSeries = plot(1:numSamples, [vRefSignal, vLsRestoredSignal, vRegLsRestoredSignal]);
set(hLineSeries, 'LineWidth', normalLineWidth);
set(get(hAxes, 'Title'), 'String', ['Restored Signals and the Reference Signal'], ...
    'FontSize', titleFontSize);
set(get(hAxes, 'XLabel'), 'String', 'Samples Index', ...
    'FontSize', axisFotnSize);
set(get(hAxes, 'YLabel'), 'String', 'Samples Value', ...
    'FontSize', axisFotnSize);
hLegend = legend({['Reference Signal'], ['Least Squares Restored Signal'], ['Regularized Least Squares Restored Signal']});

hFigure = figure();
hAxes   = axes();
hLineSeries = plot(1:numSamples, [vInputSignal, vLsRestoredSignal, vRegLsRestoredSignal]);
set(hLineSeries, 'LineWidth', normalLineWidth);
set(get(hAxes, 'Title'), 'String', ['Restored Signals and the Input Signal'], ...
    'FontSize', titleFontSize);
set(get(hAxes, 'XLabel'), 'String', 'Samples Index', ...
    'FontSize', axisFotnSize);
set(get(hAxes, 'YLabel'), 'String', 'Samples Value', ...
    'FontSize', axisFotnSize);
hLegend = legend({['Input Signal'], ['Least Squares Restored Signal'], ['Regularized Least Squares Restored Signal']});

%% Restore Defaults
set(0, 'DefaultFigureWindowStyle', 'normal');
set(0, 'DefaultAxesLooseInset', defaultLooseInset);

The code is available at GitHub - GitHub StackExchangeCodes.


Matlab's fft functions are all based on FFTW (this is confirmed here), so I guess the obvious choice for you should be FFTW. FFTW is hardware-independent but it can take advantage of some hardware-specific features.


Fourier Transformed a signal $f(t) = \cos(2\pi*4t)*5\cos(2\pi*0.4t)$ but did not get expected results. What went wrong?

Fourier Transform is linear transform. It is mean that Fourier Transforms of sum of functions will be sum of Fourier Transforms of each function. It is not valid for multiply of function.

Your example. There is very well known formula (see 7c in this link):
$$
\cos(\alpha) * cos(\beta) = 1/2 (\cos(\alpha+\beta) + \cos(\alpha-\beta))
$$
use this formula  to transform your function:
$$
f(t) = \cos(2\pi*4t)*5\cos(2\pi*0.4t) = 5/2*(\cos(2\pi*4t - 2\pi*0.4t) + \cos(2\pi*4t - 2\pi*0.4t)) = 5/2*(\cos(2\pi*3.6t) + \cos(2\pi*4.4t))
$$
So Matlab results is correct - you have 2 peaks equal amplitude with 3.6 and 4.4 frequency.


What is databending and how does it actually work?

It's pretty simple, really.

Basically, this relies on the fact that the images must be saved in an uncompressed file.

If you think about how an image is stored, it's functionally a bunch of sequential numbers:

(PX 0,0 - R:x, G:y, B:z)
(PX 1,0 - R:x, G:y, B:z)
(PX 2,0 - R:x, G:y, B:z)
(PX 3,0 - R:x, G:y, B:z)
....
(PX 639,0 - R:x, G:y, B:z)
(PX 0,1 - R:x, G:y, B:z)
(PX 1,1 - R:x, G:y, B:z)
....
(PX 639,479 - R:x, G:y, B:z)

For most bitmaps, R, G, B are either 8 or 16 bit values (1 or 2 bytes). Therefore, the "image" is really just a very long sequence of numbers, with some agreed-upon way of mapping those numbers into pixel-positions and colors.

Now, how is an uncompressed audio file represented? Very, very similarly:

(Sample 0 - R:x, L:y)
(Sample 1 - R:x, L:y)
(Sample 2 - R:x, L:y)
(Sample 3 - R:x, L:y)
....
(Sample z-1 - R:x, L:y)
(Sample z - R:x, L:y)

Most audio samples are 16 bit (some are 24, it can vary) per sample. However, again, the audio file is really just a whole bunch of numbers, with a agreed-upon method for translating the sequences of numbers back into an audio stream.

Now that you understand how the files are structured, you can start to see how the audio effects relate to the images. First, since uncompressed images are generally stored in a raster-scanned topology, audio effects are primarily going to have a horizontal "smearing" effect.

Additionally, you can see that basically the effects are going to permute the images based on the effective horizontal frequency components of each horizontal row of pixels.

The "Random EQ" setting boosts some the frequency components of the horizontal edges in the image, and the rippling that surrounds each edge is a function of the impulse-response of the filter.

The "Distortion" image looks to be adding some high-frequency noise, in addition to some apparent posterization, which makes me think it may decimate the internal sample-values in some manner.

The "Bass-Boot" is clearly exaggerating the swings between wide areas of low-values (e.g. black) and high-values (e.g. light colours).

I'm not sure what exactly the phazer is doing (I don't know what it would sound like, so I can't extrapolate to the DSP it involves).


This is a classic dynamics problem. I think you do not have to go through the convolution to solve these type of problem. Instead, you can draw the frequency response function (FRF) of the system and multiply with the FRF of the forcing function to obtain the resulting FRF. Then the resulting FRF can be transformed into the time domain using the inverse Fourier transform. You can also solve it by doing the convolution. Then you have to use the Duhamel's integral to solve this type of problem.

The link is attached.

http://people.duke.edu/~hpgavin/cee541/sdof-dyn.pdf


How exactly does sound get encoded by a cochlear implant processor?

There's an excellent non-specialist (science for the public) book on the history of cochlear implant technology:

"I Can Hear You Whisper: An Intimate Journey through the Science of Sound and Language" by Lydia Denworth

which provides a non-technical description of the use of channels in the encoding.


Is it possible to "equalise" a signal by deconvolving the impulse response of the room in which it is to be played?

Yes, the answer is here :
Compensating Loudspeaker frequency response in an audio signal ,
The filter here is your room transfer function.


Is the Sobel filter a high pass filter, and if not, what is the difference between them?

They are both highpass type filters, but used with very different intentions.

One should immediately observe the fundamental difference that the output of unsharp masking filter is an enhanced image to be viewed by humans, whereas the output of the Sobel (edge detector) filter is not an image to be viewed by humans, but rather a description of the edges to be further analysed by a machine vision system.

Digital imitation of an unsharp masking effect does not necessarily explicitly use a highpass filter in its implementation, yet effectively the operation that it performs is similar to the operation that a highpass filter performs, therefore it's equally valid to call itself as a highpass filter too.

More specifically, unsharp masking operates by subtracting a lowpass filtered version of the image from itself. Stated in impulse responses;

$g[n_1,n_2]$ being the output image and $f[n_1,n_2]$ being the input:

$$ g[n_1,n_2] = a \cdot f(n_1,n_2) - b \cdot f_L(n_1,n_2) $$

where $ a > b $ are two constants and $f_L[n_1,n_2]$ is the low pass filterd image,

$$f_L(n_1, n_2) = f(n_1,n_2) \star h_{LPF}[n_1,n_2] $$

where $h_{LPF}[n_1,n_2]$ is a lowpass filter and $\star$ denotes convolution.

Reordering the equation and noting that $$f[n_1,n_2] = f_L[n_1,n_2] + f_H[n_1,n_2]$$ we can see that

$$g[n_1,n_2] = (a-b) f_L[n_1,n_2] + a f_H[n_1,n_2] $$

or more specifically

$$g[n_1,n_2] = f[n_1,n_2] \star ( \delta[n_1,n_2] - a h_L[n_1,n_2])$$
$$g[n_1,n_2] = f[n_1,n_2] \star h_{eH}[n_1,n_2] $$

The last variable $h_{eH}[n_1,n_2]$ is the effective highpass filter that the unsharp mask algorithm is performing. So as you can see, it effectively includes a highass filter.

Unsharp masking is a picture enhancement technique borrowed from analog film processing era, which is used to improve the apparent sharpness of an image to be viewed. Since it increases the sharpness, it's effectively a high pass filter, but not very strong. It enhances edges, by making them stronger, too, and this is how it achieves improved, percepted, sharpness effect.

Sobel filter is actually an edge detection method and as such it's a nonlinear operator. However, its internal operation involves (basically) two directional impulse responses along horizontal and vertical directions that apporixmate the gradient of the image by using two LTI high pass filters. As, such it's using high pass filters internally, but its output image is obtained by further nonlinear processing (to keep edge pixels and throw out the rest) whose frequency response character does not exist, conventionally.

Sobel filter is, effectively, a very strong highpass filter used for machine vision purposes; specifically for edge detection. Since it's a strong highpass filter, it will eliminate all the low frequency content and its output will be black background and white lines denoting the edges, in addition to noisy dots and contours. It's never used to improve the sharpness of an image.

Sobel detector is fundamentally a derivate filter; note that derivative filters are always of highpass character, whereas integration filters are of lowpass character. They (the Sobel filters) compute the approximate derivative (or Gradient in 2D application) of a function at around the pixel points. This derivative is an indicator of a possible edge.


Calculating spatial width from power spectrum HWHM

Say I have a 1D (spatial) signal (resolution = $1000$) which is zero everywhere except from $x=250$ to $750$, where it equals one.

This is not "resolution". Resolution is 300 Dots Per Inch. In which case, we could say that the total physical length of your pulse is $\frac{500}{300} \approx 1.666$ inches (or any other unit of length).

  I ultimately want to calculate the spatial width of this signal using FFTs.

Why?

  I am dealing with a signal that evolves with time and wish to calculate the average "pulse" width over all the time frames, so I do not know the widths.

If there will be multiple pulses of different widths on the same signal, then by opting to detect them with the FFT you are setting yourself a very big challenge because the FFT would return to you information about the signal as a whole. So you could, for example, derive an average rate of pulses (even using the algorithm that is presented here) but not the widths of individual pulses.

If it is somehow guaranteed that within a window of 1000 samples, there will be a pulse whose length is guaranteed to be staying well below 1000 samples and all we have to do now is detect where the pulse is and how long it is, then opting for the FFT is an overkill.

The usual way to detect pulse widths is via the simple use of a threshold and a counter. Once the signal's amplitude goes above the threshold, the counter starts counting and it stops once the signal's amplitude goes below the threshold. If you are going to operate in a noisy environment, then there are a number of improvements to that such as adding hysterisis to the threshold, so that it doesn't respond to very short "bounces" of the waveform and adaptive thresholding where the threshold limit would be derived from the given window of observation (here, from the 1000 samples).

If you absolutely have to work in the frequency domain, it might be better to look into the discrete wavelet transform (DWT), whose output is a time/scale(frequency) representation. But the actual detection of the pulse width is likely to be happening (again) using some form of threhsolding on the output of the DWT. (So, again, huge overkill).


okay, so 9 hours ago i posted in a comment links to other questions/answers regarding this topic for you to maybe pre-digest.

this is the issue of audio-to-MIDI conversion, which has been discussed, particularly among the Music Information Retrieval(MIR) people.  it's not a perfectly solved problem, especially for anything polyphonic.  even picking apart a duet is hard.

for monophonic, it's not so hard to get a decent pitch contour plot (that is the detected pitch vs. time) along with a plot for amplitude vs. time and a plot for pitch confidence or tone harmonicity (or the degree of periodicity of the note waveform) vs. time.  there are other measures that one can do, such as a few different timbre parameters (like brightness, even/odd harmonic energy ratio).

what would be interesting maybe for you to do is, with your favorite monophonic test file (people often use Tom's Diner from Suzanne Vega) and see what pitch contour you get from YIN.  and compare that to the pitch contour you get from using the ASDF (turned into autocorrelation) and peak picking from my other answers:

Fast pitch recognition

How to deal with low fundamental when using AMDF for pitch extraction?

What is an AMDF?

What is the difference between pitch detection and Onset detection?

so first make sure you get a decent pitch contour for your monophonic voice and see if the "octave errors" get patched up before looking for note onsets.

for many notes, the onset is clear by looking at sharp transitions in the amplitude contour curve and in the "harmonicity" curve.  but if you're humming and not having any impulsive or sharp beginnings to notes, these curves make it hard to detect note onsets.  sometimes then you need to look at somewhat sharper pitch transitions and also when the pitch crosses semi-tone boundaries.  but somehow you don't want vibrato to trigger new notes if the vibrato happens to bring your note a little closer to another pitch in the 12 notes per octave scale.

finally, once you do effectively segment your notes, then you need to look at the pitch contour and decide, closer to the note onset, which of the equally-tempered, $\frac1{12}$ octave pitches you are closest to.  but that can have a pitch offset if the whole song is a little sharp or a little flat (but the notes are well tuned to each other).


Can someone explain why absolute is needed in power or energy signal?

Power and energy are real-valued quantities, so for a complex-valued signal $x(t)$ you need to use $|x(t)|^2$. For real-valued signals, you don't need the magnitude because $x^2(t)=|x(t)|^2$ holds. Note that this equality is not satisfied for general complex-valued signals.


You can read about shaping gain in any book about digital communication. For example, John G. Proakis, "Digital Communications".
The theory is presented in J. H. Conway and N. J. A. Sloane, "Sphere Packings, Lattices and Groups".

In SCMA case, the multidimensional mother constellation only potential can increase shaping gain. There are many different methods for design good SCMA codebooks. I conducted a study comparing known codebooks and presented method for design codebooks based on genetic algorithm. For AWGN channel, there are codebooks based on QPSK alphabet with very good performance (in comparison with codebooks based on good multidimensional constellation).

In addition to a good multidimensional mother constellation, it is necessary to optimize the rotation factors, power variation coefficients, permutation, etc... Only a good multidimensional mother constellation (it has a shaping gain) does not guarantee the good performance of SCMA. So, QPSK mother constellation without shaping with good optimization of the rotation factors may provide the good performance of SCMA system.


Bandpass filtering an audio signal with MATLAB's $\tt kaiserord$: unknown new frequency component

Your filter is a stop band attenuation of about 40 dB. Frequencies outside your passband will be attenuated by at least 40 dB but they will not go down all the way to zero.

The FIR filter is very long (2000 or samples). It takes that many samples for the filter to reach "steady state". In other words, you should throw away the first 2000 samples of your filtered signal before doing any analysis on it.


Can the (point-wise) multiplication of two discrete band-limited functions be inherently aliased?

(1) is it possible that the point-wise multiplication of two discrete, band-limited functions is aliased itself?

Yes.  Or at least it's possible that the result won't fit in a tidy manner into your sampling rate.  Because multiplying two signals together creates energy at different frequencies than the frequencies of the original (this is easiest to do by a single-frequency example: $2 \cos \omega_1 t \cdot \cos \omega_2 t = \cos (\omega_1 - \omega_2)t + \cos(\omega_1 + \omega_2)t$.  If that $\omega_1 + \omega_2$ is to great to fit into your sampling rate, then it'll alias.

  (2) would it be possible to circumvent the aliasing by outputing the 'full' convolution, thereby getting $2−1$ samples back?

Yes and no.  Yes, you could, but to be consistent the result would have to be at a sample rate $\frac{2M - 1}{M}$ times the original.

An equivalent operation would be to upsample both signals by a factor of two and then multiply -- if the upsampling were done in a way that results in correctly band limited signals, once you multiplied them they would fit into the new, bigger bandwidth made available to them by upsampling.


As pointed out by Dan Boschen in a comment, these number aren't really surprising. Usually the DFT is defined as

$$X[k]=\sum_{n=0}^{N-1}x[n]e^{-j2\pi kn/N}\tag{1}$$

where $x[n]$ is the input sequence and $N$ is the DFT length. So bin zero, i.e., $X[0]$, is just the average value of your input sequence multiplied by the DFT length. The maximum value of $|X[k]|$ can be as large as the maximum value of $|x[n]|$ multiplied by the DFT length.


How does it go about implementing pitch shifting without affecting the duration or speed on matlab?

For the speed change: Sample rate is the rate at which samples are read right? So, as per that definition, since you have a higher sample rate, you would be reading more samples in one second. The length of the file remains the same but now since the samples are read faster, the overall length of the file would become slower (with an increase in sampling frequency).

For the pitch change: Now, let us consider the case of increasing the sampling frequency. I'm sorry for giving a rough explanation for this but imagine increasing the sampling frequency would squish the signal in the time axis since now lesser samples represent the same length in samples. This would mean an increase in frequency, thus resulting in a pitch shift.


Upsample data using FFTs. How is this exactly done?

One trick, for even-length signals, is what to do with the "middle" sample.  Here, I've split it half and half between each side of the FFT.

The other trick is to ensure that you have the right amplitudes in the resampled signal. Here's it's a factor of 2.

Try this in scilab:

x = rand(1,100,'normal');
X = fft(x);
XX = 2*[X(1:50) X(51)/2 zeros(1,99) X(51)/2 X(52:100)];
xx = ifft(XX);

clf;
plot([0:199]/200,xx)
plot([0:199]/200,xx,'go')
plot([0:99]/100,x,'r.')

which appears to generate the right thing. The red dots are the original samples, and the green circles (and blue connecting lines) are the resampled data.

UPDATE

For the case of an odd number of samples, there are fewer fiddly bits (no splitting the last coefficient):

x2 = rand(1,101,'normal');
X2 = fft(x2);
XX2 = 2*[X2(1:51) zeros(1,101)  X2(52:101)];
xx2 = ifft(XX2);

clf;
plot([0:201]/202,xx2)
plot([0:201]/202,xx2,'gx')
plot([0:100]/101,x2,'r.')


How to select the center frequency and filter type for 5-band digital audio equalizer?

The answer may depend on your application. I suggest you look at some commercial audio equipments with EQs to decide for those band frequencies. For example the old analog philips casette-type stereo system in my room has 5 band graph equalizers with the frequencies: [100-300-1000-4000-10000] Hz respectively. I use it as an amplifier for my laptop :)

If you design the first and last bands as Bass/treble shelve types, be careful to limit their cutoffs so that mid bands will have enough spectral intervals.


Fractional powers of complex numbers (DSPrelated computation)

That's my article, I'm glad you are reading it.

A fundamental rule for two complex numbers (corresponding to a point on a plane) is that they are only equal if both the real and imaginary parts (corresponding to x,y coordinates) are equal.  Therefore the complex equation

$$ (a^2 - b^2) + 2abi = i = 0 + i $$

is equivalent to saying the point $(a^2 - b^2,2ab)$ is the same point as $(0, 1)$

The real parts must be equal AND the imaginary parts must be equal.

Graphically, you have the intersection of a hyperbola intersecting crossed lines going through the origin.  Substitute $x$ and $y$ for $a$ and $b$ if that helps you see it.

The important concept is that an intersection point is halfway along the unit circle, and the other one is directly opposite (the negative of it.)  This is true for taking the square root of any number on the unit circle.  One square root will be halfway and the other the negative of it.

And that is part of the exponential nature of the unit circle.  Radians is just one possible scale.  Quarter turns is another, whole turns another.

Expanding on this, the solutions to $1^{1/N}$ are $N$ evenly spaced points around the circle.  They are called the Roots of Unity.  You may know the first root better as $e^{i2\pi/N}$ which is right there in the DFT definition:

$$ X[k] = \sum_{n=0}^{N-1} x[n] \left(e^{i2\pi/N}\right)^{-nk} $$

$$ a^2 - b^2 = 0 $$

$$ a^2 = b^2 $$

$$ a = \pm b $$

And

$$ 2ab = 1 $$

For $a=b$

$$ 2a^2 = 1 $$

$$ a = \pm \frac{\sqrt{2}}{2} = b $$

For $a=-b$

$$ -2a^2 = 1 $$

No real solutions.


Sampling a low-frequency sound doesn't work at some sampling rates

Could be a few things

You may be borderline clipping. Your sine wave has an amplitude of 1, which is just at the edge of clipping (depending on how its rendered). Try it with an amplitude of 0.5
Your hardware is sloppy. For example, cheap laptop sound cards often cut corners in the anti aliasing filters and or clipping management
Your operating system is poorly configured, particular if you use Windows. If you change sample rate in the project, it may just trigger the operating system  to do sample rate conversion. If the sound card only supports 48 kHz natively, the OS will up-convert 44.1 kHz to 48 kHz and the quality of the sample rate converter can be fairly bad. On Windows all audio goes through the kernel mixer (unless you specifically bypass it) and that often messes with the audio as well.
Driver problem. The interface between the OS and the sound card could be sloppy.


Adding measuruments noises to kalman filter

For the first question, in Kalman filter there are two noises namely process and observation noises. Kalman filter assumes the noise to be Gaussian with zero mean and some variance $\mathcal{N}(0, \sigma^{2})$. Let's say we have a sensor that measures a constant quantity but random, for example distance. Let's say $x = 1$m, however the sensor gives us $z=0.9$m. In reality we don't know $x$ for certain so we rely on our sensor. In Matlab, we can simulate this sensor

x = 1; % true value
z = x + randn(); % observed value

randn() is normal distribution with zero mean and variance one $\mathcal{N}(0, 1)$. If you want to increase the accuracy of your sensor, then multiply the normal distribution randn() with some small value that you think it is the error of your sensor.

for i = 1:6
    x = 1; % true value
    z = x + 0.01*randn() % observed value
end

You will notice the observed value varies from the true value however still close to the true value because of 0.01. Remember that Kalman filter is a matter of modeling your system accurately. For kalman filter, you need the motion and observation models given in advanced. Also, your system must be linear and if this is not the case then you should use Extended kalman filter. Moreover, the noise must be Gaussian noise. Now your job to determine these parameters.


Design of Parks-McClellan Lowpass Filter

The minimal specification for a FIR lowpass filter designed by the Parks McClellan algorithm is the passband edge, the stopband edge, and the filter order. You don't need to specify passband and stopband ripples; they are the result of the design process. What you can do is specify passband and stopband error weights, which let you specify the relation between passband and stopband approximation errors (a larger weight reduces the error).


It is not a good idea to have Np and Ns as constraints, for the exact reason you mention- it will usually cause unequal spacing.  Instead, make N, the total number of points, a constraint, and then you will have the equal spacing that you need.


IIRC, etymology of "Spectro" is "of radiant energy". "Gram" is derived from the Greek word gramma, γράμμα which is a written learning.


Actually in OFDM, it is very easy. You can check Van De Beek algorithm for compensation of carrier frequency offset.


Why is the continuous time Fourier series of DC signal an impulse?

CTFS , $a_k$, of a periodic function $x(t)$ of period $T$ is :
$$ a_k = \frac{1}{T} \int_{<T>} x(t) e^{-j \frac{2\pi}{T} k t } dt \tag{1} $$
and the DC coefficient corresponds to setting $k=0$ above;
$$ a_0 = \frac{1}{T} \int_{<T>} x(t) dt \tag{2} $$
where the integration is along one period of $x(t)$. The DC coefficient cannot be an impulse, (indeed the CTFS coefficients $a_k$ form a discrete sequence, in $k$, of finite amplitude by definition of its existance.)
A signal must have finite energy per period (aka finite avg power) in order to have a valid CTFS (due to Dirichlet conditions) representation, hence the above integral in Eq.2 should be finite; i.e., the integral converges, and cannot result in an impulse.


So the point is, oversampling with a factor of $N$ requires your signal to be confined to at most $\frac1N$ of your Nyquist bandwidth. Otherwise, you're not oversampling.

If the signal is six time oversampled, then there's no loss of information/aliasing when reducing the sample rate by six.

  If I would oversample by 2 and then downsample it by 2 without AA-filter I would have stuff folding back to by final Nyquist band.

No, you wouldn't; there's, by definition of oversampling, no energy in the "upper half" of your spectrum.

Now, through discussion in the comments it became clear that "oversampled signal" might be ambiguous:

If you had your signal alone, AA filtering would do nothing, because there's no spectral component that could be folded back (aliased) into the decimated Nyquist band. However, if we're observing noisy signals, you'd have noise energy (and be it just the quantization noise) that is relatively wideband. Thus, in that situation, although the signal of interest is oversampled, the noise is not, and to avoid getting noise aliased into our band of interest, proper AA filtering is necessary.

This is really the difference between mathematically perfect signals and noisy signals; for example, if you generated a digital sine with 100 Samples per period, you could, without incurring any problems, decimate strongly – because a digitally generated signal typically has zero noise.


Possibility of temporal super resolution

Broadly similar techniques are used in a variety of applications when it's necessary to capture more samples of an event than is possible with available equipment.  The idea is really quite simple: an interpolation is just "stretching" but it does not provide any new data, which is why interpolation requires a low pass filter.  By contrast, if you add another recorder (AD converter) at a fixed delay (usually by setting up another record that is trigger by a phase shifted version of the master clock), then you are "subsampling" the analog data.

You can improve temporal resolution by using multiple ADs at phase shifts, capturing the same analog signal. You can improve spatial resolution using multiple transducers.  You can do both at the same time, although there are practical limits to how much certainty either measurement can have.

These techniques are used in probably all high end imaging: PET, airport scanners, CAT, medical ultrasound imaging, ground penetrating radar, etc.

Some of the discussion in those papers covers practical concerns that affect only video.


What could be the discrepancies between theory and practice when implementing a proportional controller?

There are no theoretical errors as you are simply describing a mathematical equation exactly. Also there are no actual 2nd order systems in the analog world due to real-world bandwidth limitations ("higher frequency" poles that we can typically ignore due to the dominant effects of lower frequency poles resulting in a reasonable estimate), and especially as your loop bandwidth increases, the effects of physical delays ($e^{-sT}$) in your system (which if you do the Taylor series expansion on that $e^{-sT}$ equation you see that it has an infinite number of poles!).

I found the following graphic useful in helping to understand the impact of the pole locations on the overall response (in particular to the point made above about "lower frequency" poles being dominant; the poles aren't actually on the frequency axis so the more accurate explanation is poles closer to the frequency axis are dominant where the frequency axis is the unit circle in the z-plane, or the imaginary axis in the s-plane as drawn below). Just look at a simple polynomial equation with two poles and you will see why it is "dominant" and how close a pole that is further out can be where it starts to matter based on the kinds of differences you are actually seeing.

When we do have a system with two complex poles which are sufficiently closer to the frequency axis than any other ones, we can make the simplifying and well documented approximations of a second order system. These are approximations.

For the graphic above which is an example Root Locus (so the poles move on the red trajectories shown as we change loop gain) I have used the following approximations:

Case 1: When the pole on the real axis is greater than 10x (in the negative direction) the real value of the complex poles, the system is "second order dominant" and it can be reasonably approximated by the two poles alone.

Case 2: When the pole on the real axis is less than 10x (in the negative direction) the real value of the complex poles, the system is "first order dominant" and it can be reasonably approximated by the one real pole alone.

Case 3: For values in between I expect to have greater errors in either approximation and likely need to model it as a third order system (depending on the accuracy of the results needed).

Here's an example from James Kamman at Western University that also used a 10x criteria but shows the significance of that in more detail:

https://web.archive.org/web/20150513061828/http://homepages.wmich.edu/~kamman/ME360DominantPoles.pdf


Zero-padding can be used to more accurately estimate the frequency of a single sinusoid in zero noise.  Since a zero-padded FFT produces a smoother looking plot, and thus shallower peaks, if there are any subtle changes in the inflections and thus peak location caused by noise, those shifts will be magnified by the zero-padded FFT, perhaps from a fraction of a FFT result bin to several bins of random shift.

So, one needs to have a good estimate of the S/N ratio to know whether and how much zero-padded will produce a usefully more accurate frequency peak estimate of a sinusoidal (or narrow band) signal.


Transfer function between voltage and decibel output for in-ear headphones

That depends to a very large degree on the type & sensitivity on the headset, the gain of the amplifier and the frequency. I don't think there is a reasonable rule of thumb, this will be all over the place.


Gardner Timing Recovery for Repeated Symbols

Your Gardner Loop would properly "Flywheel" in the absence of symbol transitions. So it would slowly drift off of synchronization, but would not spiral out of control as you describe. This is ideal for a control system with "missing updates"; to have the contributions for those updates with no additional information to contribute 0 to the accumulated error, resulting in a "flywheel" behavior. There are other detectors used to discriminate the control term in loops (such as a digital phase-frequency detector used in PLL implementations (such as this device: Maxim Phase Frequency Detector) that produce maximum error in the absence of transitions, and these systems would move off of lock at the maximum speed possible - this would be bad for timing recovery where we expect to have missing bits!

See my picture below demonstrating the operation of the Gardner Loop. In the absence of transitions, the Timing Error Detector would be zero on average as the Late and Early samples would be the same, and therefore $Late - Early = 0$.

Here is another post I made earlier with more information on the Gardner Loop:

Location of Matched Filter


An LTI system is said to be minimum-phase if the system and its inverse system are causal and stable. That's implying that all poles and zeros must be strictly inside the unit circle.
An all-pass filter's zero-pole pairs are reflections of each other across the unit circle so that an all-pass filter absorbs the poles or zeros. For a stable all-pass filter, its zeros are outside the unit circle and its poles are inside the unit circles.
A causal and stable mixed-phase system has all poles inside the unit circle and zeros which could be inside or outside the unit circle. Therefore, reflecting zeros outside the unit circle by cascading an all-pass system makes all zeros inside the unit circle, which means you get a minimum-phase system.
Once you get a minimum-phase system, the only thing you can do to prevent changing its magnitude response is to cascade an all-pass system. That brings additional zeros outside the unit circle and the system is no longer minimum-phase. Therefore minimum-phase system with a specific magnitude response is unique.


My guess is there's a problem with your wrapPhase function.

If I use the R implementation of it below (C8 shown), then I get all sorts of funny results.  That's the second image.  Removing it from the iteration (but including it in the calculation), makes for much cleaner signals (the first image).

R Code Below

#30582

wrapPhase <- function(phi)
{
  return(phi %% 2*pi)
}

# https://en.wikipedia.org/wiki/Piano_key_frequencies
cs <- c(32.7032, 65.4064, 130.813, 261.626, 523.251, 1046.50, 2093.00, 4186.01)
c1 <- cs[1]
c8 <- cs[8]

fs <- 44100

phs <- rep(0,T + 1)
phs_incr <- c8 / fs * 2 * pi

T <- 1000

sample <- rep(0,T)
sqr_wave <- rep(0,T)

for (t in seq(1,T))
{
  sample[t] <- sin(wrapPhase(phs[t]))

  if (sample[t] > 0)
  {
    sqr_wave[t] <- 1
  }
  else
  {
      sqr_wave[t] = -1;
  }

  print(sample[t])

  phs[t+1] <- phs[t] +  phs_incr
  #phs[t+1] <- wrapPhase(phs[t+1]);
}

phi <- seq(-pi, pi, 0.001*pi)
phi_wrap <- rep(0, length(phi))

for (k in 1:length(phi))
{
  phi_wrap[k] <- wrapPhase(phi[k])
}

plot(sample[1:50], col="blue", type="l")
lines(sqr_wave[1:50], col="green")


I think that the plot you display does not show the correct answer, what makes me think that is that you inject in your fft a pure sinewave with an initial phase equal to 0. Therefore your plot should display a 0 phase shift for every frequency bins.
For exemple, below you can see the phase plot of a pure sine wave with an initial phase shift of 30 degrees.

The phase spectrum will be here to tell you in which way your signals are delayed for each frequency components.


What is the difference between time complexity and sample complexity?

Some of those notions can seem quite vague, and may even depend on the field.

In algorithmics, different types of complexities come at play. One could be related to actual implementations (machine or computational complexity). Another to data access (sample, training set complexity).

Computational complexity is generally machine-, transmission-, power dependent. You have for instance time-complexity (the time need to run) or space-complexity (the memory required).

Data complexity are generally less machine dependent, and more conceptual. Sample complexity can be an estimate of the number of operations required, depending on the number of samples. In the Landau notation, an $O(N^3)$ algorithm grows as the cube of the number of samples $N$. It is a gross measure of scalabilty. Here, one is not concerned with data precision, storage, etc. If one can bound the execution time of operations, one could relate sample complexity to time complexity, but complicated computer architectures make this conversion doubtful: caching, pipelining, parallelizing,  swapping alter a lot the expectations.
In learning tasks, "data complexity" may refer to the quantity of data to be learned, on the training sets, to achieve certain prediction goals. One classical measure is the  Kolmogorov complexity.

A related SE question: Computational complexity vs other complexities.


Conceptually this works as follows: first you construct a Gaussian scale space by convolving the given image with the Gaussian filters of increasing $\sigma$. This is where the sigma comes from.  Thus you have a 3D scale volume, where each plane is an increasingly blurred version of the original image.  Then you evaluate your function $F$, which in your case would be a Laplacian, over each scale plane, i. e. at each $\sigma$.  Then you find local maxima of $F$ over $x, y$, and $\sigma$.

Note that in practice you can use the fact that convolution is a linear operation in various ways.  For example, you can generate each scale plane from the previous plane, rather than from the original image.  This allows you to use the same filter, rather using several filters with increasing $\sigma$.  Or you can generate $F(\sigma)$ directly from the original image by using the Laplacian-of-Gaussian filter of the corresponding $\sigma$


What is the average of $\DeclareMathOperator{\rect}{rect} \rect(\cos(\pi t/2)) $?

The given signal is clearly periodic. So just follow these steps and find the solution:

Find the period $T$.
Figure out the interval within a period for which the signal equals $1$.
Compute the average: $$\overline{x(t)}=\frac{1}{T}\int_0^Tx(t)dt$$

If you're convinced that you did everything right, don't worry about a given solution which is different.


Why do FIR systems need to be antisymmetric about M/2 to be linear phase?

Consider the $M$ th order FIR system, where $M$ is even
$$\begin{aligned}
H(z) &= \sum_{m=0}^{M-1}b_m z^{-1} \\
     &= b_0 + b_1z^{-1} + \ldots b_{M-1}z^{M-1} \\
H(e^{j\omega}) &= b_0 + b_1 e^{-j\omega} + b_2 e^{-2j\omega} + \ldots b_{M-1}e^{-j(M-1)\omega} \end{aligned}$$
Now if our coefficients are antisymmetric, then $b_m = -b_{M-1-m}$, or, $b_0 = -b_{M-1}, b_1 = -b_{M-2}, \ldots, b_{\frac{M}{2}} = -b_{\frac{M}{2}-1}$. This gives us
$$\begin{aligned}
H(e^{j\omega}) &= b_0 - b_{M-1}e^{-j(M-1)\omega} + b_1 - b_{M-1}e^{-(j(M-1) + 1)\omega)} + \ldots \\
&= e^{-j(M-1)\omega} \left[ 2b_0 + b_1(e^{-j\omega}-e^{j\omega}) + \ldots b_{\frac{M}{2}} (e^{-j\frac{M}{2}\omega}-e^{j\frac{M}{2}\omega}) \right] \\
&= -2je^{-j(M-1)\omega}) \left[jb_0 +b_1 \sin \omega + \ldots + b_{\frac{M}{2}} \sin{\frac{M}{2}\omega} \right] \\
&=
-2je^{-j(M-1)\omega})\left[jb_0 + \sum_{m=1}^{M/2} b_m \sin(m\omega)\right]
\end{aligned}$$
We see that $H(e^{j\omega})$ is linear phase, with its phase $ = -(M-1)\omega - \frac{\pi}{2}$. This would have not worked out if the filter coefficients were not antisymmetric or symmetric.


How can we combine 2D image with 3D image for the same sense?

Formally speaking, you would like to extrinsically calibrate the laser scanner to the 2D image. I have taken the liberty to edit your question to reflect that.

Here is how my initial approach would be:

Calibrate the intrinsics of the 2D camera. For that, just use OpenCV. You should store the intrinsic parameters: focal lengths, principal point and distortion coefficients.
Take an object which you could easily locate in 3D and 2D. Such an object could be any 3D model for instance if you can perform pose estimation both in 3D and 2D. Because it would be difficult to do this for free-form objects, I would propose to use a sphere-object which appears to be ellipse-like in the camera, and a sphere in the laser scanner. For instance, take an unused pilates-ball.
Place the pilates ball in the overlapping area of the point cloud and the 2D image. Then capture multitudes of images (at least about 30) by moving the ball in the overlapping region to various locations. At each location, make sure that the ball occupies about $1/4^{th}$ of the image and record the 2D and 3D frames.
Detect the ellipses in all the images, i.e. for each image, one ellipse. Use any open source tool, for instance 1D-hough transform to do that. I strongly suggest to refine the ellipse to subpixel accuracy, for instance via this method (source code). Record the center of the sphere (ellipse) in 2D $\{(x_i,y_i)\}$.
Detect the spheres in all of the laser point clouds. Again, many packages do exist. Even MATLAB has one method. Record the detected centers of the 3D sphere $\{(X_i,Y_i,Z_i)\}$. Again, a non-linear refinement of the sphere in 3D can help you in getting more accuracy. Check this one out for a complete pipeline on this if you like more sophisticated approaches.
Now you have an easy problem to solve: You have established 2D-3D coorespondences in the form $\{(x_i,y_i)\}\sim\{(X_i,Y_i,Z_i)\}$, have the calibration information of 2D camera and thus, all you need to do is to solve the PnP problem. Again many open and closed source software exists such as OpenCV's solvePnP or MATLAB's estimateWorldCameraPose.

The result will be a rotation matrix $R$ and a translation vector $t$ that maps the 3D scanner's coordinate frame to the image domain, thereby solving the registration problem.

I hope instructions are clear and this helps.

Note that, many things can be done to make this procedure easier, more robust, more accurate and etc. This created quite a big interest from academia, and if you dig a bit, you will come across a lot of sources. Nevertheless, for the beginning, the aforementioned procedure is pretty easy to implement.


Vivado Design Suite Reed-Solomon Decoder

$\left\lfloor\frac{n-k}2\right\rfloor$ is the maximum number of bit errors that we can guarantee a code can correct for a binary channel.
But that doesn't mean a code can't sometimes correct more! (In general: not all code words have a distance that's equal to the code's minimum distance. It is the case for minimum distance separable codes, like RS codes are, though.)
Then:
RS codes are also used in applications where erasures, not bit flips, happen. You'll find the decoder documentation contains a section on erasure decoding.
You'll also find that a RS code can correct up to $n-k$ erasures – and my guess is that they made the count of that and the count of corrected errors the same width for design consistency reasons, and maybe to make different decoder IP interchangeable.
Note that giving a unsigned number more bits than it needs has no functional or FPGA utilization downside in IP cores: if whatever takes the output uses less bits than it has, the synthesizer will optimize away all calculations resulting in the bits that weren't used.


Moving from deterministic signals to stochastic signals in s-domain (Power Spectral Density)

You have to look at the autocorrelation function of $y(t)$:
$$R_y(\tau)=E\{y(t)y(t+\tau)\}\tag{1}$$
with
$$y(t)=(h_A\star a)(t) - (h_B\star b)(t)\tag{2}$$
where $\star$ denotes convolution. If you write out $(2)$ with integrals and plug it into $(1)$ then, with the given assumptions on $a(t)$ and $b(t)$, you'll see that the mixed terms with the negative sign are zero, and the two remaining terms with a positive sign become
$$R_y(\tau)=\sigma_a^2r_{h_A}(\tau)+\sigma_b^2r_{h_B}(\tau)\tag{3}$$
where $\sigma_a^2$ and $\sigma_b^2$ are the variances of the processes $a(t)$ and $b(t)$, respectively, and $r_{h_A}(\tau)$ and $r_{h_B}(\tau)$ are the deterministic autocorrelation functions of $h_A(t)$ and $h_B(t)$, respectively. Finally, taking the Fourier transform of $(3)$ gives
$$S_y(\omega)=\sigma_a^2|H_A(j\omega)|^2+\sigma_b^2|H_B(j\omega)|^2\tag{4}$$
for the power spectrum of $y(t)$.


Using matlabs tf2cl, it may require the DSP toolbox, you can turn the transsfer function (numerator and denominator coefficients) in to a cascade allpass lattice coefficients, which is a two phase filter.

% Generate some IIR filter coefficients
[b,a]=cheby1(9,.5,.4);

% Coupled allpass decomposition
[k1,k2]=tf2cl(b,a);


How could I approach determining if this 2D system represented as a 2D summation formula is linear?

I'm actually starting this out not knowing the answer, although I have my suspicions.  I think it'll be a linear system, but I'm not sure.  Here's your starting point: you define the system $y = h(x)$ as:
$$y = h(x) \mid y(m,n) = \sum_{k_1=-\infty}^{m} \sum_{k_2=-\infty}^{n} x(k_1,k_2)$$
So let
$$\begin{aligned}
y_1(m,n) = \sum_{k_1=-\infty}^{m} \sum_{k_2=-\infty}^{n} A\,x_1(k_1,k_2) = h\left(A\, x_1 \right) = A h(x_1) \\
y_2(m,n) = \sum_{k_1=-\infty}^{m} \sum_{k_2=-\infty}^{n} B\,x_2(k_1,k_2) = h\left(A\, x_2 \right) = B h(x_2)
\end{aligned}$$
Then
$$
y_1(m,n) + y_2(m, n) = \sum_{k_1=-\infty}^{m} \sum_{k_2=-\infty}^{n} A\,x_1(k_1,k_2) + \sum_{k_1=-\infty}^{m} \sum_{k_2=-\infty}^{n} B\,x_2(k_1,k_2)
$$
Because the summations have identical limits, we can bring everything under one pair of summations:
$$
y_1(m,n) + y_2(m, n) = \sum_{k_1=-\infty}^{m} \sum_{k_2=-\infty}^{n} \left ( A\,x_1(k_1,k_2) + B\,x_2(k_1,k_2) \right)
$$
but this is just $y_1 + y_2 = h(A\,x_1 + B\,x_2)$.  You already know that $y_1 = A h(x_1)$ and $y_2 = B h(x_2)$, so linearity is proved.
(Which is what I thought -- I was going to just dismiss your original problem as "summations are linear, go away!" -- then I realized I wasn't sure).


The last expression (sum of  a geometric series) is a common abuse of notations: it should have been:

$N$ if $k= k_0$
$f(r)=\frac{1-r^N}{1-r}$ with $r=e^{-j 2\pi(k-k_0)/N}$ if $k\neq k_0$

Indeed, as you correctly remarked, numerator and denominator would vanish for $r=1$ (or $k= k_0$), so the fraction is not "theoretically" defined. However, it is consistent to the limit, as $f(r)\to N$ as $r\to 1$  since $e^{-j 2\pi(k-k_0)/N}\to 1$ when $k \to k_0$ (allowing real $k$).

Mathematically, the fraction gives one single expression, defined by continuity at $k= k_0$. It is a bit like saying that $\frac{r^2-1}{r-1}$, not defined at $r=1$, is in some way equivalent  to $r+1$ (using $r^2-1=(r-1)(r+1)$  and wrongfully simplifying the fraction. This is not correct, but makes sense.


Is there any way to differentiate between impulses?

If you are looking at vectors $\mathbf x$ and $\mathbf y$
(presumably of the same length) in which the entries are
mostly $0$s and

with occasional $1$s (unit impulses), a good measure of the difference
between the vectors is the Hamming distance between them, that is, the
number of samples in which $\mathbf x$ and $\mathbf y$ differ
with occasional non-zero entries (impulses of different amplitudes),
then we need to distinguish between impulses of different amplitudes
but the same location. If we use Hamming distance, then the measure
of the difference between
$$\cdots \quad 0\quad 0 \quad 0 \quad 1 \quad 0 \quad 0\quad 0\quad \cdots$$
$$\cdots \quad 0\quad 0 \quad 0 \quad 2 \quad 0 \quad 0\quad 0\quad \cdots$$
and
$$\cdots \quad 0\quad 0 \quad 0 \quad 1 \quad 0 \quad 0\quad 0\quad \cdots$$
$$\cdots \quad 0\quad 0 \quad 0 \quad 9 \quad 0 \quad 0\quad 0\quad \cdots$$
is the same which is not quite capturing our intuition that the two pairs
are more different than the Hamming distance is telling us. Possibilities
to consider in this case are
$$\sum_{i=1}^n |x_i - y_i|,\quad   \sum_{i=1}^n |x_i - y_i|^2,\quad
\mathrm{or}\quad \sqrt{\sum_{i=1}^n |x_i - y_i|^2}.$$
The last quantity is the Euclidean distance between the vectors. Which
measure you want to use depends on how much the measure captures the
features relevant to your application.


Why Does FISTA Algorithm Not Work for Signed Signals?

I will try explaining why would someone use such an option as defining the solution as positive.
The Proximal Gradient Method which is the hurt of the FISTA algorithm is basically a generalization of Projected (Sub) Gradient Descent Method.
As such, if one knows that the solution lies in a Convex Set and one knows the projection onto that set one could utilize it inside the method by a projection step.

Hence, if you know your solution is indeed in the convex set of the positive orthant you should use that information by projection the solution onto this set.
The projection is fairly simple:

$$ {P}_{ \mathbb{R}_{+}^{n} } \left( y \right) = \arg \min_{x \in \mathbb{R}_{+}^{n}} {\left\| x - y \right\|}_{2}^{2} = \max \left( y, 0 \right) $$

This is actually implemented in the package you linked to in the Proj files.
For instant, in the latest version of proj_l1.m you can see the code:

if opts.pos
    X = max(0, U - lambda);
else
    X = max(0, U - lambda) + min(0, U + lambda);
end

Which implements the projection onto the $ {L}_{1} $ norm ball and the positive orthant as written above.


FIR Filter - Low Pass Filter coefficients to High Pass filter coefficients

HINT:

It looks like they shifted the frequency response by half the sampling frequency, i.e.,

$$H_{HP}\left(e^{j\omega}\right)=H_{LP}\big(e^{j(\omega-\pi)}\big)\tag{1}$$

Frequency shifting corresponds to modulation (multiplication) in the time domain. Now you just have to figure out the modulation sequence that achieves the correct frequency shift.


Linear Block Codes Digital Communication Problem

In this answer, I will focus on error correction. I will use as little coding theory as possible, and I'll try to approach the problem from the point of view of a novice taking their first steps in code design.

Our objective is to find the best (8,4) linear block code possible, that is, one that corrects as many bit errors as possible. Is it possible to know in advance how good an (8,4) code can be? To answer this question, we have an invaluable tool: the Hamming bound. It says that a $(n,k)$ code that can correct $t$ errors must satisfy

$$ \sum_{i=0}^t \binom{n}{i} \leq 2^{n-k}. $$

It is easy to write a program that takes $n$ and $k$ as arguments, and outputs the largest $t$ that satisfies the bound. For $n=8$ and $k=4$, we find that the largest possible $t$ is 1. Note that, since the minimum distance is given by

$$ d_{min} = \lfloor (t-1)/2 \rfloor, $$

we still don't know whether the best code has $d_{min}=3$ or $d_{min}=4$. Both will have the same error correcting capabilites, but the code with minimum distance 4 has additional error detecting capabilities.

(Note that, even if there exist numbers $(n,k,t)$ that satisfy the Hamming bound, it is not guaranteed that a $(n,k,t)$ code actually exists. However, if the bound is not satisfied, we know for sure that no such code exists.)

Now that we have a goalpost ($d_{min}=4$), we can focus on finding a code that achieves it. We know that the generator matrix $G$ has $k$ rows and $n$ columns. In its systematic form, $G=[I_k|P]$, where $P$ has $k$ rows and $n-k$ columns. For our (8,4) code, we have to find a matrix $P$ with 4 rows and 4 columns.

Since we don't have any other tools at this point, we can try a brute-force approach, since there are "only" $2^{16}$ possible matrices $P$ (actually a bit less, since we don't want to repeat any columns). Let's try this matrix:

$$ P=\begin{bmatrix}1 & 1 & 1 & 1 \\ 0 & 1 & 1 & 1 \\ 0 & 0 & 1 & 1 \\ 0 & 0 & 0 & 1\end{bmatrix}. $$

It is not hard to write a program that generates all 16 possible code words using $P$. The final step is to verify if our code has $d_{min}=4$. This is easy, since the code is linear; all we have to do is find the non-zero code word with smallest weigth (with smallest number of ones). For the matrix $P$ given above, it turns out that $d_{min}=2$, so the code is not the best possible.

It is easy to find generator matrices with $d_{min}=3$; extending the distance to 4 takes a bit more work.

In conclusion:

The Hamming bound can be used to determine what values of $t$ and (approximate) $d_{min}$ we can expect from a code with small $n$ and $k$. For large $n$ and $k$, the achievable value of $t$ may be much less than that given by the bound. Surprisingly, determinining $t$ for given $n$ and $k$ is (for the general case) an unsolved problem.
Design by brute force (by trying many $P$ matrices) is feasible only for small codes. Even for an (8,4) code you may have to try more than 60,000 possibilities. The number grows exponentially with the code size. The amount of work required to generate the entire code, and to determine $d_{min}$ by inspection, soon grows so large that it becomes entirely unfeasible.
It is for this reason that so much research has been done into techniques that allow for less complex design of large codes. When you study other codes, such as cyclic and convolutional codes, keep in mind this motivation. You probably shouldn't try to solve your problem using the brute-force technique outlined here except to get a feel for the sheer magnitude of the computational complexity involved in the design of good codes, and also in their encoding and decoding procedures.


how can I analyze a basic sine wave for its frequency phase and amplitude at a particular point in time in matlab

If your the unknown signal $x(n)$ is modeled as:
$$x(n)=A \sin(2 \pi f n+ \phi)+w(n)$$
and you want to estimate $A$,$f$,and $\phi$ accurately, you can use least square estimation. Unfortunately the cost function is nonlinear. You can use nonlinear least square in MATLAB to find the parameters as follows:

Make a cost function:

f=@(A,f,phi) x[n]-A sin(2*pi*  f *n+ phi)

and use

p0=[A0,f0,phi0];

p = lsqnonlin(f,p0);

to find the unknown parameters.
Note that the optimizer will have a hard time finding $f$ as the problem is not convex. So it is best if you can give an initial estimate of frequency by using a method like fft.

If the frequency is known then the problem can be converted to a linear estimation as:
$$x(n)=A \sin(2 \pi f n+ \phi)+w(n)$$
$$=A \sin(2 \pi f n) \cos (\phi) +A \cos(2 \pi f n) \sin (\phi)+w(n)$$
$$=p_1 S[n] +p_2 C[n]+w(n),$$
where $p_1=A \cos(\phi)$ and $p_2=A \sin(\phi)$ are unknown parameters, and $S[n]$ and $C[n]$ are known.


Finding the coefficients of the block diagram

If you assume that the samples lie on a line then the current sample is related to the previous sample by
$$x[n]=x[n-1]+c\tag{1}$$
with some unknown constant $c$. If there are two past samples given, then $c$ is simply determined by the difference between these two samples:
$$c= x[n-1]-x[n-2]\tag{2}$$
Combining $(1)$ and $(2)$ gives
$$x[n]=2x[n-1]-x[n-2]\tag{3}$$


Since FFT treats the signal as if it is periodic you need either to apply a window function (for example hanning) on your signal or make it coherent. In the image3 you attached you may make it coherent by only using data for a number of cycles; use samples 0..31 or 0..63. If the signal is non-coherent it will be seen by the FFT as a concatenated signal with discontinuites that will cause overtones.

EDIT:
A coherent signal may be concatenated with itself without generating any discontinuites.

The blue line is the input data to FFT

And a non-coherent:


What is the right way to calculate the inverse Z-transform of $zX(z^{-1})$

The first result is correct, the second is wrong. You're right that $x_i[n]=x[n-1]$ corresponds to $z^{-1}X(z)$. And if you replace $z$ by $1/z$ you need to replace $n$ by $-n$. But $x_i[-n]$ is simply $x[-n-1]$,  and not $x[-n+1]$. Note that you don't replace the whole argument of $x[n-1]$ by its negative value, but you just replace $n$ by $-n$.


Region of the coefficients of a quadratic equation that cause the roots of it to be in the unit disk

It may be helpful to look at how the roots and the coefficients are connected to each other. Let's start with a complex root. We can write this as
$$ H(z) = (1-r/z) \cdot (1-r'/z) = 1 - (r+r')/z +r\cdot r'/z^2$$ $$ = 1 -2 \cdot Re\{r\}\cdot z^{-1}+ |r|^2 \cdot z^{-2}$$
We can see that $a1$ is simply twice the negative real part of the root and that $a2$ is the magnitude squared of the root. In this case the condition $|a_2| < 1$ is sufficient and $a_1 is automatically taken care off by the constraint that the root is complex.
Let's do the same for a pair of real roots, $r_1$ and $r_2$
$$H(z) = (1-r_1/z) \cdot (1-r_2/z) = 1-(r_1+r_2) \cdot z{^-1} + r_1 \cdot r_2 \cdot z^{-2} $$
So here $a1$ is the negative sum of both roots and $a_2$ is the product of the roots. Again, we get $|a_2| < 1$ since the product of two numbers smaller than one cannot be larger than one.  If both roots have the same sign, we can write the constraint as
$$ 0< (r_1-1) \cdot (r_2 -1) = r_1 \cdot r_2 - r_1 - r_2 +1 = a_2+a_1 + 1$$
which gives you
$$-1 < a_1 + a_2$$
Doing the same exercise for mixed signs gives you the other other inequality.


why this frequency response is added to the power of e?

$H(\Omega)$ is a complex number. $Ae^{j(\Omega t + \phi)}$ is another complex number, expressed in polar format by its magnitude $A$ and phase $\Omega t + \phi$. When you multiply complex numbers in this format, you simply multiply their magnitudes and sum their phases. Thus:

$$
H(\Omega)Ae^{j(\Omega t + \phi)} = |H(\Omega)| Ae^{j(\Omega t + \phi \angle H(\Omega))}
$$


How to predict next image in a sequence, or predict the deblurring of an image?

There is no extra information in the 720p and 480p images that is not already in the 1080p image. You can do the interpolation on the 1080p image. A recent work that may be of interest to you is Ledig et al. Photo-Realistic Single Image Super-Resolution Using a Generative Adversarial Network, arXiv:1609.04802v5, 15 Sep 2016 (v1), last revised 25 May 2017. But if your work is for scientific purposes rather than for looks, then it is probably best to use traditional interpolation methods that do not try to increase the resolution by making up new data.


how to correct minBLEP sawtooth oscillator DC bias

If your oscillator is steadily accumulating DC over time, adding a DC blocker at the output may not be an elegant solution since your oscillator will keep getting more and more DC.  For floating point this may or may not be a huge problem (I wouldn't do it), but for fixed point you will eventually start clipping or wrapping.

If the oscillator has a frequency dependent DC offset that doesn't creep up over time, a DC blocker is probably a good solution.  If you've implemented the BLEPs correctly, this should be the case.

The original reason given for BLEPs to be minimum phase was to implement oscillator sync.  To implement this correctly with linear phase BLEPs, you should use lookahead so you can start the BLEP before the aliased reset occurs.

However, apparently this author would rather not have any lookahead for some reason and therefore designed minBLEPs so that you would have to lookahead less far into the future.  If I remember correctly, in the paper they say that it's okay to just forget about lookahead if you use minimum phase, but obviously this isn't perfectly true as minimum phase filters still have some delay.  You should experiment yourself to see if the approximation works for you.

In my opinion, it's probably a good idea to use linear phase BLEPs and lookahead appropriately for sync.  We're talking a few samples of latency here for reasonable BLEPs.


Recursive, non-recursive systems; FIR, IIR systems

The logical implications are the following:
"non-recursive" $\Longrightarrow$ FIR
IIR $\Longrightarrow$ "recursive"
But the opposites are not necessarily true because a FIR system can be implemented recursively (transfer function poles can be cancelled by zeros).

Of course, when referring to "recursive" or "non-recursive" we always talk about implementations with finitely many operations per output sample. Clearly, any discrete-time LTI system can be described by a generally infinite convolution sum, but that is not what we mean by "non-recursive".


EEG signal processing with wavelet or fft?

Here are links to the coding.  I am doing something else, so I have not dug into the code to see how the algorithms work.  If you need help analyzing the code set up a private "room" here and we can discuss it there.  I am sure that it requires more chit-chat than the hall monitors like here :)
a simple Processing application
ARDUINO BRAIN LIBRARY
Your picture
Thanks for the lead!  I had an application for treating certain dementia problems using something like this but haven't followed up.