This is a common problem, and has an easy solution. The following assumes you are using a cascade of second-order direct-form-1 sections.
 For a Butterworth highpass filter with a small ratio of Fcutoff to Fs, the numerator coefficients are close to [1 -2 1], and therefore the frequency response of the numerator by itself has a large attenuation at low frequencies. This is compensated by the gain of the recursive section (the poles).
When a low-frequency sinusoidal signal is first applied, the input shift register (corresponding to the numerator) is not filled up until the 3rd clock cycle, and therefore the output of the numerator FIR may be large for the first 2 clock cycles, after which it becomes small. These initial large values get captured in the recursive section, and recirculate for some time until they die out.
So the simple solution is to delay the operation of the filter until the input shift register is filled with the first 3 samples.


Are there any algorithms that, for example, take neighboring pixels into consideration in order to determine a better estimation for each pixel?

That would essentially be a low-pass filter.
So, yes, that exists, and is commonly used.
You can see that very nicely if you take e.g. an old computer graphics sprite, and scale it as a high-color-resolution image:
original (32 x 32 px)

scaled without interpolation filter (512 px)

scaled by a simple linear interpolation filter (512 px)

While the simply scaled image has as much discrete colors as the original, the linearly interpolated one has a lot more colors, and the histogram looks much more "continuous".
But:
There's no free lunch, usually. You can't just add color information back! (not being able to add back lost info is one of the fundamental truths of Information Theory, by the way). You have to sacrifice spatial resolution. Essentially, it's the same math that underlies Heisenberg's Uncertainty Principle. You can either have perfect info state of something (the color of a pixel, the impulse of a electron), or its position (sharp edges in an image, knowledge of the position of an electron).  See this example:
original photo

original (512 x 512px) reduced to 16 colors

filtered to achieve a much more continuous color histogram, with a radius of 10px taking into consideration


Phase response of moving average filter -- how to interpret?

The frequency response of a causal length $N$ moving average filter is

$$H(\omega)=\frac{\sin\left(\frac{N\omega}{2}\right)}{N\sin\left(\frac{\omega}{2}\right)}e^{-j\omega(N-1)/2}=A(\omega)e^{j\phi(\omega)}\tag{1}$$

Note that $A(\omega)$ is not the magnitude of $H(\omega)$, but it is a real-valued amplitude function, which takes on positive as well as negative values. The phase $\phi(\omega)=-(N-1)\omega/2$, as defined in $(1)$, is obviously linear. That's also the common definition when we talk about a linear phase response.

The phase you plotted is not $\phi(\omega)$, but $\hat{\phi}(\omega)$ as defined by

$$H(\omega)=|A(\omega)|e^{j\hat{\phi}(\omega)}\tag{2}$$

The difference between $\phi(\omega)$ and $\hat{\phi}(\omega)$ is that whenever $A(\omega)$ crosses zero, a phase jump of $\pm \pi$ occurs in $\hat{\phi}(\omega)$, corresponding to a sign change in $A(\omega)$. Nevertheless, we still refer to $H(\omega)$ as a frequency response with a linear phase, because $\phi(\omega)$ is a linear function of $\omega$.

Note that in practice, a linear phase is only relevant in the pass band of a filter, i.e. in a frequency region where no zeros of $H(\omega)$ occur. In the pass band, also $\hat{\phi}(\omega)$ is linear, because it only jumps at the zeros of $H(\omega)$.


Wiener filtering is one approach. It might even be the best approach. A Wiener filter is designed to minimize the noise (in the least squares sense) and invert the effect of the impulse response, given a known signal and a signal that is known to be tainted with noise and an impulse response. Once you have a Wiener filter, you can then compare the amplitude of the filtered signal with the amplitude of the direct signal to estimate the noise (I think).

Now I've never implemented a Wiener filter myself, so I will defer to a textbook with chapter dedicated to the subject: Advanced Digital Signal Processing (Electrical and Computer Engineering) by Glenn Zelniker and Fred J. Taylor. It's out of print, but I'm sure you can pick up a used copy for cheap somewhere. I suggest it because it's more mathematical than engineering in it's approach and that might appeal to you. Many other textbooks have info on the subject.

If a more ad-hoc approach suits you, here is a suggestion. I've never tried it, but it might work:

Create a filterbank, F, (a properly windowed FFT will work) to analyze both the known music file A and the input C. You will also need some measurement function, M. I would suggest $M(x) = | x |$, rather than $M(x) = x^2$, but both will work, and you may also want to do some smoothing over time. By comparing the results of $M(F(C))$ to $(M(F(A))$ when the noise is low, you should be able to determine the effect of B. Of course, you may not know that the noise is low, in which case, you'll have to get a bit more clever, creating some estimate based on some statistics. The goal, however, is to get to a point where you can predict the output of $(M(F(C))$ from $(M(F(A)))$, to within some bounds. The lower the bounds the better you will be at detecting the noise.

To estimate the noise, you will compare how far above those bounds the signal has actually gone. This will work well for band-limited noise, but for quieter, broadband noise, it may not work well. In that case, you may want to create an aggregate statistic based not just on how far above the bounds a signal has gone, but how many bounds have been crossed, or simply apply your measurement function directly to your unfiltered signal: $M(C)$ vs $(M(A))$.


I don't think it's a good idea to add random noise to frequency response. In general, noise is added to the input or output signal. You should be clear that at which stage the noise is introduced. According to the different stages of noise introduction, there are two main transfer function estimators, $H_1$ estimator and $H_2$ estimator.

$H_1$ assumes that there is no noise on the input and consequently that all the input measurements are accurate. All noise is assumed to be on the output.

$$
\hat{Y}(\omega) = H(\omega) X(\omega) + N(\omega)
$$
where $\hat{Y}(\omega)$ is the observed output signal, which is polluted by additive noise $N(\omega)$. And you will see that there is a deviation between the estimated frequency response and the real one.
$$
\hat{H}(\omega) = \frac{\hat{Y}(\omega)}{X(\omega)} = H(\omega) + \frac{N(\omega)}{X(\omega)}
$$
This estimator tends to give an underestimate of the frequency response if there is noise on the input.

$H_2$ estimator assumes that there is no noise on the output. Noise is assumed to be only on input.

$$
Y(\omega) = H(\omega) \Big[X(\omega) + N(\omega)\Big]
$$
This estimator tends to give an overestimate of the frequency response if there is noise on the output.
Check here Transfer function estimate


Why do the OFDM training (and payload) symbols have silent sub-carriers (aka virtual carriers)?

From the first to the sixth OFDM-Symbol you have spread virtual carriers equally across [0,..,51] with a distance of 3. With this way you have spread your 'energy' loss over a wide spectrum. In difference to having a single carrier serving as virtual carrier all the time you can transmit information over every 'data' carrier and maybe have a constant mean energy level across all subcarriers.

You could use this virtual carries to detect carrier frequency offset (CFO) since you know the pattern. In the worst case you have to store 6 OFDM-Symbols to find the beginning of your pattern.

If you have an OFDM-Symbol with 64 subcarriers you can set 6 carriers at each end as subcarrier and use the other 52 with this virtual carrier pattern to transmit information with equalization.

I think this virtual carrier pattern is not connected to the Schmidl & Cox synchronsiation pattern.


How replicas are formed in Frequency domain when a signal is sampled in Time Domain?

Some people treat sampling in the time domain as the multiplication of a signal $x(t)$ by a periodic impulse train $\sum_n \delta(t-nT)$. Since the Fourier transform of a periodic impulse train is a periodic impulse train $\sum_k \delta\left(f-\frac{k}{T}\right)$ (I might be off by a factor of $T$ here, but the basic idea is correct), and
multiplication in the time domain corresponds to convolution in the frequency domain,
we get that the Fourier transform of the sampled signal is
$$X(f)\star\sum_k \delta\left(f-\frac{k}{T}\right) = \sum_k X(f)\star\delta\left(f-\frac{k}{T}\right) = \sum_k X\left(f-\frac{k}{T}\right)$$
which is a periodic function of $f$ since $X(f)$ is replicated at intervals of
$\frac{1}{T}$ Hz along the frequency axis.
If $X(f)$ is bandlimited to $\left(-\frac{1}{2T},\frac{1}{2T}\right)$, the replicas
do not overlap.

Other people resolutely refuse to consider impulse trains and come up
with other, equally valid, means of arriving at the same result.


HINT:

Answer these two basic questions, and you should be able to solve your problem:

What's the sign of the term $2\cos(2\omega)-2$? Can it be a magnitude?
What does a term of $\pm \pi$ in the phase do to the sign of an expression?


Kalman Filter - Deriving state transition function

Hi Matthias La: I hate to be critical but the first example at the link you provided is actually quite poor because they end up using exponential smoothing (their update for $\hat{x}$ can re-written as an exponential smoothing update) for the update of the state, $x_t$, without deriving it. The way to derive the exponential smoothing update for the state (They also come up with a heuristic formula for the exponential smoothing factor, $\alpha_{n}$, of $\frac{1}{n}$ which is not optimal. It actually ends up being a function of the SNR) is to write the whole system in the following way:
$y_t = x_t + \epsilon_t $        # Measurement Equation
$x_t = x_{t-1} + \gamma_t$         # State Equation
where $\epsilon_t$ and $\gamma_t$ are independent normal error terms.
With the framework above, if you derive the updating equations for the kalman filter, you will obtain the exponential smoothing update for the state ( and $y_t = \hat{x_t}$ as update for the measurement). In the statistics-econometrics literature, the KF setup above is called a "random walk + noise" model.  I'm not claiming that the  link is completely useless but, as far as an intro to KF, that is not a good first example because they pull things out of the sky and don't explain where they come from. That can be confusing when going through the later examples.
I hope the setup above helps to show how  to obtain the ES update. Oh, I recommend below for a dense look at the KF. Harvey has a terse style which can be difficult at first (it was for me) but, in the end ( after a few reads:) ), its hits you that he actually provides a very nice explanation. I've read various parts of it more than a few times.
https://www.amazon.com/Forecasting-Structural-Time-Harvey/dp/0521405734/ref=sr_1_1?dchild=1&keywords=Andrew+Harvey+Kalman+Filter&qid=1615252088&sr=8-1
EDIT: I didn't mention it above but I remembered later that you'll need a method for estimating the measurement variance and the state variance if you want to run the updating equations.  Harvey uses prediction error decomposition but my impression is that dsp people have a different way of estimating them. Prediction error decomp is the only method I'm familiar with but you may know of a better one. If you do and can explain it to me, it's appreciated.


How to calculate the instantaneous frequency of a signal?

Instantaneous frequency is a poorly defined concept over wide-band (single channel) data, but can be a useful estimation given narrow-band signals of sufficient duration.  One trick for using an instantaneous frequency estimation algorithm which is usefully accurate only in one frequency band is to resample your data to move the frequency band of interest into the bandpass of the algorithm.  You may need to bandpass filter your data before resampling so that the source band ends up entirely within the algorithm's bandwidth.

In your case, you can try high-quality resampling your data to double the sample rate to extract instantaneous frequency estimates from 50 to 100 Hz.  Triple for 75 to 150, etc., assuming the sample rate and anti-aliasing mechanism captured spectrum in those ranges.


Need help implementing DCT Type 4 algorithm

I found the solution. There was a mistake in the paper.

On the first 2 lines of Algorithm 2:

w0 ← x0
vN/2 ← xN−1

x0 and xN-1 should be scaled by 2, giving:

w0 ← x0 * 2
vN/2 ← xN−1 * 2

The following paste has working code (and an open source license, so feel free to copy it)

https://pastebin.com/TvgU4b6d


Power in a frequency band: Spectral estimation versus time domain

Personally I'd go with the bandpass filter, unless there is a strong reason not to.

Much lower implementation complexity, CPU and memory usage.
You can get a continuous output signal at the same sample rate as the input signal. The FFT method will get you only one output per frame, so it's down sampled to the hop size.
No framing artifacts based on how your signal lines up with the frame grid: You can frame it, but you don't have to.
You can design fairly flexible "ballistics" into the detector: fast attack, slow decay, peak hold, lossy peak. I.e. it's easy to customize how fast/slow it reacts to different types of energy changes.

On the downside the choices for frequency weighting functions are fairly limited. If you need fancy frequency domain shapes the FFT may be better.


Current harmonics via FFT -- do peaks matter?

You can interpolate the FFT bin results to get better estimates of the actual peaks in the spectrum.  Use a high quality interpolator, such as a windowed Sinc kernel.


How to partition collagen into separate fibers?

Your data seems to be quite noisy in a sense that there are a lot of holes in the fibers. When you know a typical minimal size of a fiber (probably its widths), you could use a binary closing filter:

Dilate the image with a structural element of the typical width of a fiber. Most likely a circle with a diameter equal to the minimal fiber widths could be appropriate.
Erode the image with the same structural element.
Use a skeletization algorithm on this data set.


Difference between filters in serial and parallel

If you're talking block diagrams, a pair of filters in series would yield
$$H(z) = H_1(z) H_2(z)$$
The only sensible way to connect filters in "parallel" in a DSP application would be to sum them:
$$H(z) = H_1(z) + H_2(z)$$
The response would be completely different (work it out with a couple of different transfer functions for $H_1$ and $H_2$ and see what I mean).
Analogies should be an aid to understanding, not a barrier.  That book must be assuming folks who are familiar with circuits.  So I'm going to give you a short course -- but finding a book that doesn't assume an intimate knowledge of circuits may be a good idea for you.
Here's two lowpass filters and two highpass filters, and two filters usin resistors and caps, and two filters using resistors and inductors.
Basically, a capacitor is an open-circuit at low frequencies, and a short at high frequencies.  So when you arrange the capacitor as a shunt to ground, at high frequencies all the signal goes to ground, but at low frequencies it all goes to $V_{out}$.  Conversely, when you use a cap as a series element and a resistor to ground, at low frequency the cap blocks the signal, and the resistor passes whatever is left to ground.
An inductor is the opposite -- it's a short circuit at low frequencies, and an open circuit at high frequencies.  So the circuits are opposite.
But I want to stress here -- if what I just said still confuses you, get a different book.  DSP is about math, not electrical circuits.  If the author is relying on an intuition about electrical circuits to give you an intuition about DSP, and you don't have an intuition about electrical circuits, then you need a book that's done a better job aiming itself at you, not the electrical engineer standing next to you.
Even if it's for a class, you may want to get a different book as a suppliment.


What do cloud-based speech-to-text services upload to the cloud?

Its probably safe to assume that all of the major companies send enough information to reconstruct the audio. This is because having that audio for training is such a valuable resource. A certain percentage of the audio segments will be listened to and transcribed by a human annotator.

Also features in these systems are more complicated than MFCCs. You obviously want to offload computing to the cloud so that's another reason to expect that they send the whole audio.


FFT: Does the result of $N \log_2(N)$ stand for total operations or (complex) adds and multiplies?

{EDIT: The original question concerned a $1920\times 1080$ FFT, not a $2040\times 2048$ FFT. So the parts of this answer that deal with sizes that are not powers of two are no longer relevant. I will leave them anyway. The last part of this answer, dealing with complex vs. real operations, remains relevant.}

A 2-D FFT of size $N_1\times N_2$ is equivalent to performing $N_1$ FFTs of size $N_2$, and $N_2$ FFTs of size $N_1$. So the cost becomes
$$C_{2D}(N_1,N_2) = N_1C_{1D}(N_2)+N_2C_{1D}(N_1),$$ where $C_{1D}(N)$ is the cost of a single $N$-point, 1-D FFT. Your proposed calculations are correct, then, if you assume $C_{1D}(N)=2N\log_2 N$, because this evaluates to $$2N_1N_2(\log_2 N_1+\log_2 N_2)=2N_1N_2\log_2(N_1N_2).$$

[EDIT: This paragraph is no longer relevant; the OP has changed the FFT sizes to a power of two.] However, your other problem is that $C(N)$ is not necessarily $kN\log_2 N$ (or close to it) unless $N$ is a power of two. In particular, if $N$ is prime, then you're going to need $kN^2$ operations. For values of $N$ that are composite prime factors, then the computational complexity is a function of the number of prime factors and their values. See this page for a discussion of efficient prime factor FFTs. Roughly speaking, if the number of prime factors is sufficiently large, then the cost is going to grow like $N\log N$.

Even for powers of two, $C_{1D}(N)=2N\log_2 N$ is not quite right for several reasons. First of all, the standard radix-2 FFT actually requires $(N/2)\log_2 N$ complex multiplications and $N\log_2 N$ complex additions. Furthermore, you really don't want to just add these two counts together. After all, one complex addition requires 2 real additions, and one complex multiplication requires 4 multiplications and 2 additions. So if you want to combine these two operation counts, you need to do it in terms of these real operations, giving you $C_{1D}=5N\log_2 N$ real operations.

You can do significantly better than this, it turns out, if you take a careful accounting of "trivial" complex multiplications. One of the most efficient algorithms from a real operation standpoint is the split-radix algorithm (which requires a power of two) requires $C_{1D}(N)=4N\log_2 N-6N+8$ real multiplications and additions for a single FFT. For a 2-D FFT of size $N_1\times N_2$, the total cost would be
$$\begin{aligned}
C_{2D}(N_1,N_2)&=N_1(4N_2\log_2 N_2-6N_2+8)+N_2(4N_1\log_2 N_1-6N_1+8)\\
&=4N_1N_2\log_2(N_1N_2)-12N_1N_2+8(N_1+N_2).
\end{aligned}$$

Finally, let me add that counting floating point operations is a crude way to measure practical performance. Memory access patterns are often a more important factor in performance. Quoting Wikipedia: "In practice, actual performance on modern computers is usually dominated by factors other than the speed of arithmetic operations and the analysis is a complicated subject (see, e.g., Frigo & Johnson, 2005), but the overall improvement from O(N2) to O(N log N) remains."


The S-transform allows you to deal with differential equations in an algebraic manner - so they become easier to solve. Since continuous/analog filters consist of integrators and differentiators the S-transform is therefore a natural way to deal with these systems.

The z-transform provides an algebraic way of dealing with finite difference systems and therefore it is a natural way to deal to discrete-time systems i.e. digital.

In the S-transform, setting $s=j\omega$ results in the Fourier transform. For continuous time systems we are interested in whether poles are on the right or left side of the $j\omega$ axis, because that determines the stability.

For the Z-transform, setting $z=e^{j\omega}$ results in the Discrete-Time Fourier transform. To determine stability, we are interested in whether the poles are inside or outside the unit circle.

So the S and Z domains are similar - they allows you solve continuous and discrete time systems, respectively,  using algebra.


Meaning of wavelet and scaling coefficients

Hope this helps:...

The fast wavelet transform (FWT) is a mathematical transform which results in a series of orthogonal coefficients of varying resolution.

Assume an input signal or time series $S$ of length $n=2^{j}$.  The length of $S$ must be equal to some power of 2, that is, $n$ equal to 8, 16, 32, 64, 128,  256, 512, 1024, etc.  To transform a signal into wavelets, at each resolution the signal is convolved using a series of scaling coefficients, $a_0$, $a_1$, $a_2$, $\ldots$ representing a low pass filter and a sequence of $d_0$, $d_1$, $d_2$, $\ldots$ of a high-pass filter.

Scaling function coefficients, $\phi$.      At the $j-1$th resolution, the convolution is taken over the signal's original length, $2^j$, in the form
\begin{equation}
\phi_k =  a_{-1}s_{2k+1} + a_0 s_{2k},
\end{equation}
where $k=0,1,2,\ldots,2^{j-1}-1$ are the indices of each scaling function, $s_0, s_1, \ldots, s_{2^{j}-1}$ are the original input signal elements, and $a_{-1}$ and $a_0$ are constants from the Haar scaling function.  An example of expanding the notation for the $j-1$ resolution is
\begin{equation}
\begin{split}
\phi_0 &= a_{-1}s_1 + a_0 s_0\\
\phi_1 &= a_{-1}s_3 + a_0 s_2\\
\vdots\\
\phi_{2^{j-1}-1} &= a_{-1}s_{(2^j-1)} + a_0 s_{(2^j-2)}.\\
\end{split}
\end{equation}
Consider a signal $S$ of length $n=1024=2^{10}$ $(j=10)$, and assume Haar scaling function values of $a_{-1}=a_0=0.5$.  Although the length of $S$ is 1024, the notation used to represent the series of values is $s_0$, $s_1$, $s_2$, $\ldots$, $s_{1023}$.   The first convolution of $S$ will result in 512 ($2^{10-1}=2^9$) scaling function coefficients, which are obtained using the relationship
\begin{equation}
\phi_k = 0.5 s_{2k+1} + 0.5 s_{2k}
\end{equation}
where $k=1,2,\ldots,511$.  This translates to
\begin{equation}
\begin{split}
\phi_0 &= 0.5 s_{1} + 0.5 s_{0}\\
\phi_1 &= 0.5 s_{3} + 0.5 s_{2}\\
\vdots\\
\phi_{511} &= 0.5 s_{1023} + 0.5 s_{1022}\\
\end{split}
\end{equation}
It is easily noticed that the scaling function coefficients are the average of each neighboring pair of signal elements.

Wavelet coefficients, $\psi$.  The wavelet coefficient is essentially based on the difference between each neighboring pair of signal elements.  At the $j-1$ resolution, these are
\begin{equation}
\begin{split}
\psi_0 &= -d_{-1}s_1 + d_0 s_0\\
\psi_1 &= -d_{-1}s_3 + d_0 s_2\\
\vdots\\
\psi_{2^{j-1}-1} &= -d_{-1}s_{(2^{j-1}-1)} + d_0 s_{(2^{j-1}-2)}\\
\end{split}
\end{equation}
where $d_{-1}$ and $d_0$ are both 0.5 based on the Haar wavelet.

For our first convolution of signal $S$ of length 1024 at resolution $j-1$, the 512 wavelet coefficients are
\begin{equation}
\begin{split}
\psi_0 &= -0.5 s_{1} + 0.5 s_{0}\\
\psi_1 &= -0.5 s_{3} + 0.5 s_{2}\\
\vdots\\
\psi_{511} &= -0.5 s_{1023} + 0.5 s_{1022}\\
\end{split}
\end{equation}


Zero Order Hold with polyphase resampler

I was able to achieve this with a square pulse impulse response of length I(interpolation rate). The square impulse response in relation to an equivalent sin(x)/x and linear interpolator(triangular) response is show in the figure below. This is still not a perfect design as I see discontinuities every I samples in the resampled output.

Edit: The discontinuities are addressed by changing the transition of high-low in the square impulse response, from two samples to three sample. So instead of the transitions being (high, low) in two samples, it will be (high, mean(high,low),low) over three sample and likewise for transition from low ot high.


Wiener Filter in Frequency Domain: What it does to a specific Frequency?

You already know the formula, so that's how the frequency response looks like. In frequency regions where the signal is stronger than the noise, the frequency response is close to $1$, so the portion of the input signal in that frequency range passes the filter with almost no attenuation, whereas in frequency regions where the noise is stronger than the signal, the filter's frequency response becomes small, and the corresponding part of the input signal is attenuated.

If you define a frequency dependent SNR

$$\textrm{SNR}(\omega)=\frac{\Phi_{ss}(\omega)}{\Phi_{nn}(\omega)}\tag{1}$$

you can rewrite the frequency response of the Wiener filter as

$$H(\omega)=\frac{\textrm{SNR}(\omega)}{1+\textrm{SNR}(\omega)}\tag{2}$$

Eq. $(2)$ clearly shows that if $\textrm{SNR}(\omega)\gg 1$ you have $H(\omega)\approx 1$ (i.e., no attenuation), whereas if $\textrm{SNR}(\omega)\ll 1$ you get $H(\omega)\approx \textrm{SNR}(\omega)\ll 1$ (i.e., strong attenuation).


Convolution of Signal with a Gaussian Filter / Kernel

If you have an implementation you want to validate the best way to do it is by validating it versus a reference implementation.

For instance, the above, in MATLAB, would be something like:

% Setting the Grid Parameters
leftBound = -10;
rightBound = 10;
numSamples = 10000;

% Signals Parameters
gaussianKernelStd   = 2;
gaussianKernelMean  = -1; %<! Basically its shift on the grid
sineAmp             = 0.1;
sineFreq            = 1;
sinePhase           = 0;

% Generating the grid
vX = linspace(leftBound, rightBound, numSamples + 1);
vX = vX(:);

samplingInterval    = mean(diff(vX));
samplingFrequency   = 1 / samplingInterval;

% Generating the Sine Signal
vS = sineAmp * sin((2 * pi * sineFreq * vX) + sinePhase);
% Generating the Gaussian Kernel
vK = exp(-((vX - gaussianKernelMean) .^ 2) / (2 * gaussianKernelStd * gaussianKernelStd));
vK = vK / sum(vK * samplingInterval); %<! Normalization

vY = conv(vS, vK, 'same');

figure();
% plot(vX, vK);
% figure();
% plot(vX, vS);
plot(vX, [vS, vK, vY], 'LineWidth', 2);
legend({['Sine Signal'], ['Gaussian Kernel'], ['Convolution Result']});

The result is given by:


Is Oversampling a Signal Same as Discretizing the Signal?

Oversampling is the case the rate the data is sampled is higher than required by the data Bandwidth and Nyquist Shannon Sampling Theorem.
It has many good properties regarding the processing yet it also generate some issues (Like coloring the quantization noise).

Regarding your questions:

In BPSK the symbol is the phase range of modulated signal (Sign). By using repmat() you indeed oversmapled the data (Symbols) yet not the modulated signal. What you needed to do is creating a grid where you sample the modulated signal. So if the the modulated signal has a frequency $ f $ than for 100 oversampling rate you should sample at $ {F}_{s} = 100 \cdot 2 \cdot f $.
What you did is replication of the symbol. Not higher sampling rate unless you imitate the effect of oversampling after Symbol inference.
In order to see the effect of oversmapling by factor of $ 100 $ generate the modulated signal over a grid where $ \Delta T = \frac{1}{100 \cdot 2 \cdot f} $ with $ f $ being the modulating signal frequency.

Regarding your question in the title.
The process of sampling data is basically discretizing data over 2 grids:

The Time Grid
Basically sampling the continuous data in specific times. Oversampling refers to the grid interval.
The Value Grid
Process we usually call quantization which means the value we sample isn't continuous but has discrete property of being a value withing a pre defined set of values (Quantization Grid).


You have the gist of it, but there are a couple of problems. I'm sorry, I don't have MATLAB, but I'll do my best to give you MATLAB-like pseudo-code to work with.

The main problem is that your Airy disk is an image, i.e. amplitude only (the result of taking abs). To reconstruct the circular aperture, you need phase as well. Here's code to first generate your complex Airy disk (PSF), from which you can then reconstruct the circular aperture (pupil function). Starting with this file, circ.jpg:

pupil = imread('circ.jpg');
pupil = rgb2gray(pupil);
% use fftshift to block-swap the fft output (move DC to the center)
complex_airy = fftshift(fft2(pupil));

% now, starting with this complex airy pattern, you were close:
pupil_reconstructed = ifft2(complex_airy);
pupil_reconstructed = abs(pupil_reconstructed);

Below are the amplitude and phase of the complex_airy pattern, obtained using abs(complex_airy) and angle(complex_airy), respectively. (I've included only the centers of the images, to show detail.) The problem with your initial approach is that it's an ill-formed inverse problem: a single Airy amplitude pattern can correspond to an infinite number of pupil functions. When you inverse transformed the amplitude only, you assume the phase is everywhere zero, which--as is clear from the phase image below--it's not. There are statistical approaches for estimating the pupil function from the Airy amplitude alone, using lots of other constraints gleaned from the light source properties and optical system, but I don't think that's what you were after.


Noise Removal from an Image Using OpenCV (Comparison to Neat Image)

NeatImage probably uses Wavelets based Noise Reduction.
You can look for methods based on that.

Today you need methods which are "Edge Aware", namely they smooth yet keep edges in tact.

Have a look here:
https://github.com/RoyiAvital/Fast-Anisotropic-Curvature-Preserving-Smoothing


The main difference between DCT and PCA (more precisely, representing a dataset in the basis formed by the eigenvectors of its correlation matrix - also known as the Karhunen Loeve Transform) is that the PCA must be defined with respect to a given dataset (from which the correlation matrix is estimated), while the DCT is "absolute" and is only defined by the input size. This makes the PCA an "adaptive" transform, while the DCT is data-independent.

One might wonder why the PCA is not used more often in image or audio compression, because of its adaptivity. There are two reasons :

Imagine an encoder computing a PCA of a dataset and encoding the coefficients. To reconstruct the dataset, the decoder will need not only the coefficients themselves, but also the transform matrix (it depends on the data, which it does not have access to!). The DCT or any other data-independent transform might be less efficient in removing statistical dependencies in the input data, but the transform matrix is known in advance by both the coder and decoder without the need for transmitting it. A "good enough" transform which requires little side information is sometimes better than an optimal transform which requires an extra load of side information...
Take a large collection of $N$ 8x8 tiles extracted from photos. Form a $N \times 64$ matrix with the luminosity of these tiles. Compute a PCA on this data, and plot the principal components that will be estimated. This is a very enlightening experiment! There is a very good chance that most of the higher-ranked eigenvectors will actually look like the kind of modulated sine-wave patterns of the DCT basis. This means that for a sufficiently large and generic set of image tiles, the DCT is a very good approximation of the eigenbasis. The same thing has also been verified for audio, where the eigenbasis for log- signal energy in mel-spaced frequency bands, estimated on a large volume of audio recordings, is close to the DCT basis (hence the use of DCT as a decorrelation transform when computing MFCC).


Where can I find information about how to implement IIR filters on micro-controllers?

Here are 3 points that you must consider :

1 - What is the minimal latency you can tolerate? This is important, as it's important to know whether to need process every sample as fast as possible or if you can perform batch processing which is typically more efficient but creates latency .

2 - No matter if you implement the IIR in fixed-point or floating-point, it is always a good idea to split the filter in sub-filters. This minimizes quantization issues that can occur even if you use double-precision arithmetic. One effective strategy is to cascade the filter in second-order sub-filters. Check out https://www.mathworks.com/help/signal/ref/sosfilt.html

3 - If you implement in fixed-point, it is a good idea to create a fixed-point model in Matlab/Python/Scilab/tool of your choice and compare it to the "ideal" floating-point filter before implementing it in your microcontroller.


Relating bin numbers to frequency non uniform fft

The magic of NUFFT is that you create a (oversampled) grid of equally spaced points by performing the interpolation and convolution of your unevenly spaced data. This is mostly done by using Gaussians. After that you apply the Fast Fourier Transform and minimise the effects of the filter you applied previously. Although this method was described by Dutt & Rokhlin in 1993. Currently you can find other, like ones based on Multiresolutional Analysis and Least Square Interpolation.

Nevertheless, you always end up with evenly spaced frequency bins - there are no unevenly spaced frequency bins due to jitter in the time domain.

If you wish to dig into math, then please refer to (among many others):

Particle Simulation Based on Nonequispaced FFT

NUFFT and MRI


Measuring length of DNA fibers from an image of single molecules

Example in Mathematica:

(* Get your image*)
img = Import["http://dl.dropbox.com/u/18072545/c_29.tif"];
(*Detect the extended minima and remove background*)
nB = ImageSubtract[img, ColorNegate@FillingTransform@ColorNegate[img]]

(*Separate RGB channels*)
cS = ImageAdjust /@ ColorSeparate[img]

(*Binarize*)
bcS = Binarize[#, .4] & /@ cS

(*Remove large elements*)
tH = TopHatTransform[#, DiskMatrix[2]] & /@ bcS

(*Detect lines using a Hough Transform*)
lines = ImageLines[#, .01, .8] & /@ tH
(*Plot them*)
Show[img, Graphics[{
   Thickness[.01], Red, Line /@ (lines[[1]]),
   Thickness[.006], Green, Line /@ (lines[[2]]),
   Thickness[.004], Blue, Line /@ (lines[[3]])}]]
(*Red and green are superimposed*)

Edit

Here you may see the Red and Green clusters apart. As you can imagine, you have to decide when a portion is red!


Beginner level : Help with terminologies : smooting length, length of channel, channel equalization delay, blind system identification, equalization

The term smoothing length (also called observation interval) is simply the length of the equalizer, i.e. its number of filter taps. In the code the filter length is L+1, so L is the filter order. Similarly, the length of the channel refers to the length of the FIR filter simulating the dispersive channel. Again, the channel filter length is actually ChL+1. What is called equalization delay in the code is the total delay introduced by the channel and by the equalizer. It is needed for aligning transmitted and received (and equalized) symbols.


I solved this by still using Otsu to find the threshold value. I then just filter the the regions with value == 1 in the thresholded signal (i.e. the region_mask; green line in plot below) by their width.
This is the result, where the rectangular peak on the right will be filtered (see below for the code):

plt.plot(signal_orig)
plt.plot(signal)
plt.plot(region_mask * np.max(signal))
plt.show()

Code:

signal = savgol_filter(signal_orig, 51, 1)  # window size 51, polynomial order 1
threshold = threshold_otsu(signal)

region_mask = signal > threshold
region_list = get_regions_from_mask(region_mask)
region_list = filter_date_regions_by_width(region_list, minimum_region_width=200)

def filter_date_regions_by_width(region_list, minimum_region_width):
    """
    Filters removes data regions from region_list with width < minimum_region_width.

    :param region_list: list of DataRegion objects, that we want to filter.
    :param minimum_region_width: minimum width of a region, so that it is not filtered out.
    :return:
    """
    filtered_list = []
    for region in region_list:
    if region.width >= minimum_region_width:
        filtered_list.append(region)
    return filtered_list

def get_regions_from_mask(region_mask):
    """
    Get list of the DataRegion objects that were found inside region_mask.

    :param region_mask: the np.array of 0 and 1 in which we want find regions of value == 1.
    :return: region_list: a list of DataRegion objects.
    """
    region_list = []  # list that will hold all regions that were found
    region = DataRegion()  # object to hold the start and end values of the region
    inside_region = False  # indicates, if we are currently inside a mask_region with value=1
    for index, value in enumerate(region_mask):
    if value == 1 and not inside_region:
        # entered a region
        inside_region = True
        region.start = index
    if value == 0 and inside_region:
        # left a region
        inside_region = False
        region.end = index
        region.width = region.end - region.start
        region_list.append(region)
        region = DataRegion()
    return region_list

Definition of the DataRegion class:

from dataclasses import dataclass

@dataclass
class DataRegion:
    start: int = None
    end: int = None
    width: int = None


Finding the output of a system where the input is a sum of complex exponentials

Okay... I managed to do it on my own. Here's the solution if you're interested.

We know that $$H\{x(t)\} = \sum_{k=0}^{\infty}a_k \  H\{e^{jw_kt}\}$$ since $H$ is a linear system.

We also know that $H\{e^{jw_kt}\}  = F(w_k) \ e^{jw_kt}$, where $F(w_k)$ is the frequency response of the system.

By definition, we have that $$F(w_k) = \int_{-\infty}^{+\infty}h(t) \ e^{-jw_kt} \ dt \\ = \frac{1}{\tau} \left(\int_{-\infty}^{+\infty}\delta(t+\tau)e^{-jw_kt} dt \ - \ \int_{-\infty}^{+\infty}\delta(t) e^{-jw_kt} dt\right)\\ =\frac{1}{\tau} \left( \int_{-\infty}^{+\infty} \delta(u) e^{\ jw_k(\tau-u)} \ du \ - \
\int_{-\infty}^{+\infty} \delta(t)e^{\ jw_k(0-t)}\ dt \right) \\ = \frac{e^{ \ jw_k\tau} - e^0}{\tau}$$

Therefore we can simply write $$H\{x(t)\} = \sum_{k=0}^{\infty} \frac{e^{\ jw_k\tau}-1}{\tau} \ a_ke^{\ jw_kt}$$


You could simply run your signal through some sort of high pass filter or shelf and compare the energy before and after the filter. If the ratio is high, you have a lot of high frequencies in your signal and should probably do less smoothing. If the ratio is low, there are not a lot of high frequencies and you can do more aggressive smoothing.

The frequency selectivity can be dialed in by the shape of the filter

As a simple example, the R code below runs two filters: a high pass one with unity gain at $f_s/2$ and a low pass one with unity gain at $0$.  Two test signals are used: one with more low pass content, one with more high pass content.

The output is:

[1] "X: 344.6108  LPF1: 32.1535  HPF1: 15.1717  LPF2: 17.4986  HPF2: 27.7821"

The low pass filtered test signal x_test_lpf retains more energy from the low pass filter than the high pass one (32.1535 vs 15.1717).  The high pass filtered test signal x_test_hpf retains more energy in the high pass output than the low pass one (17.4986 vs 27.7821).

R Code

#26668

x <- runif(1000,-1,1)

Nfilt <- 20

x_test_lpf <- c(x, filter(x,rep(1/5,5), circular = TRUE))
x_test_hpf <- c(x, filter(x,(-1)^(0:4)/5, circular = TRUE))

x_hpf_1 <- filter(x_test_lpf,(-1)^(0:(Nfilt-1))/Nfilt, circular = TRUE) # Unity gain at fs/2
x_lpf_1 <- filter(x_test_lpf,rep(1/Nfilt,Nfilt), circular = TRUE) # Unity gain at 0

x_hpf_2 <- filter(x_test_hpf,(-1)^(0:(Nfilt-1))/Nfilt, circular = TRUE) # Unity gain at fs/2
x_lpf_2 <- filter(x_test_hpf,rep(1/Nfilt,Nfilt), circular = TRUE) # Unity gain at 0

print(paste("X:", sum(x^2), " LPF1:", sum(x_lpf_1^2), " HPF1:", sum(x_hpf_1^2), " LPF2:", sum(x_lpf_2^2), " HPF2:", sum(x_hpf_2^2)))


Comparing an audio signal with a previous 'taught/profiled' signal

FFT won't work because it captures too many details which are rather irrelevant to the task. You need to extract a representation of sound which is more or less invariant by pitch and amplitude changes (properties of sound which are irrelevant to your problem). So you need a "blurred" version of the spectrum - and MFCCs are exactly that.

The textbook-ish, baseline system for speaker verification is to extract MFCCs (first coefficient excluded) on the "profile" recording; train a Gaussian Mixture Model on this data. The more data we have in the profile, the better (ideally, the most common phonemes should all be represented).

To score a new recording (during the "recognition" phase), extract the MFCCs, and compute the likelihood of the model given this data (you'll have to normalize by the length of the recording). Compare the score to a threshold - set according to the false positive rate you want to achieve.

For speaker recognition, you train one model per speaker to recognize; and select the speaker associated with the model with the higher score.


Causal LTI system having exponential input

Any LTI system with transfer function $H(s)$ reacts to an input signal $x(t)=e^{s_0t}$ with an output signal $y(t)=H(s_0)e^{s_0t}$, if $s_0$ is inside the region of convergence (ROC) of $H(s)$. This is equivalent to saying that the functions $e^{s_0t}$ are eigenfunctions of LTI systems, and $H(s_0)$ is the corresponding eigenvalue.

Note that $s_0$ is an arbitrary complex constant inside the ROC of $H(s)$, so the relation mentioned above is also valid for $s_0\in\mathbb{R}$. Also note that causality is not required for the above relation to hold.

However, if your question is about input signals of the form $x(t)=e^{at}u(t)$, where $u(t)$ is the unit step function, i.e., for input signals that are zero for $t<0$, then the output is generally not equal to $H(a)e^{at}u(t)$, because the functions $e^{at}u(t)$ are no eigenfunctions of LTI systems.


what is the output of BoW after an image has been trained with SIFT algorithm and k-means

You essentially got it right: the final purpose of the whole BoW clustering algorithm is to somehow produce a single image descriptor for every image.

In case of BoW clustering (either K-means, or hierarchical K-means, or some other clustering), this image descriptor is a histogram of visual words for that image often normalized by the number of local features for each image.

Let me offer a step-by-step explanation of the process:

first you get all the local features from all the images (or just from your learning dataset), and cluster all those features by a clustering algorithm of choice
after completing the clustering, for every cluster we calculate the "representative member" by some kind of mean or average over all the local feature samples from that cluster

This "representative member" is actually your visual word -- it's like the "core" of all the features of the cluster (if you would be working with actual words, the samples "work", "working" and "workaround" might all have a common representative member "work").

This step might actually be skipped -- we don't care about exact visual word, as long as we can distinguish between the clusters (e.g. knowing how to classify a new local feature in to "cluster 1", "cluster 2" or "cluster 3" is enough).
now, for every image in our database we construct the "image descriptor" by making the histogram by clusters for that image.

E.g., if there is 4 clusters in final clustering, and our image has 4 features in cluster 1, 10 in cluster 2, 4 in cluster 3 and 2 more in cluster 4, the descriptor might be: 4|10|4|2.
finally, we can normalize the histograms

Following with the same example from above, the image had a total of 20 features. We could thus normalize the 4|10|4|2 histogram by dividing everything by 20 and getting: 0.2|0.5|0.2|0.1.

With this, you ended your "database description" -- you have a numeric descriptor for all your images. Now, when a new image arrives, you do the following:

extract local SIFT (or any other kind) of features
determine for all the local features of the new image to which cluster they belong
make a histogram for the new image (+ normalize histogram)
compare the histogram of the new image to the histograms of database images

You have methods other than BoW to create the "image descriptors". One alternative to BoW is VLAD (Vector of locally aggregated descriptors), which is a bit more compact but has the same purpose.


looking for simplified version of motion capture for simplified situation where objects are solid and have no joint (no internal degree of freedom)

Fast way could be to convert your image to HSV format and then apply color thresholding as you have mentioned. Here is a short snippet which could get you started quickly.

import cv2 as cv2
import numpy as np

fn = 'image_or_videoframe'
# OpenCV reads image with BGR format
img = cv2.imread(fn)
# Convert to HSV format
img_hsv = cv2.cvtColor(img, cv2.COLOR_BGR2HSV)

# Choose the values based on the color on the point/mark
lower_red = np.array([0, 50, 50])
upper_red = np.array([10, 255, 255])
mask = cv2.inRange(img_hsv, lower_red, upper_red)

# Bitwise-AND mask and original image
masked_red = cv2.bitwise_and(img, img, mask=mask)

in this case, the red is filtered in the image and masked_red would contain only the red pixels in the image.

After the filtering, depending on the requirement you can use houghlines to detect the lines. It is not trivial but it is definitely doable.

Hope this helps and all the best :)


Please, tell what is not enough "layman".

Because standard linear smoothing (like averaging) tends to filter noise but reduces the contrast. Image composition is not often linear: cast shadows, occlusions. So linear methods tend to be suboptimal.
When you want to preserve sharp contrasts, for instance with multiplicative noise (whose intensity varies with image intensity), when you have a slow variation of the global illumination across the image (like progressive shadowing effect), when you want to both denoise and enhance
Images taken under varying illumination conditions, with noise varying with intensities, for instance related to temperature (infrared, visible of hot objects).


(python+opencv) haar cascade training using cloud computing

I created a google cloud platform compute engine. It gives you the option to choose the number of CPUs and ram size to choose from. So I experimented with a couple of combinations.

4 CPU 26gb ram - after starting the training I realised it used only ~7 GB of ram with all 4 CPUs running at 100%
8 CPU ~25gb ram - used around 9gb ram
24 CPU ~35gb ram. - used around 9gb ram.

All three instances ran at 100% CPU usage. With the increase in the CPU count, the faster the output files were generated.


YIN Pitch detection Algortithm ( how do I improve my results )

If you are not doing this in low-latency real-time, you can work backwards from the stable portion of the pitch estimate to the transient attack portion of the waveform at the beginning.

The sound of a plucked guitar string evolves in a possibly predictable pattern over time (e.g. more so than voice).  If you can estimate the onset time and/or have neighboring pitch estimates, you can adaptively set window sizes and threshold levels over time to more optimal values, as potentially determined by experimentation on some data set of guitar notes.  You can also use statistical decision theory to determine if any local pitch estimate fits the history of a reasonable spectral evolution of any guitar note, and reject outliers as noise, transients, harmonics or octave errors (harmonic and octave errors potentially being correctable errors).  This is especially useful working backwards in time, as the attack is usually noisier than the sustain/decay portion of a note's evolution.

Some improvement in steps 4) and 6) can also be acquired from psycho-acoustic experimentation with human listeners.  For instance: For an octave difference to be perceived above a certain error rate by typical human listeners, how much normalized difference in octave pitch estimate peaks over what amount of time is required?  Any difference over a smaller amount of time might be imperceptible.

ADDED: A window size of 1024 at a sample rate of 44100 only allows correlating/using a little more than 2 periods of the pitch of the lowest string of a guitar (E4 = 82.40Hz).  A 3X or 4X longer window might be more reliable for the lowest guitar notes, but shorter windows would be more responsive or provide better time locality for the higher fretted notes.


As Jason R has already mentioned in his comment, the noise reduction algorithm obviously introduces some amplification which renders the SNR estimates we proposed in our earlier answers useless. This becomes clear if we look at how the SNR is calculated.

Let $x_k$ be the speech signal without noise and let $n_k$ be the noise. Then the SNR before noise reduction is
$$
\gamma_{before} = \frac{\sum x_k^2}{\sum n_k^2}
$$

Now the noise reduction (NR) is applied to both signals yielding the pure speech signal after NR $\tilde{y}_k \neq x_k$.  Applying NR to the noisy signal produces some new signal $y_k + \tilde{n}_k$. As can be seen from your screenshots $y_k \neq \tilde{y}_k$. If we apply dspGurus method for calculating the SNR after NR we obtain
$$
\gamma_{after} = \frac{\sum \tilde{y}_k}{\sum [y_k - \tilde{y}_k + \tilde n_k]^2}
$$
The denominator not only conatins the residual noise but also an additional error $y_k - \tilde{y}_k$. This leads to a wrong estimation and a much lower SNR value.

To overcome this problem, you could try to normalize the signals. This will be difficult, because obiously the NR algorithm scales the noisy and the noise free signal differently. If you've implemented the algorithm yourself, you should try to fix this.

Another possible solution is to estimate the noise power in a region where the original signal contained silence and to calculate the total signal power in a region that contains noise and signal:

Take a region from the noisy signal after NR that only contains noise and calculate the mean power denoted as $\sigma_{n}^2$
Take a region from the noisy signal after NR that contains noise and speech and calculate the mean power denoted as $\sigma_{n,s}^2$

Now the SNR after NR can be estimated as
$$
\gamma_{after} = \frac{\sigma_{n,s}^2 - \sigma_{n}^2}{\sigma_{n}^2}
$$

How good this estimation is, depends on how good the following prerequisites are fulfilled in your case:

A large number of samples should be used to calculate the mean powers. You used 1401, that's not much.
Noise and signal after NR have to be uncorrelated. The NR algorithm possibly introduces noise that is a function of the signal. So this condition might not be fulfilled.
The mean noise power should be constant over time. I think this is fulfilled.


Continuous wavelet transform subsampling at lower frequencies?

The core of continuous wavelet transform

$$ C_s(\tau,a) = \int s(t) \psi^*_{\tau,a}(t) dt$$

discretization is that (this is only one of the numerous possibilities) one can find constants $a_0$ and $b_0$, whose product   $a_0 b_0 <C_\psi$ ($C_\psi$ depending on the wavelet), such that the following sampling pattern

$$c_{j,k} = C_s(k b_0 a_0^j,a_0^j), \quad (j,k)\in \mathbb{Z}$$

keeps all the information contained in $C_s(\tau,a)$. The subsampling at low frequencies is bigger, i.e. one need less samples since large wavelets smooth the signal.

The varying sampling space from low (top) to high (bottom) is illustrated below:


A couple of things jump out when looking at your diagram:

The first low-pass filter has a decimation of 4, so the sampling rate at its output is 256 kHz. However, you have set the sampling rate of the FM receiver and the second LPF to 224 kHz. These should be 256 kHz.
The cut-off frequency of the first LPF is 150 kHz. This means a total passband of 300 kHz. This is larger than necessary, so you'll end up picking up more noise than the optimum.
You multiply the audio output by a constant factor of 2. This may be too little, or too much. I'd suggest using a slider to control the volume. If you connect a scope sink to the audio sink input, you'll be able to easily get the optimum gain: you want the signal to be withing [1,-1].


There's a lot of literature on MFCCs on the web, so it would be a bit easier if you could be more specific as to which part of the processing you don't understand. But I'll give an overview of what needs to be done, hoping this is helpful for you:

compute the squared magnitudes of the FFT bins
weigh the bins using triangular windows; usually the windows are chosen such that the centers of the triangles are equidistant on a mel-frequency scale, and such that each triangle begins and ends at the centers of the two adjacent triangles. The mel-frequency scale is defined by
$$m=2595\log_{10}(1+f/700)$$ where $f$ is the frequency in Hz. Look at the figure to see how it works:

Note that the kernels are normalized such that the sum of the weights per triangle equals 1. Usually around 20 such triangular windows are used.
Take the logarithm of the weighted coefficients.
Compute the DCT.


Is a "Spectral Tilt Filter" the best way to get an arbitrary slope and cutoff LPF?

yeah, that's the way we've been doing pink noise filters (or "pinking filters") since the 80s.

if you ignore edge effects it's made by alternating these pole functions:

$$ \log|H_{2k}(j\omega)| = -\tfrac12 \log\left( 1 + (\tfrac{\omega}{p_k})^2 \right)  $$

and these zero functions:

$$ \log|H_{2k+1}(j\omega)| = \tfrac12 \log\left( 1 + (\tfrac{\omega}{q_k})^2 \right)  $$

and the total frequency response is

$$ \log|H(j\omega)| = \sum_n \log|H_{n}(j\omega)| $$

(this is not done yet)


How to deconvolve matrix when a model of noise exists?

This isn't (as far as I am aware) a standard approach to deconvolution, but it seems to address your problem:

Assume that your original signal $s$ is distorted by a circulant matrix $\mathbf{C}$ made up of the (shifted) components of your measured signal $m$:
$$
m = \mathbf{C} s
$$
Then just find the generalized inverse of $\mathbf{C}$, denoted $ \mathbf{C}^\dagger$ and form your estimate $\hat{s}$ as:
$$
\hat{s} = \mathbf{C}^\dagger m
$$

The image below shows:

$\color{black}{\tt black}$ :  the original signal in your question.
$\color{red}{\tt red}$ :  the distorted signal in your question.
$\color{green}{\tt green}$ :  100 realizations of inverting a noisy version of your distorted signal.

As you can see, the green estimates bounce around the original noiseless signal reasonably well.

R Code Below

# 31682

# http://stackoverflow.com/a/15796694/12570
circ<-function(x) {
  n<-length(x)
  matrix(x[matrix(1:n,n+1,n+1,byrow=T)[c(1,n:2),1:n]],n,n)
}

original <- c(0,0,0,1,0,0,0)
distorted <- c(0.001,0.005,0.1,0.5,0.2,0.001,0)

inverse_matrix <- ginv( circ(distorted[c(4,3,2,1,7,6,5)]) )

Nruns <- 100

output <- inverse_matrix %*% distorted

plot(original, type="l",
     ylim=c(min(c(original,distorted,noisy_distorted, output)),max(c(original,distorted,noisy_distorted, output))),
     xlim=c(1,length(output)), lwd=5)
lines(noisy_distorted, col="red", lwd=5)

for (runNo in 1:Nruns)
{
  noisy_distorted <- distorted + rnorm(length(distorted), 0, 0.01)
  output <- inverse_matrix %*% noisy_distorted
  lines(output, col="green")
}

lines(original, lwd=5)

title('Original, distorted, and estimated originals (100 realizations) ')


I don't see the question here.

But you're right.
On the DFT map of the image, the magnitude / amplitude, etc. are powers of the length of the vector of the complex number.

On more useful visualization to help working with large Dynamic Range is using the log(1 + abs()) on the DFT.


What Is the Difference between RLS, LMS and Wiener Filter? When Is One Preferred Over Another?

All three are Estimators / Predictors.
All of them try to estimate the coefficients of Linear Filter which minimizes an MMSE Cost Function.

The Wiener filter assumes all data is given and sets the way to calculate the optimal solution.

The LMS and RLS are sequential / on line methods to solve the same problem and given the data is stationary they all will converge to the same solution.

The LMS / RLS can adjust their solution according to data and hence are suitable in cases the model changes in time.

The easiest to implement with the least computing resources needed is the LMS and it is the one you should start with. The LMS works on the current state and the data which comes in.
The RLS, which is more computational intensive, works on all data gathered till now (Weighs it optimally) and basically a sequential way to solve the Wiener Filter.

So, I'd start with the LMS.
If the model is constant and the LMS Convergence rate is good enough for my problem I'd stick with it.
If it seems I need faster convergence rate then I'd go with the RLS.


Would the locations of SIFT features generally agree with features detected by Shi-Tomasi method?

Probably not. The SIFT detector finds centers of blob-like features. Shi-Tomasi detector finds corners. Furthermore, SIFT detector operates at multiple scales, while the classic Shi-Tomasi does not.


What is a good distance measure for matching SIFT descriptors depending on the distribution of their noise?

For obtaining distributions or modes in your, a basic approach is to cluster the data. That will of course bias the distribution towards the metric, or in other words the distribution will be tied to the metric. For example, it is possible to quantize the SIFT descriptors by using vocabulary trees to visual words and then plot the 1D histogram of the individual words. That's only one possible way to do this from a large collection of SIFT descriptors. Another approach would use nonparametric methods (Parzen Window, Meanshift...). In either case the distribution will be data-dependent and quite complex.

I think this question also deserves to mention the Root-Sift paper, where authors propose to use the Hellinger distance instead of Euclidean. The implementation is very easy (just take a sqrt) and the results are consistently improved. K-means can still be used. Apparently, Hellinger is a better metric than Euclidean for SIFT.


Should I consider rgb values of a pixel as one value?

For example if I have r=255, g=255, b=255, can I take pixel value as (in binary), (r<<16)+(g<<8)+b ?

While this number would identify the colour, it would not help with a mean or standard deviation. e.g. pure blue would be 255 but pure red would be 16711680

The paper is unclear so you have a few choices:

Process each channel individually
Convert the image to grayscale, find the mean and std dev, and use this for each channel.

There may be other methods that use different colour spaces.


estimation of ARMA(1,1) with MA parameter greater than 1

I suspect that the reason is that many ARMA fitting algorithms require that the ARMA system be stably invertible.  In your example, the MA process is not minimum phase as the zero will be outside the unit circle. That makes the inverted system unstable (the pole will be outside the unit circle).

I don't know the details of the sm.tsa.ARMA(y, (1, 1)).fit algorithm, but I suspect that it will also be so restricted.

See, for example, this lecture by Florian Pelgrin, which contains the slide below.

In particular, note the second condition.

This means that I suspect it will be harder to find an algorithm that does the estimation for the case you're  interested in.


Standard deviation of Gaussian corresponding to perfect low pass filter [1,1]

Suppose you want a low pass filter supported on $f \in [-1,1]$. The Fourier transform of a Gaussian is a Gaussian. So let's assume our Gaussian filter has a frequency response of the form $G_\sigma(f) = \frac{1}{\sigma \sqrt{2 \pi}} e^{-f^2/\sigma^2}$. For the best Gaussian fit in least squares sense, let us minimize the squared error ($\mathsf{SE}$):
$$
\mathsf{SE}(\sigma) = \int_{-1}^1 |G_\sigma(f)-1|^2 \;\;\; d\,f
$$
by choice of $\sigma$. This calls for some numerical evaluation.

I used Mathematica:

Integrate[ (1/s/Sqrt[2 Pi] Exp[-f^2/s^2] - 1)^2, {f, -1, 1}]

which gave me

2 - Sqrt[2] Erf[1/s] + 1/(2 Sqrt[2 Pi] s) Erf[ Sqrt[2]/s ]

Then plot the result as a function of $\sigma$ (i.e. s).

I see a minimum around $\sigma =0.7$.


Feature matching of images without corners

I have referred to Stage I.D of this tutorial. Hope this helps.
http://www.robots.ox.ac.uk/~vgg/practicals/instance-recognition/index.html#stage-id-improving-sift-matching-using-a-geometric-transformation

When the features have scale and orientation assigned (e.g. SIFT features have these properties), you can compute similarity transform between each corresponding feature pair. The similarity transform has 4 parameters: Tx, Ty - translation, R - rotation and S - scale.

The tentative correspondences are then filtered using RANSAC-like algorithm such that only correspondences consistent with common similarity transform (up to a predefined threshold) are used.


high & low pass correlated dither noise question

To achieve high-pass shaping , the “DC term” must be 0. Note that in the 2-dice high-pass case, every roll appears twice, once with a positive sign and once with a negative sign, which enforces the DC-free condition.  Any scheme involving more than 2 dice would still need to meet this condition.
The triangular PDF condition is satisfied by the fact that each sum or difference combines 2 rectangular PDF’s. I don’t know how you can get a triangular PDF by combining 3 rectangular PDF’s, regardless of how they are weighted.

Bob Adams


Frequency-domain zero padding - special treatment of X[N/2]

Let's look at the frequencies of the bins in your 8-point DFT:

$$
\begin{array}{c}
\omega_A = 0,\\
 \omega_B = \pi/4,\\
 \omega_C = \pi/2,\\
 \omega_D = 3\pi/4,\\
 \omega_E = \pi = -\pi\ ({\tt mod}\ 2\pi),\\
 \omega_F = 5\pi/4 = -3\pi/4\ ({\tt mod}\ 2\pi),\\
 \omega_G = 3\pi/2 = -\pi/2\ ({\tt mod}\ 2\pi), \\
 \omega_H = 7\pi/4 = -\pi/4\ ({\tt mod}\ 2\pi)
\end{array}
$$
So when you interpolate by a factor of 2, point $E$'s frequency becomes $-\pi$ or $+ \pi$.

At a first glance, I can't see what the problem is with your approach as it's not clear whether $E$ should be put into the bin associated with $\pi$ or $-\pi$.

On Julius O. Smith III's page, he states a condition:

  Furthermore, we require $x(N/2) = x(-N/2) = 0$  when $N$  is even, while odd  requires no
  such restriction.

And his example there is for an odd $N$, which avoids the problem.

Not sure it's required, but here's the full reference to Julius's work:

Smith, J.O. Mathematics of the Discrete Fourier Transform (DFT)
with Audio Applications, Second Edition,
http://ccrma.stanford.edu/~jos/mdft/, 2007, online book,
accessed September 28, 2011.


Is it only me or matlab's periodogram's function is confusing?

As the documentation states, periodogram provides a power spectral density estimate pxx:

[pxx, w] = periodogram(x);

meaning that it shows how the total variance of the signal, var(x), is distributed over the frequency w. This implies that the integral of the estimated spectral density over the frequency axis

trapz(w, pxx)

should correspond to that total variance.

This also means that the psd is not normalized, neither to a total variance of 1, nor to a maximum peak of 1 (= 0 dB). That the result is plotted (not computed!) in dB merely means that the vertical axis shows $10 \, \log_{10}($ pxx $)$, not that it shows a ratio. There is therefore no reason to expect a peak to align with 0 dB; depending on the range of values of the signal it can be both much higher and much lower.

To produce a plot from pxx and w, use

plot(w, 10 * log10(pxx))
xlabel('normalized frequency   [rad / sample]')
ylabel('power spectral density   [(rad/sample)^{-1}] in dB')

Reducing the level of the total signal by 3 dB means to divide the signal x by a factor of 2. If you do this numerically

periodogram(x)
hold all
periodogram(x / 2)

you will see that the second psd estimate has exactly the same shape as the first, it's just 3 dB lower. If this is not what you find with your measurement data, then this means that either your volume regulation or your ADC do not work correctly.

If you are really sure you want a normalized spectrum, do the following: First, compute the psd estimate:

[pxx, w] = periodogram(x);

Then, normalize the psd, e.g. with respect to the highest peak:

npxx = pxx / max(pxx);

You can then plot the normalized psd on a logarithmic scale like this:

plot(w, 10 * log10(npxx))
xlabel('normalized frequency   [rad / sample]')
ylabel('normalized power spectral density')

Note that npxx does not have a unit.

A note about units: The unit of numbers returned as pxx by periodogram & friends depends on the units of the input signal  x. Let's say these are dimensionless numbers; then the unit of the numbers in the pxx output is $\rm \frac{1}{rad / sample}$, because it is a density over w which is in units of $\rm rad / sample$. Or let's say these are voltage measurements coming from a calibrated ADC, so that the numbers in x are in volts. Then the unit of the numbers in the pxx output is $\rm \frac{V^2}{rad / sample}$.

If you additionally specify the sampling frequency in the call to periodogram:

[pxx, f] = periodogram(x, [], [], fs);

and fs is in Hz, the pxx is a density over frequency f, and its unit is $\rm \frac{V^2}{Hz}$.

Now if pxx is plotted, it is often plotted logarithmically, like stated above. If the function $10 \log_{10}$ is used, then it is said that the result is in decibel, dB. Then sometimes the unit of 10 * log10(pxx) is written as $\rm \frac{dB}{rad / sample}$. Strictly speaking, this expression doesn't make sense because dB is not just a unit (a different measurement stick) but denotes the use of a different quantity: the logarithm. That's why above I wrote the unit as [(rad/sample)^{-1}] in dB.


What phase rotation occurs when you take the derivative of an audio signal?

Simple Fourier Transform propery for differentiation is
$$\mathcal{F}\{\frac{\partial}{\partial t}x(t)\}=j\omega \cdot\mathcal{F}\{x(t)\}$$
Differentiation in time corresponds to multiplication with $j\omega$ in frequency. Hence the +6dB/octave slope. The phase shift is a constant 90 degrees for all frequencies.


What kind of power spectrum implies uncorrelation in time?

Since this does seem to answer the question, I'm reposting my comment as an answer:

I guess you may have to distinguish continuous (analog) and discrete (digital) signals. For analog signals, the autocorrelation function is a Dirac if and only if the spectrum is flat for all frequencies. Typically though what we care about is whether or not our digital samples are uncorrelated. It is easier to safisfy this. If your low-pass filter before your ADC is, say an ideal one, your ACF will be a sinc. If you get things right, your samples are in the zero crossings of the sinc and the discrete ACF looks like a discrete delta...

As the OP mentioned, the trick is to sample the signals in the zero crossings of the sinc function. For an ideal (rectangular) lowpass filter with (two-sided) bandwidth $B$ this means sampling exactly at a rate $1/B$.


Designing a 2d low pas filter in fourier domain for an image

It's been said so many times before that multiplying the Fourier transform coefficients of an image with a mask of ones and zeros in the frequency domain in order to achieve a short-cut to the ideal brickwall filtering in the time-domain is not a recommended method; and almost never preferred due to unavoidable, and unforeseen, errors it would introduce.

That being said, there are algorithms in computer geometry that describe how to approximate a cirle (or ellipse) on a discrete grid like that of a computer monitor.

The following is a circle algorithm that I remember having taken from wikipedia a decade ago. I have no idea on how it works, but just put it here for your sake.

% MIDPOINT CIRCLE ALGORITHM
Nx=256;          % image pixel size
Ny=256;
I=zeros(Ny,Nx);  % image onto which the circle will be drawn

x0 = 128;     % circle centers x-y
y0 = 128;
r = 50;       % circle radius

% ALGORITHM :
% -----------
f = 1-r;
ddfx = 1;
ddfy = -2*r;

x = 0;
y = r;

I(Ny-y0-r+1,x0) = 1;
I(Ny-y0+r+1,x0) = 1;
I(Ny-y0+1,x0+r) = 1;
I(Ny-y0+1,x0-r) = 1;

while x<y

    if f >= 0
       y = y-1;
       ddfy = ddfy + 2;
       f = f + ddfy;
    end%if

    x = x + 1;

    ddfx = ddfx + 2;
    f = f + ddfx;

    I(Ny-y0-y+1,x0 + x) = 1;
    I(Ny-y0-y+1,x0 - x) = 1;
    I(Ny-y0+y+1,x0 + x) = 1;
    I(Ny-y0+y+1,x0 - x) = 1;
    I(Ny-y0-x+1,x0 + y) = 1;
    I(Ny-y0-x+1,x0 - y) = 1;
    I(Ny-y0+x+1,x0 + y) = 1;
    I(Ny-y0+x+1,x0 - y) = 1;

end%for

figure,imshow(I)

The output will be your mask :

I think you can modify it to work for ellipse and some fill in effects...


Assuming that you understand the left-hand side of Eq. $(8.47)$, for understanding the right-hand side you need to know that $-1=e^{j\pi}$, and that $e^{j2k\pi}=1$. So in order to obtain all $2N$ roots of $(-1)^{\frac{1}{2N}}$ you rewrite $-1$ as
$$-1=e^{j\pi}e^{j2\pi k}\tag{1}$$
from which you get
$$(-1)^{\frac{1}{2N}}=e^{j\frac{\pi}{2N}}e^{j\frac{2\pi k}{2N}}=e^{j\frac{\pi}{2N}(1+2k)},\qquad k=0,1,\ldots,2N-1\tag{2}$$
And since $j=e^{j\pi/2}$, the final result is
$$\Omega_cj(-1)^{\frac{1}{2N}}=\Omega_ce^{j\pi/2}e^{j\frac{\pi}{2N}(1+2k)}=\Omega_ce^{j\frac{\pi}{2N}(1+2k+N)},\qquad k=0,1,\ldots,2N-1\tag{3}$$
If this is all new to you, you should really read up on complex numbers. They are fundamental in understanding many concepts in signal processing.

PS: There is a typo in the equations in the book. The numerator on the right-hand of Eq. $(8.46)$ should have $\Omega_c$ instead of $\Omega$. Also, the left-hand side of Eq. $(8.47)$ should have $\Omega_c$ instead of $\Omega$.


Create an input so that poles show up as outputs?

You should always be able to achieve this after $N-1$ samples. For this I will denote the given transfer function as

$$
H(z) = \frac{A(z)}{\prod_{k=1}^N (z - p_k)}. \tag{1}
$$

Now by starting with an impulse and filter it with

$$
F(z) = \frac{\prod_{k\neq i} (z - p_k)}{z^{N-1}} = \sum_{k=0}^{N-1} \alpha_k\,z^{-k} \tag{2}
$$

and use the output of this filter as the input to the system. The combined transfer function looks like

$$
H(z)\,F(z) = \frac{A(z)}{(z - p_i)\,z^{N-1}}. \tag{3}
$$

Therefore this filtered signal should only excite the mode of the system associated with the pole $p_i$. And the $z^{N-1}$ term only adds a delay of $N-1$ samples, so it only takes a finite amount of time after which only $p_i^n$ is visible in the output.

However it can be noted that this does require you to have perfect knowledge of the poles you cancel.


Information i can take from coeficient intensity map in DCT transform

The DCT is used in many places for its "energy compaction" properties. This is mostly useful in compression. For example, if I have a signal and I take its wavelet transform, then I can "compress" the signal by throwing away certain coefficients. These coefficients are not known before-hand. I have to look at the wavelet coefficients and then throw away all coefficients below a certain threshold to achieve some amount of compression with an acceptable loss in reconstruction accuracy.

For the DCT however, most of the energy of a signal is collected in its "low" frequency coefficients. What that means is that if you only keep the first lets say N coefficients (you don't need to store the locations of which coefficients you kept), you can get decent compression. In 2D, a "zonal" or "zig-zag" pattern is used to choose the coefficients to keep.

On a more fundamental level, the DCT is basically a linear combination of your signal as a bunch of discrete cosine waves. If you think of each cosine wave as an elementary building block, then the DCT tells you how much contribution to your image comes from a discrete cosine of a particular frequency.


I have a partial solution, without filtering the black spots.

These are my suggestions. They are pretty low tech, so feel free to play around with them.
First get rid of the uneven illumination using a Bandpass Filter: These are the specific operations that I applied to your image recorded as a FIJI / ImageJ macro:

//convert to 32 bit image
//manually by: Image->Type->32 bit
run("32-bit");
//Apply bandpass Filter:
//Manually by: Process->FFT->Bandpass Filter...
run("Bandpass Filter...", "filter_large=75 filter_small=0 suppress=None tolerance=0");
//Enhance Contrast
run("Enhance Contrast", "saturated=0.35");

What that does is get rid of the slowly varying background in your image. Now you can manually adjust the threshold using Image->Adjust->Threshold. Press Apply after adjusting the threshold to your needs. This is what it looks like:

Again, this is all very low tech, but try to play around with this. Finally you might want to run Process->Binary->Skeletonize on the image.


How to find correlation between two signals regarding peaks?

Your premise is due to the physical nature of the device under test it is possible that occurrence of higher frequencies vibrations can be predicted by lower frequencies and this is important as it will allow you to use a lower bandwidth sensor that can only respond to lower frequencies.

Interesting problem and I would suggest this approach to find evidence of any such correlation:

You will be simply correlating “waveforms” to determine this and in this case the suggested waveform to correlate would be frequency versus time. You can establish these waveforms using a sliding windowed DFT. Decide on a window length (further details on that below) and starting at the beginning of your much longer data captures window a block of the data (more details below) and then use the FFT to compute the DFT for that block. Save the data as time index n=0 and then shift the window by one sample and repeat to get the data to time index n=1. Repeat this to the end of the data set such that the result for each DFT bin represents frequency versus time (a form of the STFT).  Once you have this result with frequency bins as columns and time index as rows you can cross correlate all columns to assess if there is anything significant. Anything that does not stand out significantly represents a noise sample, and you can rss all noise samples to get the noise floor of your measurement approach. I would say that anything that is <4 sigma is not significant and anything <2 sigma would reasonably be noise only. You can also do a histogram of the results which should be a standard bell curve but will also help show significant outliers.

Considerations for window size and window choice: with the length of the window you are trading frequency resolution with time resolution: the wider the window the narrower will be the frequency resolution. Specifically if you used a rectangular window the equivalent noise bandwidth of the result is 1 bin wide so the resolution in this case is $f_s/N$ where N is the number of samples or bins in the DFT and $f_s$ is the sampling rate. However the result will be the average frequencies over the time duration of the window so the trade space between time and frequency should be clear.

Next comes the window choice. With a rectangular window (no further scaling of the data prior to taking the DFT) the equivalent noise bandwidth is one bin as mentioned above. This is actually the best resolution you can achieve with the approach I describe as any other window will make the resolution wider. However most importantly is that with the rectangular window the shape of the equivalent filter in frequency is a Sinc function centered on each bin with first null at 1/T or 1 bin away. This means that there will be significant spectral leakage from other frequencies resulting in a false “correlation” if further windowing is not used. Windowing provides a trade of frequency resolution with dynamic range. The rectangular window has the best resolution (narrowest main lobe in its frequency response) but poorest dynamic range (highest sidelobes in its frequency response) compared to all other windows. It is important given your goal to window so that you don’t get a false sense of correlation due to those high sidelobes. I recommend the Kaiser window since you can adjust its Beta term to dynamically make the trade of resolution and dynamic range.

My guess is this would show that there is no significant correlation beyond integer harmonics but without having the mechanical model I can’t be certain of that so this will be interesting to see. If you do proceed with this,please post your results of what you find.


Why do signal have mirror signal after pure frequency shift?

Your code is (probably) ok. But your way of plotting is slightly wrong.

The complete spectrum of the baseband PSK signal would include two peaks one at positive frequency $f_0=1000$ and the other at negative frequency $f=-1000$. When you IQ modulate this signal with $f_m=3000$ you will create two twins one at $f_{11}=2000$ and $f_{12} = 4000$ and the other at the negative frequencies at $f_{21}=-2000$ and $f_{22} = -4000$.

Since you are taking only the first half of the FFT you do not observe the peaks at negative ones, which causes the confusion. If you plot the full spectrums, you can see the missing components at the second half that you omit...

t = 0:0.000001:0.1;
fs = 1/(t(2) - t(1));
%ss = exp(j*2*pi*100*t) + exp(j*2*pi*200*t) + exp(j*2*pi*300*t) + exp(j*2*pi*900*t) + exp(j*2*pi*600*t); %Generating I-Q signal *.bin file for AWG

N=8;
bit_stream = round(rand(1,N));
f = 1000;
P1= 0;
P2= pi;
PSK_signal = [];

for ii = 1:1:N
PSK_signal = [PSK_signal (bit_stream(ii)==0)*sin(2*pi*f*t + P1)+...
        (bit_stream(ii)==1)*sin(2*pi*f*t + P2)];
end

ss = PSK_signal;
%figure(1); plot(ss); title('The PSK baseband signal');

ff = linspace(-fs/2,fs/2,length(ss));
figure,plot(ff,abs(fftshift(fft(ss))));
title('The PSK baseband SPECTRUM ');
% sp = fft(ss);
% N = length(ss);
% fp = ((0:N-1)/(N-1))*fs;
% fp = fp(1:N/2);
% sp = abs(sp(1:N/2))/(N/2);
% figure(2); plot(fp,sp); %Spectre of primordial signal
                     %title('The PSK baseband SPECTRUM ');

w = 3000*pi*2;
t = (1:length(ss))/fs;

%Real outputs from AWG
I = real(ss); % I Input on FPGA
Q = imag(ss); % Q Input on FPGA

%Adding frequency SHIFT in Simuation FPGA plate
ss_i_real = I.*cos(w*t) - Q.*sin(w*t); % I component from FPGA
ss_i_imag = I.*sin(w*t) + Q.*cos(w*t); % Q component from FPGA

sp = fft(ss_i_real+1j*ss_i_imag);  % computing the Spectrum of the ANALYTIC signal here! Hence the NEGATIVE frequencies are deleted.
figure,plot(ff,abs(fftshift(sp)));
title('IQ modulated signal FSK baseband');

% sp = fft(ss_i_real+1j*ss_i_imag);
% N = length(ss_i_real);
% fp = ((0:N-1)/(N-1))*fs;
% fp = fp(1:N/2);
% sp = abs(sp(1:N/2))/(N/2);
%figure(2);  plot(fp,sp);
%figure(3),plot(fp,sp); title('IQ modulated signal FSK baseband');

w = 40000*pi*2;
t = (1:length(ss_i_real))/fs;

%ss_c = ss_i.*exp(1j*w*t); %Modulation that carrying signal from FPGA on Carrier Frequency (40KHz)
ss_c_real = ss_i_real.*cos(w*t) - ss_i_imag.*sin(w*t);
ss_c_imag = ss_i_real.*sin(w*t) + ss_i_imag.*cos(w*t);

ss_c = ss_c_real+ss_c_imag*1j; % Real part of Signal onto carrier <- this is the signal that will be transmitting from modulator

sp = fft(ss_c);
N = length(ss_c);
fp = ((0:N-1)/(N-1))*fs;
fp = fp(1:N/2);
sp = abs(sp(1:N/2))/(N/2);
figure(4); plot(fp,sp); %Spectre of signal onto carrier freq
title('IQ modulated and up-mixed FSK signal');

% figure(4); plot(ss_c_real); %I - component of signal onto carrier freq
%
% figure(5); plot(ss_c_imag); %Q - component of signal onto carrier freq


why do we use $X(e^{j\omega})$ instead of $X(j\omega) $ in Discrete Time FT

The argument $e^{j\omega}$ emphasizes the $2\pi$-periodicity of the discrete-time Fourier transform (DTFT) of a sequence. Furthermore, if the $\mathcal{Z}$-transform $X(z)$ of the sequence $x[n]$ exists and if the unit circle $|z|=1$ is inside the region of convergence, then the DTFT of the sequence is simply obtained by evaluating $X(z)$ on the unit circle $z=e^{j\omega}$.

That's the same thing as with the continuous Fourier transform and the Laplace transform, referred to in your previous question.


Sample and reconstruct a real exponential (just one period)

First, since $t>0$ in your case, you can write your function as

\begin{equation}
C(t)=1.6925\left(\exp^{-0.136t}- \exp^{-1.192t} \right) u(t)
\end{equation}
where $u(t)$ is a unit-step function. Then, denote the FT of $C(t)$ as $C(F)$, which is given by
\begin{equation}
C(f)=1.6925\left(\frac{1}{0.136+2\pi jf}- \frac{1}{1.192+ 2\pi jf} \right)
\end{equation}

The spectrum $C(f)$ is plotted below. We can see that $C(f)$ is almost zero after 5 Hz. You can choose the sampling frequency according to the Nyquist criterion as $f_{\max} \geq 10$ Hz. You have taken a very high sampling rate (1000 Hz) that is 100 times more than the required rate. When I took 10 Hz to sample the signal above and took the FFT I got the results as follows  (in order time domain $C(t)$, frequency domain signal $C(f)$, frequency domain domain signal obtained via FFT, that is, FFT(C(t)). Note the inherent periodicity of the FFT spectrum)

And here is the MATLAB/Octave code for your reference..

clear all; clc;close all;

%% Time-domain signal
Ts=0.1; % Sampling time (obtained from the frequency domain plot)
t=0:Ts:60; % in sec
C=1.6925* (exp(-0.136*t)-exp(-1.192*t));
figure,plot(t,C);
title('Time domain signal');
xlabel('Time in sec')

%Frequency domain signal
f=0:.1:10;
F1=1.6925*((1./(0.136+(j*2*pi*f)))-(1./(1.192+(j*2*pi*f))));
figure,plot(f,abs(F1));
title('Frequency domain signal');
xlabel('Frequency in Hz')
% Via FFT
F2=fft(C);
L=length(F2);
figure,plot((0:L-1)*(1/(L*Ts)),abs(F2));
title('Frequency domain signal via FFT of C(t)');
xlabel('Frequency in Hz')


How do I go from LPC coefficients to a filter polynomial?

The prediction polynomial is the representation in the $z$-domain of the first equation you wrote.

When performing LPC, one assumes that the filter consists of only poles and the transfer function is $\frac{1}{A(z)}$ where

$$A(z)=1 - \sum_{i=1}^{M} a_i z^{-i}$$

That's the polynomial you are talking about. Note that it should have only negative powers of $z$, not positive.

As you wrote in your question, one estimates the signal $h(n)$ based on past values of that same signal. The $z^{-i}$ factors in the polynomial represent those delays in the $z$-domain, while the coefficients $a_i$ are just constants and, due to the $z$-transform linearity, they just stay as they are.


if this is a first-order HPF DC blocking filter, this is the recommendation i have for it in fixed-point, and i think it would be a good idea in floating point, too.

while (n<something)
  {
  y[n]   =  x[n]  -  w[n]*( 2^(1-N) );
  w[n+1] =  w[n]  +  y[n]*( (2^(N-1))*(1-p) );
  n = n+1;
  }

where p is the pole (in the z-plane and inside the unit circle) and $0<1-$p$\ll 1$.  N is a big integer.  with fixed-point you would perform the w[n]*(2^(1-N)) multiplication with arithmetic shift right (ASR) or the (w[n]>>(N-1)) operation.  in float you can just multiply it and that simply changes the exponent in the floating point value.  the constant coefficient (2^(N-1))*(1-p) should be precomputed and N chosen so that this is between $\frac12$ and $1$.

the -3 dB corner frequency $\omega_0 = 2\pi\frac{f_0}{f_\mathrm{s}}$ (as normalized angular frequency) is a function of the pole p and is:

$$\begin{align}
\cos(\omega_0) &= 1 - \tfrac{(1-p)^2}{2p} \\
 &= 2 - \tfrac12 \left(p + \tfrac1p \right) \\
 &= 2 - \cosh(\ln p)
\end{align}$$

the pole p is

$$\begin{align}
p \ &= \ e^{\operatorname{arccosh}(2-\cos(\omega_0))} \\
\\
 &= \ 2 - \cos(\omega_0) - \sqrt{\big(2 - \cos(\omega_0)\big)^2 - 1} \\
\\
 &= \ 1 - 2\sin\left(\tfrac{\omega_0}{2}\right) \left(\sqrt{1 + \sin^2\left(\tfrac{\omega_0}{2}\right)} - \sin\left(\tfrac{\omega_0}{2}\right) \right) \\
\end{align}$$

so
$$ 0 \ < \ 1-p \ = \ 2\sin\left(\tfrac{\omega_0}{2}\right) \left(\sqrt{1 + \sin^2\left(\tfrac{\omega_0}{2}\right)} - \sin\left(\tfrac{\omega_0}{2}\right)\right) \ \ll \ 1 $$

and the shift or exponent change N is an integer and should be about

$$ N \approx 1-\log_2(1-p) \approx 1 - \log_2(\omega_0) $$

someone can check my math.

this is an improvement by Tim Wescott to my DSP guru trick.  it deals with both the cosine problem and the limitations of fixed point arithmetic.

latest thinking:  i am not sure, but this might work fine for floating-point operation x, y, and w are all single-precision floats (double-precision floats should be no prob).  (1-p) is a coefficient calculated in advance and is very close to zero.

while (n<something)
  {
  y[n]   =  x[n]  -  w[n]*(1-p);
  w[n+1] =  w[n]  +  y[n];
  n = n+1;
  }


Finding the amplidude of a single frequency in signal

Goertzel in fact is the matched filter for a single frequency – so, it's, in the presence of uncorrelated noise, the estimator that gives the best estimation variance under fixed observation.

But: for Goertzel to work, you need to use a number of samples that is

an integer and
a multiple of $f_\text{sample}/83\,\text{Hz}$.

Unless your sampling rate hence is a small multiple of 83 Hz, that's usually a bad deal.

You can remedy that with a number of things: often, a narrowband (often: even resonant) filter can filter out a frequency range, but you'll need to deal with the skirt of that filter. That skirting will be different than that of Goertzel's algorithm (unless you happen to have also windowed the signal before applying the Goertzel, and now use the same Window as filter design prototype; sounds unlikely).

If you know more about the signal than you're telling us (it's result of AR processes, so you can apply Yule-Walker to get an appropriate model of what's happening, for example. Or, you know the signal is the sum of white noise and a known number of oscillations – then, autocorrelation subspace methods such as MUSIC are very efficient), then you might have methods that are very exact, but use even fewer samples than minimal Goertzel.

Anyway, the beauty of Goertzel really lies in the fact that it's pretty efficient to begin with. Honestly, if your frequency of interest is 83 Hz, and you still can't detect it in realtime even on the weakest of microcontrollers, then you're doing something wrong, with 99.5% probability¹. You really might want to ask a different question on here, describing what you overall want to do (where does your signal come from? What's at 83 Hz? Why do you want to determine the amplitue of that? How exactly do you know that frequency? What's the amplitude estimate reliability that you really need? So many things to talk about...).

¹ This number, like 22.36% of all study results, is totally made up.


How to properly time-shift a function in relation to another function

See that if you are going to go from $y(t)$ to $y(t-\tau)$ then you should replace every $t$ in the first one with $t-\tau$. This means that, in your case,

$$y(t-\tau) = x(t-\tau-2) + x(2-(t+\tau)) = x(t-\tau-2) + x(2-t+\tau)$$

In other words,

$$y(t) = x(-t)\implies y(t - \tau) = x(-(t - \tau))= x(-t + \tau)$$


Downsampling/Decimation - When is well done?

Ignoring your original question, and referring to your comment:

  I am observing the FFT of the downsampled signal presents an attenuation. I don't understand why.

The output of FFT (or DFT in general) isn't normalized - it is relative to the signal length (the total amount of samples being transformed). When decimating a signal (say, by a factor of M), you end up with a signal with 1/M length the original signal's length.

So, because DFT/FFT results are relative to signal length, a decimated signal with a 1/M length would have DFT/FFT components "attenuated" by a factor of M as well, even if the original signal was properly downsampled (i.e. bandlimited before being decimated).

To properly compare between the DFT of different signals, you should divide the output of DFT by the length of the transformed signal, to normalize the transform.


How to compute the power bands of an eeg signal using python?

Here is some code that may solve your problem:

from scipy.io import loadmat
import scipy
import numpy as np
from pylab import *
import matplotlib.pyplot as plt

eeg = loadmat("mydata.mat");
eeg1=eeg['eeg1'][0]
fs = eeg['fs'][0][0]
fft1 = scipy.fft(eeg1)
f = np.linspace (0,fs,len(eeg1), endpoint=False)
plt.figure(1)
plt.plot (f, abs (fft1))
plt.title ('Magnitude spectrum of the signal')
plt.xlabel ('Frequency (Hz)')
show()

You can also check this other link:

http://forrestbao.blogspot.pt/2009/10/eeg-signal-processing-in-python-and.html


Find the Differential Equation of an Analog Band Pass Filter

I'm afraid your first equation is wrong because it assumes a highpass filter without any load. However, there is a lowpass filter connected at its output and for this reason the current through the resistor leading to ground does not equal the current through $C_1$.

What you need to do is apply the complex voltage divider rule in the frequency domain, which will give you

$$\frac{V_{io}(\omega)}{V_{in}(\omega)}=\frac{Z(\omega)}{1/j\omega C_1+Z(\omega)}\tag{1}$$

where $Z(\omega)$ is the impedance of $R$ in parallel with the series connection of $R$ and $C_2$. I'm sure you can figure this out. The relation between $V_{out}$ and $V_{io}$ is simply

$$\frac{V_{out}(\omega)}{V_{io}(\omega)}=\frac{1/j\omega C_2}{1/j\omega C_2+R}=\frac{1}{1+j\omega RC_2}\tag{2}$$

The total response $V_{out}/V_{in}$ is simply obtained by multiplying Eq. (1) with Eq. (2). The time domain differential equations are obtained by multiplying Eqs (1) and (2) by the denominators of both sides and noting that a multiplication by $j\omega$ corresponds to differentiation.


Allan Variance vs Autocorrelation - Advantages

I currently work with the inner design details of Atomic Clocks so use the Allan Variance and Allan Deviation (ADEV) extensively. The primary point is that it can be used for non-stationary processes (which frequency noise is). For non-stationary signals where the autocorrelation or power spectral density can't be used with consistency, the Allan Deviation will converge to a consistent metric for comparison of the frequency accuracy of two different clocks (and generally for the comparison of the noise in non-stationary signals).
Non-stationary signals will often have a shorter time duration over which the signal can be sufficiently assumed to be stationary. Oscillator phase noise is an excellent example of this: for larger frequency offsets (typically 1 Hz and above for lower frequency oscillators such as 10MHz) a power spectral density (PSD) can be computed and given as the oscillator’s phase noise performance. The issue is this PSD cannot be computed consistently at lower frequency offsets as the non—stationary contributions start to dominate and this is where the ADEV computation will shine: in computing the longer term frequency accuracy of the clock over time durations where non-stationary contributions dominate (1/f noise, drift). ADEV is computed as an rms error given by the difference between the average frequency accuracy over interval tau seconds compared to the average frequency accuracy tau seconds ago. The rms of this error is computed over a very long duration (this error difference is a stationary signal and thus the standard deviation can be consistently computed).
One specific utility for me that the ADEV has provided well beyond the world of clocks is to determine the time constraints under which a random process can be sufficiently assumed to be stationary: if the random process is stationary we can continue to average to get a better estimate of the underlying mean— therefore we can use the ADEV to determine the optimum averaging time after which no further improvement can be gained (and through further averaging our estimate can get worst). This applies to everything from channel estimation to financial markets, wherever we may be observing past results to make an estimate of an underlying process.
I explain further details of the Allan Deviation and its potential use and application well beyond oscillator stability in these other posts:
How to interpret Allan Deviation plot for gyroscope?
What determines the accuracy of the phase result in a DFT bin?
Averaging data from 2 sensors
Is it possible to add some kind of noise to a non stationary (having unit root) signals to make them stationary?


Determine a sequence whose N point DFT is given without computing inverse

Hints only:

Think about time-frequency duality of DFT pairs. What is the DFT of sequence $-1^k$?

And then there is a phase term in the DFT, so there should be a shift in discrete time sequence. Use time-delay property of DFT.

These 2 properties should give you the solution.


Intuitive explanation of the Fourier Transform for some of the functions

What might seem intuitive differs greatly between individuals.  But let's start with a few basic things about a Fourier transform that people who have studied it might know about it.  Intuitively or not.

Basic concept 1:
Something symmetric around t=0 in the time domain is strictly real in the frequency domain (as only cosine functions and DC are purely symmetric in the time domain).

Basic concept 2:
The duality principle of the FT implies that if any transform looks intuitive in one direction, the transform in the opposite direction should also.

Transform 1:
A DC level in the time domain corresponds to an impulse at 0 in the frequency domain.

Transform !B:
Getting away from impulses and infinities, but staying with the same picture, a strictly real and really skinny Gaussian pulse in the frequency domain at DC must convert to a really really wide Gaussian envelope in the time domain, also symmetric around t=0.  And vice versa.  Going in between those two (and depending on your graph scaling), a medium width Gaussian in one domain should correspond to a medium width Gaussian in its transform domain.

Possible intuitive transform concept 2 (if you have studied signal modulation concepts):
Adding strictly real and symmetric sidebands to a carrier in the frequency domain corresponds to AM modulation with a cosine sinusoid in the time domain.  The farther away the sidebands from the carrier center, the higher the frequency of the modulation.  The closer the sidebands, the lower the frequency of the modulating sinusoid.

Thus if we add symmetric wide sidebands to a wide Gaussian, that will correspond to a narrow Gaussian that is high frequency AM modulated in the other domain.  If we add only narrow sidebands close in to a narrow Gaussian, that will correspond to a really wide Gaussian with a low frequency AM modulation in the other domain.

Possible intuitive transform concept 3 (if you have done much filter design):
If you have had to design  FIR filters, you might know that the sharper the cut-off and the flatter the pass-band desired in the frequency domain, the longer the FIR filter becomes to meet those requirements, corresponding to a longer and wider impulse response in the time domain.  Or, in the other direction, the sharper-edged and flatter-topped a pulse becomes in the time domain, the wider the bulk of its frequency response spans in the frequency domain.

So, let's take that Gaussian pulse, and add some similar shaped side-bands within the bulk of its main lobe so that its top become flatter. The Gaussian in the other domain will get amplitude (AM) modulated with sinusoidal modulators, mostly of frequencies where the tops of the added sidebands are much higher that the original Gaussian roll-off. (To make it even flatter, perhaps remove the some or all of carrier (original Gaussian) to produce a modulation over 100%, e.g. that crosses zero.)

Let's also chop off the tails of the Gaussian with a brick wall filter and increase the flatness of the pass-band top even more to “square it up” in one domain.  In the other domain, that will cause lots of added ripples extending wider than before, or a “fattening up” of the tails of the original Gaussian.

So as we “square up” a Gaussian in one domain, it gets amplitude modulated plus fatter tails in the other domain.

Thus, a skinny “squared up” Gaussian pulse in one domain might corresponds to a wide slowly AC modulated pulse with fat tails in the other domain,  And a wide but “squared up” Gaussian might correspond to something like a narrow pulse with fast AC modulation, but also with fat tails, in the other domain.

An AC modulated hump with fat tails starts to approach looking more similar to a Sinc waveform than the original Gaussian waveform, and with a width inversely proportional to the width of the "squared up" looking waveform in the other domain.

Intuitive?  Who knows…


how to stop matlab from rounding my values and keep them at same precision during calculations?

The precision is kept. It is just the display values that are truncated.

Type format long to change the displayed precision.


This should be impossible in most cases due to different modulation and coding schemes.

For example, DVB-S2 uses single carrier phase modulation wheras DVB-C2 uses OFDM modulation. There is no way to convert single carrier to OFDM without going to the digital domain.


Is there just a carrier in each frequency band in FDMA?

That's undefined.

FDMA just says that different users use different parts of the spectrum.

What these users do within "their" frequency ranges is up to the system to define.


How to plot the discrete time signal $u[2k - 3]$?

$u[n]$ is the unit step function with the definition
$$ u[n] = \begin{cases}{ 1 ~~~, ~~~n \geq 0 \\ 0 ~~~,~~~ n < 0 } \end{cases} $$
Then you will have for $u[2k-3]$ as
$$ u[2k-3] = \begin{cases}{ 1 ~~~, ~~~ 2k-3 \geq 0 \\ 0 ~~~,~~~ 2k-3 < 0 } \end{cases} $$
$$ u[2k-3] = \begin{cases}{ 1 ~~~, ~~~ k \geq 2 \\ 0 ~~~,~~~ k< 2 } \end{cases} $$
I think this last line is easy to see. You can plot it wrt the variable $k$.