What's the difference between using DFT, IDFT or DCT to calculate cepstrum of a power spectrum?

They are all somewhat different, in a family of potential candidates,  depending on the system you use and the features you are interested in.

Cepstra play on anagrams related to standard Fourier terms:  quefrency alanysis, liftering, cepstrum, saphe (for the phase)  etc. Not surprisingly, one can find some mixing in cepstrum definitions. You can think of it as a spectrum of a spectrum, scrambled.

The original cepstrum was a squared magnitude of the inverse Fourier transform of the logarithm of the squared magnitude of the Fourier transform of a signal (Bogert et al. 1963). But people have developed generalized versions (related to homomorphic systems). This is where the apparent inconsistency arises.

From a very mundane point of view, take a transformation $F$ (often orthogonal) that gives you some kind of spectrum. Separate (or not) a phase component ($\angle\left( F\left( x\right)\right)$, which could be the sign, some official phase or a mere constant) and an amplitude component $|F\left( x\right)|$. Put an exponent $p$ to the latter, add a little constant $k$ to it to avoid it from vanishing. Pour everything in a companding function $L$: the  $\log$, or even a Box-Cox transform: $(x^\lambda-1)/\lambda$. Now, transform the result  with $F^*$, which can either be $F$ or its inverse $F^{-1}$:

$$F^*\left(\angle\left( F\left( x\right)\right).L\left(k+ |F\left( x\right)|^p\right)\right)\,.$$

Now you can look at the absolute value of this, or a power again, or some inverse of the function $L$.

The key features are:

what is a meaningful frequency information for you or your system? If your system natively uses a DCT, then it is natural to compute a spectrum in a DCT domain.
computations: indeed, a few orthogonal transforms have been used for cepstra: sine, cosine, or Hadamard, which can be quite fast.
features looked at: from basic to invertible.


Why use a 1-bit ADC in a Sigma Delta Modulator?

First of all, because it's easy to build a 1-bit ADC. It's a comparator. It's literally the easiest ADC you can build. The $\Delta\Sigma$ ADC was invented (or, rather, published) in 1962¹ !
The 2-bit ADC is more than twice as complex as that,  you need some window decision: so if you have the choice of making your 1-bit ADC run faster or building a somewhat exact 2-bit ADC, there's a solid chance you'd go for the 1-bit ADC, simply because there's less analog semiconductor design to do, and, more importantly, to go wrong!
Imagine the 1-bit ADC doing a sign decision ("is the analog voltage > 0V?"). No matter how you scale the analog voltage (multiply it with a factor $\alpha$), the result will always be the same.
Thus, the 1-bit ADC-based $\Delta\Sigma$-ADC can be scaled to fit any range simply by changing (scaling) the threshold (which effectively is the analog part of your 1-bit ADC) in the forward chain (picture from wikipedia):

It's absolutely not obvious how you'd change above diagram to fit a 2-bit ADC! But, above diagram is very close to the original publication¹, so abstraction yields that the threshold is just a single-bit quantizer.
However, it's been done, there's multiple patents² covering higher-quantization  methods. You'd have to limit the bandwidth of the forward chain³, and, depending on what you want to achieve (class-D amplifier, faster ADC, lower power consumption…), you'd then have a bit of clever logic that translates the ADC output into different kinds of feedback "pulses" and different kinds of counter increments.
So, your original claim "we don't use more than one bit in the quantization within the loop of a $\Delta\Sigma$ ADC is wrong; it is only right for the very classical implementation of the original inventors.
To totally disprove your point: There's a lot of $\Delta\Sigma$ ADCs on the market that actually do higher-order quantization; for example, the Analog Devices AD9267 uses a nine-level quantizer⁴ (3.something bit, which was a pretty impressive feat on its own to integrate, considering the speed of 640 MS/s of that quantizer).

¹  Inose, H., Yasuda, Y. and Murakami, J., 1962. A telemetering system by Code Modulation – $\Delta$-$\Sigma$ Modulation. IRE Transactions on Space Electronics and Telemetry, (3), pp.204-209.
² WO2008028142, which already cites multi-bit feedback as "prior art".
³ otherwise, you'd always only get the lowest ADC output before you could get any of the higher ADC outputs, and with infinite bandwidth, that would mean you'd already emit a pulse to "reset" the integrator, and you'd never use any but 1-bit of the ADC. That bandwidth-limiting typically happens "involuntarily" because electronics are limited by physics.
⁴ AD9267 Datasheet, p. 13, "Theory of operation"


Beginner Attempting FFT Signal Filtering With Numpy

Two potential problems with your approach:

You are computing the FFT on your whole signal, which will be terribly inefficient if your input data gets too large, and prevents you from using windowing. If you want to do frequency-domain modifications of your signal, consider using a Short-Term Fourier Transform, processing the resulting FFT frames, and resynthesizing a time-domain signal using overlap-add.
FFT -> zeroing FFT coefficients -> IFFT, especially without windowing of the input data, is seldom used for filtering as it will yield a filter with many unwanted characteristics (side-lobes + non-causal).

A potential problem in your code:

I don't understand why you zero the filtereddata[len(data) / 2 + f] coefficients. I assume it's intended to zero the "negative frequencies" but then the array index is wrong, and you are using numpy.fft.rfft which has a different output layout from numpy.fft.fft (FFT coefficients for negative frequencies are not returned). I'm surprised this runs at all, your code should actually cause an index out of bounds error.

My suggestion would be for you to learn about FIR / IIR filters, design a suitable filter matching your constraints, and apply it to your data. There is a well-established body of work, software tools and practices on FIR / IIR filters and their design. Once you have the coefficients, for small sequences and/or if complexity is not an issue, there's scipy.signal.lfilter to do the job. As suggested by Jason R, applying a long FIR can be made very efficient computationally with FFTs + overlap-add.


Periodicity of discrete time signal: $ x\left [ n \right ] = \cos (\frac{\pi n^{2}}{8}) $

You want to prove that $\frac{nN}{8}$ is an integer for any integer $n$. Consider $n=1$. Clearly, $N$ cannot be less than 8.

You also need to prove that $\frac{N^2}{16}$ is an integer. This means that $N$ is a multiple of 4. So, the smallest $N$ that meets both conditions is 8.


It is easier to do this in HSV color space than RGB color space.
Assume the first image(the pinky one) is imported into matlab workplace as variable A,and the second(the yellowish one) as variable B.
First, I removed the white margin of both A and B with a same Rect:

Rect=[82.51 32.51 511.98 509.98]; % I found this with [~,Rect]=imcrop(A);
A_c=imcrop(A,Rect);
B_c=imcrop(B,Rect);

I am going to find a simple way to transform image A_c to B_c in HSV color space;

A_c_hsv=rgb2hsv(A_c);
B_c_hsv=rgb2hsv(B_c);

Read the wikipedia page of HSV.. you will find that most different parts  between the two images must be in the Hue parts,i.e,  between A_c_hsv(:,:,1)  &  B_c_hsv(:,:,1) ;

And the average difference between their hues is:

ave_dif_in_hue=mean(mean(B_c_hsv(:,:,1)-A_c_hsv(:,:,1)));

I just used the hue parts in my transform:

A_c_hsv_transformed=A_c_hsv; % make a copy of A_c_hsv
A_c_hsv_transformed(:,:,1)=A_c_hsv_transformed(:,:,1)+ave_dif_in_hue; % change the hue parts only

let's check our transform:

figure;
subplot(1,3,1);imshow(A_c);title('original A_c');
subplot(1,3,2);imshow(B_c);title('original B_c');
subplot(1,3,3);imshow(hsv2rgb(A_c_hsv_transformed));title('transformed from A_c');

here is the result:

you can optimize the transform further, with dealing with the other two channel in HSV color space..
Or.Maybe something more needed.


How does filter works in noise reduction?

The illustration is showing your data (labelled as "original") and the result of convolving your moving average filter with your data (labelled as "smoothed"). The filter/kernel is [1 1 1 1 1]/5 and gives each of the five samples equal weight. Effectively, it takes the average of the five samples it overlaps with during each step of the convolution process. It can also be written as [0.2 0.2 0.2 0.2 0.2] or [1/5 1/5 1/5 1/5 1/5]. The values in this filter are the same, but you can also choose different values depending on what you wish to accomplish. In the next slide, they give the filter [1 4 6 4 1]/16. This filter gives more weight or importance to the samples that are near the center of the filter.
I hope this helps.


$$ H(s) = {K*\frac{s-z}{s-p}}\\
\\
$$

when |z| < |p|, you have a lead compensator. It adds phase between a certain band, and can help you improve your phase margin. It can also increase your bandwidth. You typically use it to improve your transient as you mentionned. You can think of it as a PD controller cascaded with a low-pass filter. The gain K is optional, but if you set it to $\frac{p}{z}$, the DC gain is 1 which is usefull if you cascade the lead compensator to an existing controller.

I actually never used a PID with a pure derivative in real-life. I have always used a PI controller cascaded with a lead compensator, which is almost the same as a PID but as the advantage of limiting the high-frequency gain, thus improving stability at high frequencies.

When |z| > |p|, you have a lag compensator, which is almost the same as a PI controller. However, you won't have an infinite DC gain. Think of it as a PI controller with a soft reset. I have also seen the term "leaky integrator" used. It could be useful as an alternative to a PI to prevent integrator windup, or if you want a PI controller with a soft reset to slowly reset the integrator.

I have never really used one in real-life. Because, I use appropriate anti-windup strategies. And in the few cases where I needed a "soft reset" for my integrator (for a PLL), I implemented a conditional soft-reset that was only enabled in certain conditions.

If I have time, I will add the Bode plots of both.


Yes.

However, if you are dealing with a discrete-time signal, if your A/D is proceeded by an ideal low-pass filter with cutoff frequency $F_s/2$, with unity passband gain, the output power spectral density before sampling is:
$$
S_y(f) = |H(f)|^2 S_x(f) = \frac{N_0}{2} \operatorname{rect}(f/F_s)
$$
This means that the autocorrelation of the filtered noise will be:
$$
R_y(\tau) = \mathscr F^{-1}\{S_y(f)\} = \frac{N_0}{2T_s}\operatorname{sinc}(\tau/T_s)
$$
where $T_s = 1/F_s$. If the noise is zero-mean, Gaussian, and sampled at time instants $t_n = n T_s$, then each sample will be uncorrelated, and since they are Gaussian, will also be independent.   In this case, in the digital domain, the power spectrum will be flat for all frequencies.

This is not the only way to generate noise that is white in the discrete-time domain, but is one way, albeit a physically unrealizable one.


Confusion understanding phase shift/delay?

Delay is the negative derivative of phase with respect to frequency. Phase that increases linearly negative as frequency increases will result in a positive (later) delay in time. This is given by the Fourier Transform property of time delay:
$$\mathscr{F} \{f(t-\tau)\} = e^{-j\tau \omega}F(\omega)$$
And is easy to prove from the Fourier Transform:
$$\mathscr{F} \{f(t-\tau)\} = \int_{-\infty}^\infty f(t-\tau)e^{-j\omega t}dt$$
$$ = \int_{-\infty}^\infty f(t-\tau)e^{-j\tau \omega}e^{j\tau \omega}e^{-j\omega t}dt$$
$$ = e^{-j\tau \omega}\int_{-\infty}^\infty f(t-\tau)e^{-j\omega (t-\tau)}dt$$
If we substitute $u=t-\tau$, then $du = dt$ and we get:
$$ = e^{-j\tau \omega}\int_{-\infty}^\infty f(u)e^{-j\omega u}du =e^{-j\tau \omega}F(\omega) $$
Also without going into the details of the Fourier Transform, consider how we subtract phase in the time domain to shift a sine wave in time to the right:


Linkwitz-Riley Crossover Sum as Allpass Filter

Just adding to the answer. The sign for recombining L/R filter alternates, so it's '-' for second order '+' for 4th order and so forth.
The algebra is isn't all that pleasant so I will only go through the 2nd order. A second order L/R lowpass is simply the cascade of two first order Butterworth lowpass. The first order butterworths are
$$L = \frac{1}{1+jx}, H = \frac{jx}{1+jx}$$
where x is the normalized cutoff frequency $x = \omega/\omega_0$
That makes the difference of the L/R filters
$$H = L^2-H^2 = \frac{1^2}{(1+jx)^2} - \frac{(jx)^2}{(1+jx)^2} = \frac{1+x^2}{(1+jx)^2}$$
The numerator is real, so we just have to proof that the magnitude of the denominator is the same.
$$|(1+jx)^2| = |(1+jx)|^2 = (1+jx) \cdot (1-jx) = 1 + jx =ix + x^2 = 1+x^2$$
which indeed equals the numerator.
Another interesting class of "recombining to allpass" are odd order Butterworth filters. Here the lowpass and the highpass are 90 degrees out of phase so both the sum AND the difference of the filters will be create an allpass, although one has significantly more group delay than the other and the "optimum" choice of the sign also alternates with order.
Odd order Butterworth also have only -3dB of gain at the crossover whereas L/R have -6dB.


How to find H(z) and H(k) from a given causal function

It is a Homework problem, so only Hint will be provided. You need to show some work of your own following the hints.

First of all, take $\mathcal Z$-transform of the equation you have. Take care of the delays in time-domain and apply appropriate property of $\mathcal Z$-transform.

Then, try to get the transfer function $H(z)$ from there.

Once you have $H(z)$, figure out a way to convert $H(z)$ into $H(e^{j\omega})$ and then try to proceed for $H[k]$.

Show your work by making updates in your question and we will be clarify the doubts and be able to help you understand.


Computing Discrete Convolution in terms of unit step function

As a first remark, I would really suggest to also familiarize yourself with the graphical way of convolving two functions or sequences. This will also make your calculations a lot easier because you get some intuition on the nature of the solution.
Second, your definitions of $x[n]$ and $h[n]$ in terms of the step sequence $u[n]$ are not correct. E.g., the given $x[n]$ equals $1$ for $n=8$, but $u[n-3]-u[n-8]$ is zero for $n=8$.
It is formally correct to represent the solution as the sum of four individual convolutions of shifted step sequences, so what you could do is derive a general solution of the individual convolutions of the form
$$u[n-n_0]\star u[n-n_1]=(n-N_0)u[n-N_1]\tag{1}$$
I leave it up to you to derive the correct values of $N_0$ and $N_1$ in terms of $n_0$ and $n_1$. Once you have obtained that result, you can just compute the solution to the original problem as the sum of four terms of the form $(1)$.


Frequency magnitude distribution of noise

Let me try to elaborate a bit on your questions: First, the Wiener-Khintschin-Theorem relates the autocorrelation of a stationary process to its PSD, with the relation of the Fourier Transform, i.e.
$$\mathcal{F}\{r_{nn}(\tau)\}=S_{nn}(f).$$
The autocorrelation of a signal describes the correlation between two noise samples at distance $\tau$. It is just the second moment of the (joint) distribution of both noise samples. Hence, it does not really tell you something about the actual distribution. However, mostly, the noise is assumed to be Gaussian distributed with zero mean. Then, the autocorrelation completely describes the noise process.

Furthermore, in case you have AWGN, then $r_{nn}(t)=N_0\delta(t)$, i.e. the noise samples are uncorrelated and the PSD is a constant with $S_{nn}(f)=N_0$.

Your question is: What is the distribution of the noise magnitude in the frequency domain? So, you want to know $Pr(S(f)|_{f_0})$, where $S(f)$ is a realization of the noise process.

To make this more accessible (at least to me), let me switch to a linear algebra setting (i.e. the finite discrete-time case). The results can be straight-forwardly taken over to the continuos-time case. Let $\vec{n}$ be a realization of the noise process with autocorrelation $R_{nn}=E[\vec{n}\vec{n}^H]$. Let $\mathbf{F}$ be the unitary DFT matrix. So, $\vec{N}=\mathbf{F}\vec{n}$ is one noise realization in the frequency domain. Assuming that $\vec{n}$ consists of Gaussian distributed elements with zero mean, $\vec{N}$ also consists of Gaussian elements (A linear combination of two gaussian random variables remains to be Gaussian). What is the autocorrelation of the frequency-domain noise? It's given by

$$\mathbf{R}_{NN}=\mathbf{F}R_{nn}\mathbf{F}^H.$$

On the diagonal of $R_{NN}$ there is the variance of the Gaussians in the frequency domain.

So, the distribution of noise in the frequency domain remains to be a Gaussian, with variance given by diagonal on the above expression. I believe the diagonal is actually the PSD that follows from the Wiener-Khintschin Theorem. So, the final answer is:

  The noise in the frequency domain is Gaussian distributed, with zero mean and variance given by the PSD of the noise. (Assuming zero-mean Gaussian time-domain noise)


A first short shot, because I could waste hours of your time on that.

The reality (whether it exists or not) looks, sounds, appears complicated
As our senses are limited, and our brain is too slow, we need to shrink our perceptions, enough to understand most of what we need to grasp and share
So with our limitations, we need (and are only capable of) approximations of the above reality, here some signal $s(t)$
Depending on the time we have, the purpose, the extent, different kinds of approximations are required.

Convergence is the science of 'how approximations can be expected to be close enough to the reality'.

To draw a circle by hand, you only need to seize that $\pi \approx ~ 3$ (uniform 5% error). With a compas on a gridded sheet graph paper, perhaps $\pi \approx ~ 3.15$  is sufficient. For the Jet Propulsion Lab (JPL) "for interplanetary navigation, we use 3.141592653589793", see details in  How Many Decimals of Pi Do We Really Need?

Convergence is about how more complex questions are answered (see for instance how can a function can be approximated by sines or cosines) :

can a calibration function be approximated by a parabola
if I want  raise of 20 % in two years, can I have an increase of only $9.6 \%$ each of the two years (th answer is yes).

In more academic terms: suppose that there is a function $s(t)$. One has little knowledge about it, except that it wiggles. So you can choose a set of functions $\psi_\omega (t)  $ depending on a parameter $\omega$, that wiggles like exponential sines. Can we approximate $s(t)$ with a certain combination of  a finite number of $\psi_\omega (t)$? In the standard linear case, we seek a linear combination of the $\psi_\omega (t)$:

$$\sum a_\omega \psi_\omega (t)$$

The question rephrases as, for a given interval for $t$, and an increasing set $\Omega_\lambda$, $\lambda\in \Lambda$, of $\omega \in \Omega_\lambda$. The causal ideas, is: the bigger the number of terms (so the bigger the cardinal or the number of terms in  $\Lambda$). In other terms, is

$$s(t)-\sum_{a_\omega \in \Omega_\lambda} a_\omega \psi_\omega (t)$$

very close to $0$, and how, depending on $t$ and $\Omega_\lambda$. Or: how fast does

$$\Delta\left(s(t)-\sum_{a_\omega \in \Omega_\lambda} a_\omega \psi_\omega (t)\right)$$

tends to zero, where $\Delta(\cdot)$ is a measure of closeness (a distance, a norm, a divergence).


Haar Wavelets vs STFT for graceful audio quality degradation (newbie question)

For some non-stationary 1D signals, containing sudden changes, or a lot of noise, wavelets can have a close to optimal performance. Here, by wavelets I mean short oscillating shapes spanning several frequency ranges or scales.
Wavelets are somehow better than Fourier to that respect, as they can deal with peicewise regular + jumping signals.

However, this covers a huge range of possibilities for continuous and discrete tools.

Hence, I am confused with your "Haar Wavelets vs STFT" title. Haar wavelets can be orthogonal, STFT very redundant. So one cannot compare them directly in terms of sparsity and quality of approximation.

I will thus hypothesize that you performed a Hamming window of some audio signal frame of length $L$. Then applied either an FFT or an orthogonal Haar discrete wavelet transform on the result. Then performed a nonlinear approximation by keeping the $N$th largest coefficients, with $L \ll N$.

The result you evoke is a bit surprising to me, and I am a wavelet supporter: it seems a little too good for Haar wavelets. Haar wavelets are often  not the best fit for audio, and few audio coders use wavelets. What makes me wondering about your results is:

the length of your frame, or is the signal sufficiently nonstationary to win by KO on Fourier?
"starting from low to high energies": I have problem understanding what it means with wavelets. Are you performing a fair treatment on Fourier ?
magnitude and complex symmetry: Fourier coefficients are usually "doubled" on real signals. Did you take that into account? Are you sure you have reconstructed Fourier coefficient fully, and not from magnitude only?

All in all, a sketch of your code could help us deciphering the results.


How i can create ghost on wavelet or seismic trace and then Deghosting that in MATLAB?

YOu can use SeisLab, which apparently has an option for that:

  0,6.1,'s_ghost_response','Apply effect of source/receiver ghost on seismic dataset';


How do I scale Gaussians for Edge detection

Does this means I have to manually adjust the filter sizes for any new image?

yes! You need to know the kind of edges you want to detect. And that "kind of edge" is a "width of the edge in pixels" in disguise. So, you need more "intelligence" than just a DOG if you want to detect your favourite kind of edge in pictures of any zoom / resolution.

  Is there an automatic way to do so?

Plenty! You could, for example, simply try out some parametrizations for your DOG and automatedly select the parametrization that leads to an edge image that fits your signal model best – for example, if you want to find square boxes in images, looking for the "cleanest" straight lines.

Again, all depends on what you want to detect, and how you model your signals. There's no "universally" best edge detector. (or else, that would be the only one everyone uses, right?)


How to plot ROC curve for weka multilayerperceptron?

From University of Waikato WEKA wiki:

  The class used for calculating the ROC and also the AUC (= area under the curve) is weka.classifiers.evaluation.ThresholdCurve.


what are and why are sine and cosine modulated integrals used?

As someone made me notice in the comments, the integrals are just the coefficients of the Fourier Series for that frequency component (when the series is expressed with the cosine and sine parts separately).

The integral is performed on a period much bigger then the period corresponding to the frequency of the disturbance so that a better averaged is obtained.

To get back the signal in the time domain it is sufficient to multiply the two components by the cosine and sine respectively. Of course in this way the signal has only one frequency component.


Capacity loss of fast fading regarding AWGN

Standard Rayleigh fading is $C\mathcal N\left(0,1\right)$ then $E\{\lvert h \rvert^2\} = 1$. Using Jensen's inequality $$\begin{align}\Delta&=E\left\{\log\left(\lvert h \rvert^2\right)\right\}\\ &\ge \log(E\{\lvert h \rvert^2\}) \\&= \log(1) \\&= 0 \\&> -0.83\end{align}$$ and this $\Delta$ should not be $-0.83$? I think I made some mistakes but I don't know what it is.

Rayleigh fading is when the amplitudes of the real and imaginary part of $h$ are normal distributed; so $|h|$ follows a Rayleigh distribution, which has a probability density function  of

$$
p(x=|h|) = \frac {x}{\sigma^2}e^{\frac{x^2}{2\sigma^2}}\text.
$$

Hence, $E\left\{\lvert h\rvert^2\right\}$ is not the variance of the Normal distribution with variance $\sigma^2$.

Instead, the value you're looking for is
$$\begin{align}
E\left\{\log_2\left(\lvert h\rvert^2\right)\right\} &= \int_{-\infty}^\infty \log_2\left(\lvert x\rvert^2\right)\frac {x}{\sigma^2}e^{\frac{x^2}{2\sigma^2}}\,dx\\
&= \int_0^\infty 2\log_2\left(\lvert x\rvert\right)\frac {x}{\sigma^2}e^{\frac{x^2}{2\sigma^2}}\,dx\\
&= \frac2{\sigma^2}\int_0^\infty x\log_2\left( x\right) e^{\frac{x^2}{2\sigma^2}}\,dx
\end{align}$$


Convolution in Spatial Domain Is Multiplication in Frequency Domain

For finite size arrays, multiplication in the frequency domain is equivalent to circular convolution in the time domain, not linear convolution.  (Is that toroidal convolution in the 2D case?)

Note that with sufficient zero-padding, the results of circular convolution and linear convolution end up identical.


2D Image Convolution: Spatial Domain vs. Frequency Domain Convolution in the Computational Complexity Sense

Vanilla implementation of each method for image of size m x n and kernel of size k x l will yield:

Spatial Domain Convolution - O(mnkl) as for each pixel in the image we do kl multiplications (Additions are discarded).
Frequency Domain Convolution - O(mn log(mn) + mn) as the complexity of the FFT is mn log(mn) and we add the multiplication (You could add factor of 4 for doing it twice and back and forth).

But in practice there many nuances:

Many times the kernel is Separable or can be approximated by 1-2 separable filters (For instance, Gaussian Filter).
If the kernel is very small compared to the image, better use direct method.
If you apply the same filter on many images you may even consider using Matrix Form (As done in CNN's).
In Frequency Domain you apply Convolution with Circular Boundary Condition. It means that if you're after different boundary conditions you'll need to pad and then complexity is higher and many memory operations are done.

I can go on and on...

You may have a look at Filtering n×n Images by Separable m×m Filters. Computation Time for Filtering Using FFT, 2D Convolution and Two 1D Convolutions.

In order to say which is better, one must be very specific about the case.
In most cases, use efficient spatial domain code. This is my choice in most cases.


Assuming from your screenshot that the frequency magnitudes of "key clicks" are always higher (brighter) than the normal signal, you could just employ a basic Limiter? In other words, if a frequency magnitude value exceeds a certain threshold, set that frequency value to its last (valid) value. If you wanted to keep the plot from looking unnatural, you could simulate the data where the clicks are or use data from previous samples.


Why Z-transform is considered as separate transform?

There is indeed a transform called discrete Laplace transform and it is of course closely related to the $\mathcal{Z}$-transform. The (unilateral) discrete Laplace transform of a sequence $f_n$ is defined by (cf. link)

$$\mathcal{L}_T\{f_n\}=\sum_{n=0}^{\infty}f_ne^{-snT}$$

with some $T>0$. The discrete Laplace transform can be interpreted as the Laplace transform of a sampled function $f(t)\cdot\sum_n\delta(t-nT)$  with $f_n=f(nT)$.

In practice it is not convenient to have the factor $e^{sT}$ appear everywhere. If one substitutes $z=e^{sT}$, the discrete Laplace transform is called (unilateral) $\mathcal{Z}$-transform:

$$\mathcal{Z}\{f_n\}=\sum_{n=0}^{\infty}f_nz^{-n}$$

The same can obviously be done for the bilateral versions of the transforms, where the integrals and the sums start at $-\infty$.


Applying frequency-domain filters on a centered Fourier transform

Note that an FFT of strictly real data is conjugate symmetric.

For a length M FFT of strictly real data, the data in FFT result bin M-1 has the same magnitude as in bin 1 (for a low frequency), but complex conjugated.  So to maintain symmetry after filtering, which is necessary for the IFFT-ed image to remain strictly real, a filter has to modify bin(M-i) by the same ratio as bin(i).


Constant modulus algorithm - performance poor?

The constant modulus algorithm does not work for QAM because QAM is not a constant modulus signal (amplitude at the symbol decisions is not constant), such as BPSK and QPSK.

SIR is the "Signal to Interference Ratio" and the greater the SIR the better the performance, not the other way around as you describe. With a higher SIR, the signal is stronger relative to the interference.


period detection for a series of photons

As discussed in the comments to your question the recorded data that consists of arrival times of photons must first be converted to photon counts $n$ for this is the light intensity. The result of this conversion will be a vector $n(k)$ (with discrete time $k$) that represents the evolution of light intensity over time and that can be analyzed in more detail with an FFT, for example. The Lomb-Scargle method won't help here as you'd like to analyse light intensity, not photon arrivals. It could be of use if your data was a sequence of unequally space light intensity values. But as I understand your question this is not the case.

In order to obtain $n(k)$ count the number of photon arrivals in every time interval $t_0 + kT$, where $t_0$ is some constant time offset (e.g. first arrival time) and $T$ is the length of the equidistant time intervals. $T$ determines the sampling frequency $f_\mathrm{s} = 1/T$ of $n(k)$ and consequently the range of any spectral analysis, that is $0 \ldots f_\mathrm{s}/2$. On the other hand $T$ should be sufficiently large to obtain a reliable photon count in each time interval.

To analyse the periodicity of $n(k)$ calculate its FFT and look for distinct peaks. Functions like periodogram also come handy sometimes.


The DFT of a discrete-time sinusoid will only have components at the sinusoid's frequency if there are an integer number of periods inside the window. I.e., you must figure out the period $M$ of the sinusoid and choose the number of points $N$ as an integer multiple of $M$ ($N=kM$). I'm sure you can figure out the period of the given sinusoid.

Also have a look at this answer to a related question.


In edge detection, what value is being compared to the threshold value?

Its simple. The gradient matrix is computed while performing the Canny/Sobel edge detection. No the gradient varies according to the edge. When you digitize this gradient function, you get a value which tells the magnitude of the gradient at that point.

When you threshold. These magnitudes are compared to the threshold value. If you want weak edges, you keep low threshold values. In such a case all minor edges will show up as well.
If you want only strong 'clear' edges, then you keep higher threshold values.


Find transfer function from root locus and step response diagram?

The transfer function is

$$H(s)=\frac{K(s - 4)(s + 5)}{(s + 3)(s + 6)(s + 10)}=\frac{K(s^2 + s - 20)}{s^3 + 19s^2 + 108s + 180}$$
So we need to find $K$.

The LT of input step is $\frac{3}{s}$. Using the FVT for the step response $g(t)$:

$$\displaystyle\lim_{t\to \infty}g(t)=\displaystyle\lim_{s\to 0}s\cdot \frac{3}{s}\cdot\frac{K(s^2 + s - 20)}{s^3 + 19s^2 + 108s + 180}=-4$$
or
$$\frac{3\cdot K\cdot (-20)}{180}=-4$$
which gives $K=12$.
Then use the following code to plot the step response

K= 12;
sys = tf(K*[1 1 -20],[1 19 108 180]);
step(3*sys)

So there is no feedback.


Discrete-time Fourier transform of $a^{|n|} u[n]$

There's an obvious typo in that DTFT pair, because - as you've probably noted yourself - we clearly have $\alpha^{|n|}u[n]=\alpha^nu[n]$.
The given expression for the DTFT corresponds to the two-sided time-domain sequence $\alpha^{|n|}$, if $|\alpha|<1$ holds.


Constant Magnitude Chirp across audio band (Parseval's Theorem?)

If you change the "X = chirp(t,fo,10,f1,'logarithmic');" to "X = chirp(t,fo,10,f1,'linear');" you will get the response that you want.

To have the energy be equal it must spend an equal amount of time at each frequency.  With the linear sweep it does that.  With the logarithmic sweep it spends more time at the lower frequencies and much less time at the higher frequencies, resulting in unequal energy.


Separation of signals without loss of information about the signals themselves

This is probably going to take multiple steps.  As a start, here is what I'd do:

You want a filter that retains pulse information, but filters noise.  The most extreme you can do with this is a rectangular filter that's $100\mu\mathrm{s}$ long, or possibly a raised-cosine filter that's about the same length.
Assuming that the pulse timing is very consistent, phase lock to each pulse train.  This means you have one 1kHz PLL, and one 1.4kHz PLL.  When you get ready for this step you may want to ask a separate question, but if your pulses always start with the same timing (i.e., if a "long" pulse just ends $100\mu\mathrm{s}$ later than "normal") then you can lock onto the leading edge of each pulse.
Any time the PLL phases indicate a pulse will overlap, ignore both -- basically, output "null" on both channels.  You probably need to allow room for those $200\mu\mathrm{s}$ pulses you speak of.
Detect long pulses by looking for unexpected energy $100\mu\mathrm{s}$ later than the expected center of the filter output.
Detect pulse height by looking at the height of the filtered pulse where the PLL tells you.

The very beginning part of development for this would be to record your data into a file so you can play with it, run it through your best guess at a matched filter, and look at it.


Derive velocity and acceleration from position

There are several ways to calculate the derivative of a stochastic signal. I will show you two methods.

First of all, adjust the output of a second order transfer function to your data: $${\omega_n^2}/({s^2+2\zeta\omega_n s+\omega_n^2})$$ Once you are satisfied with the fit, put an additional $s$ in the numerator and you will get the derivative of your signal: $${\omega_n^2\times s}/({s^2+2\zeta\omega_n s+\omega_n^2})$$ You can use a bilinear transformation and convert this solution into a discrete one, to implement in a software easily, if you like.

Another solution is to use the spectrum and calculate the fft of your signal. Then filter it vanishing the small values in the higher frequencies. Now, take the result and multiply it by $j\omega$. This is taken from the Fourier Transform Properties Table to calculate the derivative of a signal. Finally, redo the signal with ifft.

I added here a code to exemplify these aforementioned algorithms. The output picture is added as well at the end of the post. I hope it helps.

% \brief: this code exemplifies two techniques suitable for
%         calculating the first-order derivative of an
%         stochastic signal.
%
% \algorithm 1: LTI transfer function integrated with Tustin.
% \algorithm 2: FFT derivative property.
%
% \author:     Luciano Augusto Kruk
% \web:        www.kruk.eng.br

% some constants:
fs = 1000; % [Hz] sample rate
T  = 1/fs;
t  = 0:T:0.5;
nT = length(t);
% signal + colored noise:
[B,A] = butter(10,250*T*2);
s     = sin(2*pi*20*t)+filter(B,A,0.1*randn(1,nT));

figure(1); clf;
subplot(3,1,1);
plot(t,s);
title('signal s(t)');
grid on;

% first technique: second order transfer function
qsi  = 0.7;
wn   = 2*pi*60;
wnT2 = (wn*T)^2;
b    = 2*T*wn*wn;
a1   = 4 + (4*qsi*wn*T) + wnT2;
a2   = (2*wnT2)-8;
a3   = wnT2 - (4*qsi*wn*T) + 4;
dsdt = zeros(size(s));
for i = 3:nT
    dsdt(i) = (1/a1) * ...
        ((-a2*dsdt(i-1)) - (a3*dsdt(i-2)) + (b * (s(i) - s(i-2))));
end
subplot(3,1,2);
plot(t,dsdt, t,2*pi*20*cos(2*pi*20*t))
title('algorithm 1');
ylabel('ds/dt');
grid on;
set(gca, 'ylim', 150*[-1 1])
legend('algorithm', 'real')

% second technique: fft/ifft properties
s_fft = fft(s);
idx   = 1:(nT/2);
f_fft = (idx-1) * (fs/2) * (1/length(idx));
idx   = round((nT/2) + ((-nT/2.5):(nT/2.5)));
s_fft(idx) = 0;
w     = 2*pi*(((-nT/2):((nT/2)-1))/nT)*fs;
dsdt  = ifft(s_fft .* sqrt(-1) .* ifftshift(w));

subplot(3,1,3);
plot(t,dsdt, t,2*pi*20*cos(2*pi*20*t))
title('algorithm 2');
ylabel('ds/dt');
grid on;
set(gca, 'ylim', 150*[-1 1])
legend('algorithm', 'real')


After going through the literature regarding HHT and EMD, I found that the "Huang" part of HHT comes from the fact that he is the one who proposed EMD in the first place. That explains the name of the method...

For more developments regarding EMD and HHT, I recommend the papers by Rilling et al. "On empirical mode decomposition and its algorithms". For the lucky ones who speak French, Rilling's PhD. thesis on EMD no longer seems to be available online; the document seems very complete and contains a very detailed mathematical analysis of EMD. Related articles are also available in English here on Google Scholar.

One may sum-up HHT this way:

EMD: decomposes the initial signal as a list of Intrinsic Mode Functions (IMF);
Hilbert tranform computes the instantaneous frequencies associated to the IMF (which are precisely well suited for such transformation)
Hilbert spectrum, meaning a representation of the amplitude of the IMF in a frequency/time domain using the instantaneous frequencies.

Here is a simple example based on the time signal: $y(t)=t^2sin(t^2)$ for $t\in[0;~16]$s. This signal features a quadratic growth of its amplitude combined with a linear growth of its frequency. Here is what a HHT of the signal leads to:

signal of interest

IMF obtained from EMD

$y(t)$">

instantaneous frequencies

Hilbert spectrum (white to black from 0 to maximum amplitude)

Some peaople see the HHT as a generalized Fourier Transform in the sense that the decomposition of the signal interest by HHT leads to both amplitude and frequency varying time signals.

A signifiant drawback of HHT lies in its sensitivity to edge effects (what happens for the signal close to its left and right boundaries in 1D). Several techniques exist in order to mitigate these effects. Articles by Rilling mentioned above go with mirroring strategies while other engineering oriented strategies may involve waveform matching techniques.


8*8 block matrix in JPEG image compression?

The lossy JPEG compression does not merely remove small coefficients in higher frequencies. It encodes them with a precision relative to a (relatively crude) visual perception model; most notably, horizontal and vertical frequencies are not quantized with the same precision. And as in many compression formats, it essentially assumes that the data is locally stationary.
If you apply the  DCT over the whole image, and  you quantify DCT coefficients, this quantization will effect the whole image. Imagine an image with a background checkerboard pattern, and a small zebra in foreground. With a whole-DCT compression, the zebra is likely to loose its stripes, because their energy is negligible with respect to that of the checkerboard. Even more, as JPEG applies DCT on chrominances as well, with down-sampling, color coefficient quantization is likely to produce false colors at places where they do not belong. With   block size larger than  about $8 \times 8$, several meaningful image gradients (parts of object edges) are more likely to happen simultaneously in the same patch.
Meanwhile, one can get slightly better results with $16\times 16$ blocks. That is for the perceptual part.
There is a computational part too. Small non overlapping blocks are easier to process, and require less memory, a still expensive part of electronic devices. Early JPEG 2000, applying wavelets on the whole image, have failed adoption despite better results, partly because of their memory footprint. Now, with the  advent of GPUs, processing blocks is becoming quite attractive. Last, a power of two ($8=2^3$) and clever tricks make the block DCT very efficient. In modern coders, standard often use different $2\times 4$, $4\times 4$... block sizes.
Last, but not least, it is not modular, as you need DCT implemetations for each image size.  Although it is possible with software,  chip makers are unlikely to like the hardware part.
So, applying DCT over the whole image:

is detrimental to  non-stationary images, and waster some local orientation,
can be costly in terms of memory, and hardware,
is not quite modular nor computationally efficient.

Thus being said, very stationary images would be better compressed with whole-size  DCT.


Your integrator has a negative sign:

$$H(s)=-\frac{1000}{s}\tag{1}$$

That's why its step response also has a negative sign:

$$y(t)=-1000\int_{-\infty}^tu(\tau)d\tau=-1000\cdot t\cdot u(t)\tag{2}$$

where $u(t)$ is the unit step.


How to find smoothed estimates of the derivative and second derivative of a signal?

It probably depends more on your data. Just know, since differentiation is a linear operation, if you choose any linear filter to smooth f' and f'', it is equivalent to smoothing f using that same filter, then taking its derivatives.

Can you post some pictures or more information about the signal you want to differentiate? Probably what you're looking for is some sort of lowpass filter to smooth the signal. A couple really simple options include a single-pole recursive filter like $y(n) = a \cdot x(n) + (1-a) \cdot y(n-1)$, or a Hann filter, which is just convolving the signal with a Hann window. The Hann filter option is nice because it's linear-phase. If you know the frequency range you care about, you can just design a suitable lowpass filter in the frequency domain.


I believe that if you will keep the phases and not convert it into magnitude (power spectrum) the IFFT will work. When you use the center frequency, the phase is 0, in other words, the real part is also the magnitude. As you move from the center frequency the error in magnitude growth as the phase growth. As result your error in recovery by means of IFFT growth.


Suggest an image segmentation technique to find the paste radii?

This problem can probably be solved easier if there is a known relationship between the camera position and the metallic device that is over the paste. If the device and the camera are both fixed, for example, it would be easy to use a picture of the device without the paste to eliminate it from the other pictures.

But here is one suggestion for measuring the diameter of the paste using only this given image. I'm using DIPimage, because I know it best. I'm sure the same can be accomplished with the Image Processing Toolbox.

% Read in image:
img = readim('https://i.stack.imgur.com/5FLCs.jpg');
img = img{2}; % Use the green channel
% Alternating sequential filters to keep only the dark big circular shapes
paste = opening(img,15,'elliptic','add max');
paste = closing(paste,30,'elliptic','add max');
paste = opening(paste,60,'elliptic','add max');
paste = closing(paste,90,'elliptic','add max');
paste = opening(paste,120,'elliptic','add max');

The sequence of opening and closing are meant to remove specular reflections and much of the metallic device. Here is what the image looks like after this step:

There are some bulges to the left still, I didn't put any effort into trying to remove those, because there is sufficient outline for correct measurement.

% Find paste outline:
labs = label(~threshold(paste,'otsu')); % threshold and label
D = traceobjects(labs,[],2,'convex hull'); % get object outlines
[~,index] = max(cellfun(@length,D)); % find the longest outline
D = D{index}; % keep only the longest outline
% Display:
dipshow(img)
hold on
plot(D(:,1),D(:,2),'r')

This shows the convex hull of the outline over the original image, to show that the filtering didn't move the edges of the paste:

Next, I'll use code to list antipodal pairs that I get from this blog post. These are points on opposite side of the convex hull. This is the reason I used traceobjects with the 'convex hull' option, the algorithm to list antipodal pairs only works with a convex hull.

% Signed area of triangle
area = @(pt1,pt2,pt3) ( (pt2(1)-pt1(1))*(pt3(2)-pt1(2)) - ...
                        (pt2(2)-pt1(2))*(pt3(1)-pt1(1)) ) / 2;
% Height of triangle, pt3 is the top vertex
height = @(pt1,pt2,pt3) ( (pt2(1)-pt1(1))*(pt3(2)-pt1(2)) - ...
                          (pt2(2)-pt1(2))*(pt3(1)-pt1(1)) ) / ...
                                               norm(pt2-pt1);
% Next point along the polygon
next = @(p,N) mod(p,N)+1;

diameters = [];
N = size(D,1);
p0 = N;
p = 1;
q = 2;
while area(D(p,:),D(next(p,N),:),D(next(q,N),:)) > ...
      area(D(p,:),D(next(p,N),:),D(q,:))
   q = next(q,N);
end
q0 = q;
while p~=p0
   p = next(p,N);
   listq = q;                       % p,q is an antipodal pair
   while area(D(p,:),D(next(p,N),:),D(next(q,N),:)) > ...
         area(D(p,:),D(next(p,N),:),D(q,:))
      q = next(q,N);
      listq(end+1) = q;             % p,q is an antipodal pair
   end
   if area(D(p,:),D(next(p,N),:),D(next(q,N),:)) == ...
      area(D(p,:),D(next(p,N),:),D(q,:))
      listq(end+1) = next(q,N);     % p,q+1 is an antipodal pair
   end
   disp(listq)
   for ii=1:length(listq)
      v = D(p,:)-D(q,:);
      d = norm(v);
      diameters(end+1) = d;
   end
end

Finally, we can use these distances between antipodal pairs to find the average diameter of the paste. We'll assume that about 1/3 of the boundary is affected by the bulges on the left that I didn't completely remove. Therefore we'll use the 2/3 smallest values, and take the mean of those:

diameters = sort(diameters);
diameters = diameters(1:round(length(diameters)*2/3));
disp(mean(diameters))

This gives me 773.25 pixels. You will need to know the size of a pixel to translate this value to physical units.


To your first question, yes, each symbol is sent on a separate subcarrier; that's the basic idea of OFDM.

I'm not sure where your confusion lies with your second question. It's well-known that the Fourier transform of a system's impulse response gives its frequency response. In this case, the system is the communication channel. If you choose the sample rate of your impulse response measurement properly, and use the appropriately-sized FFT, you can get a frequency response measurement whose bins line up with the subcarrier frequency grid. This can be used for equalizing the channel.


When sampling theory condition fails, what should I do?

Because you've sampled at 11kHz the Frequency spectrum becomes periodic with a period of 11 KHz. So you're 15kHz will not only appear at 15 kHz but also at every frequency that is a multiple of 11 kHz away from 15 kHz. So the frequencies at which the 15 KHz tone are given by $ f= 15\times10^3 +N \cdot 11 \times10^3 $, where $N$ is an integer.

By convention you are looking at the frequency bins from -5.5 kHz to 5.5 kHz. Therefore you're 15 kHz tone will appear at 4 KHz.

Note that if your original signal was real then there is also a freqency component at -15kHz. This also gets extended periodically every 11 kHz - so this will produce a corresponding peak at -4 kHz.


Calculate the power of a discrete-time signal

Apart from possible dimensional problems as discussed in the comments, your result is correct and the given solution is wrong. In general with such exercises, don't forget to check if the cross-term really cancels out (in this case it does).


peak-notch block is a two-way component. You can set the response type to peak or notch. Here is the properties page:

int the link: http://www.mathworks.com/help/dsp/ref/peaknotchfilter.html

Using this filter you can obtain desired output by summing input signal and output of the peak filter.


Log Distance-Normal Shadow Model and Friis

I verified all the steps and they are all correct.


The exact definition of dominant frequency?

As dominant frequency it is usually meant the one that carries more energy w.r.t. all the other frequencies in the considered spectrum.

Since in the first example  the two frequencies have the same absolute value I suppose they are both dominant frequencies. However the concept is not really important in these case, since you are just looking at the combination of two pure sinusoids of equal amplitude.

In the second case, the lowest frequency is both the fundamental and the DOMINANT frequency since it has the highest peak of the considered spectrum.

Reference: https://arxiv.org/pdf/1306.0103.pdf


What this gain actually means?

The factor between in- and output amplitudes for a signal in the passband

  Does it compensates for the amplitude scaling inherent to upsampling?

Exactly!

  If I would hypothetically have an "interpolation" filter which is not attenuating spectral replicas (let's call it unity filter) would this filter still have gain L, or in other words is gain L something what is by the definition within the interpolation filter and then the filter coefficients which define how well the replicas will be attenuated is something which depends on the design/quality of the filter?

You design a filter to have a certain gain; setting gain is easy, just multiply your filter with a scalar.


Interpreting the amplitude of signals in fourier transform

"I totally understand the concept of Fourier transform"

Lucky you if you really do. Some of us (me, in first place) don't (in totality). The Fourier transform (and its avatars) is a prototype for duality. Duality here means that you can represent a signal on some primal domain (time) onto a dual domain (here frequency). This transformation is meant to possess useful properties, to preserve signal information, and to add insight, like a somewhat easier interpretation.

Signals are traditionally represented as "amplitudes of something for each time". These amplitudes are not always absolute. But we can hope they are linearly related to the actual physical values. And that the linearity factor does not change over time. Under the above premises, sines and cosines, or complex exponentials (cisoids), are a very good way (perhaps the best) to model  the system. And in turn they can be used to represent the data, as a linear combinaison of shifted sines.

Since the sines are orthogonal, the squared magnitude of the coefficients weighting the sines is proportional to the energy of each corresponding sine frequency.

The proportional factor depends on how the coefficients are computed, often up to a scaling factor. The amplitude of a given frequency component can be directly computed with correlation. Note however that this coefficient should be taken with care, as outer conditions play a role: the finite horizon observation of a function/signal, amplitude discretization, noise, etc.


Difference between order and length of a filter

In your example order is $2$ and length is $3$.

Its has more to do with convention. Consider a polynomial in $x$:
$$ a_n x^n + a_{n-1} x^{n-1} + ... a_0 x^0 $$
The order of this polynomial is $n$, whereas, number of coefficients is $n+1$. In your case an FIR filter of order $L$, is a polynomial you write as:
$$ h_L \delta(n-L) + h_{L-1} \delta(n-L+1) + ... h_0 \delta(n)$$
The $z$-transform of which is:
$$ h_L z^{-L} + h_{L-1} z^{-L+1} + ...+ h_1 z^{-1} + h_0 z^{0}$$
Note this transform, which you might not have studied yet, just places the shift in $\delta(n)$ on power of $z$, e.g., $z^{0}$ and $\delta(n)$ correspond to no-delay, just scaling the input signal. $z^{-1}$ and $\delta(n-1)$ correspond to order $1$ and one delay.

I think real confusion comes from MATLAB, as it does not allow $0$ as first index value. You cannot write x(0) in MATLAB. The first index is $1$. Therefore, we have to start from x(1).


To ask "how fast can you sweep a sinewave?" is to ask "How fast can you hit $\frac{F_s}{2}$?".

Because, sweeping a sinewave is equivalent to modulating it.

A sine wave looks like:
$$x(n) = \sin(n 2 \pi f / F_s)$$

Let's focus on that instantaneous phase. I have put $n$ first, on purpose, to show that phase, in discrete time, is a line that departs from $0$ at a slope of $\frac{2 \pi f}{F_s}$. Even more descriptively, $2 \pi \frac{f}{F_s}$.

The sharpest slope for a given $F_s$ that can exist is at $\frac{2 \pi \frac{F_s}{2}}{F_s} = \frac{2 \pi F_s}{2 F_s} = \pi$.

Anything even slightly higher than that rate will start "catching" the sinusoid at its "past", or in other words, create aliasing.

So, how fast can you sweep a sinusoid?

Consider a cartesian plane. It's $x$ axis corresponds to time and it's $y$ axis corresponds to rate of phase change. We mark a red line on the $y$ axis at $\pi$. That's our limit. We can't cross that. The time axis changes in terms of $T_s = \frac{1}{F_s}$. So, theoretically, we could go $0 \ldots \frac{F_s}{2}$ in just one step, in just one $T_s$.

Is that going to look like a sweep?

It is going to look like a sweep, as much as a sinusoid looks like a sinusoid when it is sampled at some $\frac{F_s}{2}$. That's two samples per period and it looks exactly like a triangle and nothing like a sinusoid, but when you are sampling at $F_s$, that's the best sinusoid you can "see" at $\frac{F_s}{2}$.

This cartesian plane representation though, helps convey another point too. That is, how many frequencies did the sweep end up going through anyway? In the extreme example given above, instantaneous phase was so rapidly modified that effectively just two frequencies showed up, 0 and $\frac{F_s}{2}$. Sometimes, when you are trying to obtain the frequency response of a system, you might "miss" a part of it because your sweep is so fast that it "jumps" over the mode you should have excited.

Therefore, a more meaningful way to determine your sweep rate is to accept an average power per Hz (or frequency band) that you want your chirp to deliver and then find the rate of phase change that "allows enough samples" to deliver that power.

Hope this helps.


Linear Predictive Coding for general signals

There may be a misunderstanding about LPC here. For one, it is mainly used for speech because it was designed to fit the model of the human vocal system. The figure below is taken from Purdue's VISE lab in one of their speech processing exercises. The figure shows speech modeled as a filter that is excited by either an impulse train or white noise. LPC works well for speech because it gives you a concise way of describing the filter part. You also need to estimate the driving frequency of the pulse train, detect if a time segment is voiced or unvoiced, and control for the power of each time segment.

In the frequency domain, the a speech signal will look like

and all the LPC is expected to do is to fit the spectral envelope.

  I see why LPC is often used for speech compression - because the error e(n) is roughly white and so we can roughly model it by its variance.

The significance of a white, particularly additive white Gaussian, error signal is that it means that you have used all of the information. If your error is highly structured, it means that you could do a better job at modeling the system.

  Do people use LPC for things other than speech, and if so, how do they go about modeling e(n) to achieve compression?

I'm not aware of it being used for anything apart from speech. As stated above, the vocal system is being modeled by LPC, not the error. The vocal system is modeled using the source-filter model, and LPC is a good way of describing the filter part of the model.

I'm not sure about your second question, but I hope the rest of this help. Speech encoding with LPC is a really interesting, particularly because it works so well. If you want to do some matlab experiments, you could try Lab 9b on Purdue's VISE lab.


Control design: under what conditions can closed-loop poles be placed arbitrarily?

Unnecessary additional information that might help

From a pragmatic point of view, there might be a problem when trying to control a given plant using only its transfer function $G_p(s)$. In general, you test the plant putting a series of inputs and seeing how it reacts, measuring the respective outputs. From that, you can model the plant as a transfer function in the Laplace domain. However, the intrinsic dynamics of the plant may have hidden unstable modes (for instance, poles that cancel with zeros, which means you are not able to notice their presence in the system) and thus the control system you design wouldn't work as expected.

But that's in real life. In Control Theory, we are often offered the transfer function of the plant with all the poles and all the zeros. So let's assume that $G_p(s)$ doesn't show this weird behaviour I described above. In that case, one must notice (at least) three advantages that state spaces (SS from now on) offer:

Working with SS instead of transfer functions can be way easier. If the order of the system is too high, operating with matrices can be the only way to deal with the system.
SS provide an easy way to deal with MIMO systems. Using transfer functions for each pair of input-output would be too tedious.
SS let us see those hidden unstable modes I referred to before (looking at the eigenvalues of the matrix $A$), while in the transfer function they cancel out.

The actual answer to your question

Now getting to your question, you've asked what happens when the system is of low order (in your example, the plant is of second order) and SISO.

You can definitely find your transfer function performing the variable isolation you've written for $G_c(s)$. The problem with that is that some things could go wrong. For example,

$G_c(s)$ could be impossible to implement.
$G_c(s)$ could be unstable, thus making the control system non-sense.
$G_c(s)$ could have zeros in the RHP that might cancel out with poles of $G_p(s)$.

If none of those things happen (maybe there are some other aspects that one should take into account, but I couldn't think of any additional one), then go ahead: you are free to find your control system, $G_c(s)$, without using full-state feedback at all.

The thing is that SS offer us a lot more information than transfer functions do. This is easy to notice, as when working with SS we are looking at not only the input and output variables, but the other state variables as well, which are totally left apart when working with transfer functions and block diagrams alone. When doing full-state feedback, we are using this additional information to control the system (i.e. you don't just use the output, but each one of the state variables). This means that there are other tools we can use to stabilize the plant and we were not using before. If the design using full-state feedback is good, then the closed-loop will work as a charm and none of the problems I've described above (unstability, unstable poles and zeros cancelling out, etc.) will appear.


How do I show if this system is time-varying?

In order to show that a system is time-varying, it is sufficient to find one counter-example, where the response to a shifted input signal does not equal a shifted version of the original output signal.

In your case you might want to consider the signal

$$x(t)=\begin{cases}&1,& 0<t<1\\
&-1,& -1<t<0\\
&0,&\text{otherwise}\end{cases}$$

Now try to figure out the corresponding output signal $y(t)$. Next, take a shifted version of $x(t)$, e.g. $x(t-1)$, as the input signal and see what output signal you get. Is it equal to $y(t-1)$? If not, you've just shown that the system is indeed time-varying.


actually, the diode and transistor are also time-invariant components.  the transistor is even reasonably linear in a range of operation between cutoff and saturation.  even, for very small displacements, a diode can be sorta linear, but a diode in normal use is clearly not linear.


Difference Between Dirac delta and unit impulse function

The Dirac delta $\delta(x)$ is a continuous argument (generalized) function where $x$ is a continuous variable and can be time, space, frequency etc. Unfortunately, a simple definition like

$$\delta(x) =
\begin{cases}
\infty ~~~&,~~~ x=0 \\
0 ~~~&,~~~ x \neq 0
\end{cases}
$$

does or can not properly describe its behaviour. You need to define its properties under the integral sign, and that requires advanced mathematical stuff for a formal verification. But  a simplified Riemann definition of $\delta(x)$ is given as:

$$\int_{-\infty}^{\infty}f(x)\delta(x-a) dx = f(a)$$ for a sufficiently smooth test function $f(x)$. All properties of $\delta(x)$ follows from this Riemann (seeming) integral operation.

On the other hand, the discrete counterpart, the unit-impulse or the Kronecker delta is a sequence of integer index $m$ with very well defined behaviour as:

$$\delta[m] = \begin{cases}
1 ~~~&,~~~ m=0 \\
0 ~~~&,~~~ m \neq 0
\end{cases}
$$
These two functions serve the same purpose for discrete and continuous systems, hence share the same name impulse. Yet they are two different, distinct entities; the Dirac impulse is a very problematic concept, whereas the discrete  impulse is a mere simplicity.


Why should the last point be excluded when performing a least-squares fit of a periodic discrete time signal?

First (wrong) answer (for integrity) The $y$-value of the last point
is the same as the first one. As you apparently know the frequency,
this point comes in excess of the "fundamental period". It sounds like
this additional point comes like an implicit double-weight to the
first  point of the period.

Second take: I have tried to fit the data, with or without the last point. It seems to fit well, both cases.

%https://dsp.stackexchange.com/questions/71164/why-should-the-last-point-be-excluded-when-performing-a-least-squares-fit-of-a-d
clear;close all
%% Settings
T = 1.5; w=2*pi / T;
t = (0:0.15:1.5)';
y = [2.200 1.595 1.031 0.722 0.786 1.200 1.805 2.369 2.678 2.614 2.200]';

%% Fitting all points
ft = fittype(@(a1,a2,a3,x) a1+a2*cos(w*x)+a3*sin(w*x),'coefficients',{'a1','a2','a3'},'independent', {'x'});
f = fit( t, y, ft );
% Plot fit
subplot(2,1,1)
plot( f, t, y )
axis([t(1) t(end) 0.5 3])
grid on
title('Whole point set')
%% Fitting all points but the last
ft = fittype(@(a1,a2,a3,x) a1+a2*cos(w*x)+a3*sin(w*x),'coefficients',{'a1','a2','a3'},'independent', {'x'});
f = fit( t(1:end-1), y(1:end-1), ft );
% Plot fit
subplot(2,1,2)
plot( f,  t(1:end-1), y(1:end-1) )
axis([t(1) t(end) 0.5 3])
grid on
title('Whole point set minus 1')


Is this the correct solution for finding the samples of phase delay at a given frequency of an IIR BiQuad filter?

You almost got it, but the final result is wrong because the phase of the ratio of two complex numbers is the difference of the phases of the numerator and denominator, respectively, not the ratio of the phases:

$$z=\frac{u}{v}=\frac{|u|e^{j\phi_u}}{|v|e^{j\phi_v}}=\frac{|u|}{|v|}e^{j(\phi_u-\phi_v)}\Longrightarrow \phi_z=\phi_u-\phi_v$$

Furthermore, for finding the phase delay, you need to divide the phase by the (normalized) frequency in radians (as shown in this answer). The result is the phase delay in samples.


Relationship between the IDFT of a sampled DTFT and its discrete-time domain signal

For simplicity, I prefer a 1D notation without losing generality.

The IDFT associated with the properly taken uniform (in frequency) samples of a given valid DTFT can be shown to be the following.

Case-1: Assume a finite length signal $x[n]$ of length $N$, and consider its DTFT $X(\omega)$. Now if you take $M \geq N$ uniform (frequency) samples of $X(\omega)$ according to $$ X[k] = X( \frac{2\pi}{M}k) ~~~, k=0,1,...,M-1 $$ then the M-point inverse DFT of $X[k]$ will exactly be $x[n]$ padded with $M-N$ zeros; i.e.,
$$ y[n] = IDFT_M\{ DFT_M\{ x[n] \} \} = x_M[n] $$ So the sequence $x[n]$ and $y[n]$ have the same values for the first $N$ samples (i.e., all samples of $x[n]$ are retained in $y[n]$), but $y[n]$ have further (effectively padded) zero samples.

If $ M < N$ , then time-aliasing in $x[n]$ will happen and the first $N-M$ samples of $y[n]$ will be in error. But the remaning $2M-N$ samples of $y[n]$ will be identical to the original $x[n]$ at the corresponding index positions.

Case-2: $x[n]$ is infinite length (and stable) sequence and $X(\omega)$ is its DTFT. Then according to $M < N$ subcondition of case-1, all samples of $y[n]$ will be aliased. No samples of $y[n]$ will be equal to $x[n]$, except other than (pure) coincidence.


Eroding noisy borders of polygons in an image

I have a mask image of shape 500x500x1 that contains polygons.

. . .

These polygons [...] are labeled

These two statements are incompatible. Assuming that the third number is the color depth, then 1-bit implies a binary mask, but each class being assigned to a different level implies a color depth >1.
With a 1-bit mask, binary erosion is applied. Otherwise grayscale erosion is applied, which is why you see those "out of class" values at the borders of the mask.
The simplest way to deal with this is to apply the erosion at each class within the image iteratively. So, collect the unique values from the mask in unique_classes. Then, for every level in unique classes, threshold the mask image (this returns a binary mask) apply the erosion and then bring it back to the level it is supposed to be.
Alternatively, duplicate the mask into two images, hold one back to be the one that describes the classes and threshold the second one so that it becomes a binary mask. In this setting, you apply erosion once, to the second mask, but use the first mask to decide which class a given pixel is a member of.
Hope this helps.


Why does my zoom audio, when there are nearby microphones, seem to usually take an open loop format?

This effect is usually undesirable and is prevented by (acoustic) echo cancellation algorithms. These algorithms estimate the distortion of the incoming signal on its way to the microphone (via the loudspeaker and the room acoustics). That distortion is modeled by a filter and the output of that filter is subtracted from the microphone signal, such that the contribution of the incoming signal to the outgoing signal is minimized.
Zoom has an option to turn off echo cancellation, so you could try playing around with that.


A straightforward way of estimating the PSD is to compute magnitude squared DFT (a.k.a. periodogram).

Suppose we have an $N$ point sequence $x(n)$ and we would like to estimate its power spectral density (PSD) $P_{xx}(f)$.

Method 1: Periodogram

We compute the magnitude squared DFT. No averaging.
$$
\hat{P}_{xx}(f) = \frac{1}{N}\left |  \sum_{i=0}^{N-1} x(n) e^{-j2\pi fn/N} \right |^2
$$

Method 2: Periodogram Averaging

We instead use $M$ different realizations of the $N$-point sequence $x(n)$ giving us $M$ sequences $x_i(n)$, $0\leq i\leq M-1$. We compute the periodogram estimate $\hat{P}_{x_ix_i}(f)$ of each $x_i$ and average them to get the final PSD estimate:
$$
\hat{P}(f) = \frac{1}{M}\sum_{i=0}^{M-1} \hat{P}_{x_ix_i}(f).
$$

It can be shown that the variance of this averaged periodogram estimate is lower:

$$
\mbox{Var}[\hat{P}(f)] = \frac{1}{M} \mbox{Var}[\hat{P}_{xx}(f)].
$$

Note that it is a little unfair to compare Methods 1 and 2 because Method 2 uses more data ($MN$ data samples, instead of just $N$). If all you have is the original $N$ point sequence $x(n)$, you end up making a tradeoff between PSD estimator's variance and frequency resolution.

The PSD of the underlying signal is what it is. A computational routine will produce an estimate of this PSD. Whether a PSD estimate is good or bad will depend on your definition of good or bad. One commonly used metric variance. In that sense, periodogram averaging estimate $\hat{P}$ is a better estimator of the PSD $P_{xx}$ than $\hat{P}_{xx}$.

cf. Bartlett Method, Welch Method

Reference: Digital Signal Processing, Proakis and Manolakis, Ch. 14 (4th ed.)


The integer $N$ is the number of samples of $x[n]$ in one period. There need not be any relation to absolute time whatsoever, simply because $x[n]$ are just discrete data that are not necessarily samples of a function of time.

However, if the data $x[n]$ are samples of some continuous-time signal, then the time difference between two adjacent samples is, by definition, the sampling interval $T_s$, which is the inverse of the sampling frequency $f_s=1/T_s$.

In that case, the time corresponding to one period of $x[n]$, i.e., the time difference between two samples $x[n]$ and $x[n+N]$, equals $N\cdot T_s=N/f_s$.


How to estimate the frequency characteristics of a nonlinear chirp?

How about just using $\dfrac{d\phi}{dt}$ to compute the instantaneous frequency?

GNU Octave code (I just made up a sample rate of 2 Msps):

pkg load signal;
x = csvread('nonlinear_freq_fit.csv');
y = hilbert(x);
z = y - mean(real(y));
dphi_dt = diff(unwrap(arg(z)));
Fs = 2e6;
f_inst = (Fs/2)/pi * dphi_dt;
k = [0:length(f_inst)-1];
t = k/Fs;
plot(k, f_inst);
axis([5000 5800 -2e5 2e5]);
title('Instantaneous Frequency');
xlabel('Sample index');
ylabel('Frequency (Hz)');
grid on;


What is 'energy compaction' in simple terms?

Yes, I believe that your understanding of energy compaction is correct.

  Does that mean, for any x, if A is DCT matrix, energy of y will be
  more concentrated when compared to the x.

No, it does not.  All that is needed to prove that such is not the case is to show that there is an $x$ that is not compacted by the DCT.  White noise, for example, would not be compacted by the DCT.  The DCT is useful because in many real life situations the "signal" (e.g. audio, images, videos, etc.) tend to be "pinkish", i.e. tend to have most of their energy in the lower frequencies and so there is a natural assymetry that can be used to our advantage.

Another way to look at this is from an information theory perspective.  If a signal is not "white" (i.e. it doesn't have the same power at all frequencies) that implies that there is some correlation between samples in the time domain, which means that the samples have "information" about the value of other samples.  This mutual information implies that there is redundancy, and thus should be able to reduce the amount of data without losing information.

  Does this happen to every x or x has to satisfy some properties to get this energy compaction?

For frequency-based transforms like the DCT it is clear from the above that what is needed is for the signal to be non-white.  The more non-white it is, the more compaction can be achieved.

There are other kinds of transforms though, which presumably could compact signals based on other features.  That is pretty much the whole point of compressed sensing.


Computing shifted signal without first reconstructing

Because the sampling frequency is above the Nyquist rate, the original signal can be written in terms of its samples $y(n)$:

$$x(t)=\sum_{m=-\infty}^{\infty}y(m)\frac{\sin[\pi(t-mT)/T]}{\pi(t-mT)/T}\tag{1}$$

By setting $t=nT+T/2$ you get from (1)

$$z(n)=\sum_{m=-\infty}^{\infty}y(m)\frac{\sin[\pi(n-m)+\pi/2]}{\pi(n-m)+\pi/2}$$

and since

$$\sin[\pi(n-m)+\pi/2]=(-1)^{n-m}$$ you finally get

$$z(n)=\sum_{m=-\infty}^{\infty}y(m)\frac{(-1)^{(n-m)}}{\pi(n-m)+\pi/2}$$

So for computing $z(n)$ you do not need to reconstruct $x(t)$ but you need to compute a sum over all values $y(m)$. In practice a few samples $y(m)$ around the value $n$ will be sufficient because the weighting factors in the sum decrease with increasing distance $|n-m|$.


Is there a linear filter whose step response is an S-curve?

Sums of harmonic cosines

The Hann window is a sum of a rectangular function and a truncated cosine. As an impulse response it has the qualities you asked for:

Figure 1. Hann window as impulse response, and the corresponding step response. The step response is valued $x - \frac{\sin(2\pi x)}{2\pi}$ over $0 \le x \le 1$, and is constant outside that range, because the impulse response is zero outside that range. The step response and its first two derivatives are continuous, because the impulse response and its first derivative are continuous.

A filter with Hann impulse response can be constructed as the (normalized) sum of a truncated integrator filter and the real part of the lowest-non-zero-frequency bin of a sliding discrete Fourier transform, both of which can be calculated recursively as truncated IIR filters.

The filters have poles on the z-plane unit circle (they are marginally stable) and require periodical refreshing of their state to cut accumulation of round-off errors. For fixed point input the truncated integrator can be implemented with extended bit depth so that no rounding errors accumulate and refreshing its state is not necessary.

Octave source for the composite filter:

r = 20;
w = 2*pi/r;
b0 = zeros(1, r+1);
b0(1) = 1;
b0(r+1) = -1;
a0 = [1, -1];
b1 = zeros(1, r+2);
b1(1) = -1; b1(2) = cos(w); b1(r+1) = 1; b1(r+2) = -cos(w);
a1 = [1, -2*cos(w), 1];
length = r+10;
impulse = eye(1, length);
ir = filter(b0, a0, impulse) + filter(b1, a1, impulse);
normalizationfactor = 1/sum(ir)
normalizationfactor = 1/r #Alternative calculation
ir *= normalizationfactor;
plot([0:length-1], ir, "o");

Figure 2. Hann impulse response from the first Octave script.

It is possible to get rid of the unnecessary one-sample delay without additional cost by multiplying the complex transfer function of the sliding DFT bin by $e^{-i\omega}$ before taking the real part:

r = 20;
w = 2*pi/r;
b0 = zeros(1, r+1);
b0(1) = 1;
b0(r+1) = -1;
a0 = [1, -1];
b1 = zeros(1, r+2);
b1(1) = -cos(w); b1(2) = 1; b1(r+1) = cos(w); b1(r+2) = -1;
a1 = [1, -2*cos(w), 1];
length = r+10;
impulse = eye(1, length);
ir = filter(b0, a0, impulse) + filter(b1, a1, impulse);
normalizationfactor = 1/sum(ir)
normalizationfactor = 1/r #Alternative calculation
ir *= normalizationfactor;
plot([0:length-1], ir, "o");

Figure 3. Undelayed Hann impulse response from the second Octave script.

There are other window functions with harmonic cosine terms, and filters with them as impulse responses could be implemented similarly.

B-splines

Piece-wise polynomials such as B-splines are another alternative for the impulse response shape. Cascaded integrator–comb (CIC) filters can be used to implement them. A triangle function would be the minimal piece-wise polynomial function (and B-spline) to meet your demands. The implementation is very simple: A truncated integrator in cascade with an identical truncated integrator. This comes from the favorable property of B-spline functions that they are each a convolution of an impulse by rectangular function multiple times. The number of continuous derivatives is incremented by each convolution. The earlier notes about stability apply.

Figure 4. A triangle impulse response and the corresponding step response.

Octave source:

r = 20;
b = zeros(1, r+1);
b(1) = 1;
b(r+1) = -1;
a = [1, -1];
length = 2*r+10;
impulse = eye(1, length);
ir = filter(b, a, filter(b, a, impulse));
normalizationfactor = 1/sum(ir)
normalizationfactor = 1/(r*r) #Alternative calculation
ir *= normalizationfactor;
plot([0:length-1], ir, "o");

Figure 4. Triangular impulse response from the Octave script.


Determining the pattern direction of a bitmap image

Solution is using structure tensors

from skimage.feature import structure_tensor

Axx, Axy, Ayy = structure_tensor(gray_image, sigma=1)
tan2t = (2*Axy)/(Axx-Ayy)
t=np.arctan(tan2t)/2
data = t.ravel()

# find the mode of `data` array (can use histogram). it is the direction.

another question and answers from other community


what should be the min number of samples and sampling frequency while playing sound

Your sampling frequency should be at least two times higher than the highest frequency in your signal. (Read about Nyquist frequency).
This has nothing to do with clipping.
I can imagine that your clipping occurs because you are adding multiple signals without scaling (normalizing) the amplitudes.


For transmission of audio by amplitude modulation of for example light, the audio can be made unipolar in order to make silence correspond to no transmission which for some receivers means no noise during silence. Just a guess but your audio may have gone through such processing. From the look of the waveform it may have worked like this:

Start with bias = 0. For each time step:

bias = bias + constant, where constant determines the slope seen during silent segments in your graph
If bias + input > 1, then set bias = 1 - input to avoid clipping.

It is not possible to exactly recover the original audio from the output of such processing, but you could use a similar method that keeps shifting the bias towards the direction that centers the audio, and changes the bias appropriately any time the audio would be clipped.


given a mixture of two signals, how to separate it to its base components?

I would suggest you take a look at 'independent component analysis', and 'blind source separation' tasks. In general, given a signal of mixtures of components to be separated to its contents well addressed by optimization in complicated cases. Here you show one example, but probably in each experiment you have something very different. I will also take a look at my files to see if still I have relevant codes to solve it or not.


Suppose you have a sinusoidal that has a whole number of cycles ($k$) in your DFT frame containing $N$ sample points.  It can be parameterized like this:

$$ x[n] = A \cos \left( \left( k\frac{2\pi}{N}\right)n + \phi \right) $$

If you take the $1/N$ normalized DFT of this (FFT is a DFT that is computed efficiently), all the bins will be zero except for bins $k$, and $(N-k)$.  With MATLAB, bin $k$ occurs at index $k+1$.

$$ X[k] = \frac{A}{2} e^{i\phi} $$

and

$$ X[N-k] = \frac{A}{2} e^{-i\phi} $$

So, you can see, in the ideal case of a pure tone with a whole number of cycles in the frame, the phase angle of the DFT bin corresponds directly to the phase argument in the signal.

The values from $-\pi$ to $\pi$ are by convention and are measured in radians.  This range covers every possible angle.

If you don't have a whole number of cycles, you can find my simplified bin value formulas here:  https://www.dsprelated.com/showarticle/771.php

https://gizmodo.com/pentagon-ordered-to-tell-congress-if-it-weaponized-tick-1836391549

In the time domain, a $2\pi$ change in the phase represents the shift of an entire cycle, which can also be considered one wavelength.  Shift it by $\pi$ (popularly known as a 180 degree shift), and you effectively invert it.  Shift it back and forth by $\pi/2$ and sine becomes cosine and vice versa.

Maybe this will spark an understanding:

$$ cos( ft + \phi ) = \cos( ft ) \cos( \phi ) - \sin( ft ) \sin( \phi ) = a \cos(ft) + b \sin(ft) $$

So, fiddling with the phase adjusts how much cosine vs the sine is in the tone within that reference frame.

What I described above is the bedrock connection between the phase in the time domain and the phase in a DFT bin for real valued signals.

The relationship between this shifting in the time domain and the angle in the corresponding bin is one-to-one for sinusoidals with a whole number of cycles in a DFT frame.

$$ a = \cos(\phi) $$

$$ b = -\sin(\phi) $$

$$ \frac{b}{a} = -\frac{\sin(\phi)}{\cos(\phi)} = -\tan(\phi) $$

Mix in an $i^2$ and you got bin interpretation.


How to calculate the plane at infinity? either in Matlab or Python

The plane at infinity has to be calculated as following:
function plane = computePlaneAtInfinity(P, K)

%Input
%   P - Projection matrices
%   K - Approximate Values of Intrinsics
%
%Output
%   plane - coordinate of plane at infinity

% Compute the DIAC W^-1
W_invert = K * K';

% Construct Symbolic Variables to Solve for Plane at Infinity
% X,Y,Z is the coordinate of plane at infinity
% XX is the scale
X = sym('X', 'real');
Y = sym('Y', 'real');
Z = sym('Z', 'real');
XX = sym('XX', 'real');

% Define Normal to Plane at Infinity
N = [X; Y; Z];

% Equation of Dual Absolute Quadric (DAQ)
Q = [W_invert, (W_invert * N); (N' * W_invert), (N' * W_invert * N)];

% Select Any One Projection Matrix
M = P(:, :, 2);

% Left hand side of the equation
LHS = M * Q * M';

% Solve for [X, Y, Z] considering the System of Linear Equations
% We need 4 equations for 4 variables X,Y,Z,XX
S = solve(LHS(1, 1) == (XX * W_invert(1, 1)), ...
          LHS(1, 2) == (XX * W_invert(1, 2)), ...
          LHS(1, 3) == (XX * W_invert(1, 3)), ...
          LHS(2, 2) == (XX * W_invert(2, 2)));

plane = [double(S.X(1)); double(S.Y(1)); double(S.Z(1))];

end
```


Bandpass Filter Order vs Sample Numbers of Noisy Signal

You are dealing with a time frequency uncertainty in estimation problem.  It takes longer than 3 seconds for there to be enough information in a signal to tell the difference between 3 clearly separated sinewaves between 49 and 51 Hz, and only 1 fat (slightly modulated) one.  And you can't filter out a signal unless there is enough information to tell whether it's even there or not.


Detecting sloppy hand-drawn arrows in images

Well, the easiest way to approach this problem would be to obtain the Hough Transform and look for the "signature" of the arrows. In addition and because these lines are drawn on a piece of paper, some thinning might be required because the trace left by the pen will not be a perfect straight line but rather a set of somewhat straight, parallel lines.

OK, but what does all this mean?

The Hough Transform obtains the sum of all pixels in the original image along a profile ($\rho$) set at some angle ($\theta$). For instance, to produce the sum along the top-bottom direction, the algorithm would first add all the pixel values of the first column and produce one value. Then it would sum all values of the second column and produce yet another value and so on until a sum for each column has been produced. Each column's sum now represents one line of data in the transform's output image. This line of data corresponds to some angle ($\theta$). Then, the algorithm rotates the image by some $d\theta$ and repeats the process thus creating a second line of data at some angle $\theta + d\theta$ and so on. This is how this image is produced.

So, the end result of the Hough Transform is to transform lines to points.

Therefore, looking for an arrow now becomes the task of finding 1 "strong" pixel value at some angle $\theta$ and displacement $\rho$ (the middle line of the arrow) that denotes the angle that the arrow is pointing and 2 less "strong" pixel values at angles (approximately) $\theta \pm 45$ (these are the adjacent lines that form the point of the arrow). So, the Hough Transform of an image with 3 arrows is an image with 3 distinct triplets of points.

BUT! A scanned image of those arrows, is not going to be composed of just three lines. Having hand-drawn arrows means having imperfections. These imperfections "spoil" the output of the Hough Transform and instead of a sharp point where there is a line, what you get is a nice fat blob.

This is where thinning will help. Thinning is a morphological operator that will reduce a "fat" line to its thinest representation which will help a lot later on with the Hough Transform.

With this technique, you would have to scan the output of the Hough Transform once to find all the arrows in the image.

Please note that this assumes that the arrows are not co-linear. If you have more than one arrows pointing towards the same direction, along the same line, they will be contributing to the sum of the same accumulator and you would have to make assumptions on the length of each arrow to distinguish how many arrows are there along a line from its pixel value.

The Hough Transform is all about polar coordinates...and summation. The thinning is an improvement.

I guess that there might be ad-hoc ways of scanning the image until an "ink" area is found and then tracing that area with a set of constraints to decide if it is an arrow and which direction it is pointing towards but these would not be as straightforward as the Hough Transform approach.


Calculating Signal Approximation Error (Theoretical signal - Real hardware signal)?

There are lots of different error statistics you depending on what you want to know. If what you care about is the maximum error then what you suggest is fine.

Probably more common is to use the RMS error i.e. $ \sqrt{\frac{1}{N}\sum(Y-Y_{real})^2} $. This is related to the variance of the system. I'm not sure if this is exactly what you want for your case as I suspect there will be strong spatial correlation across the filter. This is something you should look at.

Another issue is that the error is dependent on the input signal. You may wish to look at the frequency response of both systems to try and remove this issue.


How to calculate power of a signal at one specific frequency?

Ultimately if your noise is white then the matched filter would be the best approach (multiply by one frequency and sum-- if it is complex you would multiply by the complex conjugate- both a form of correlation). This is equivalent to computing a single bin in the DFT (assuming your frequency is on bin center, otherwise equivalently the operations in the DFT can be modified to be centered on your particular frequency). This approach has bandwidth which for white noise will have the equivalent power of a similar brickwall filter with a bandwidth of $1/T$ where $T$ is the duration of your correlation; but importantly the frequency response is a Sinc function (or Dirichlet Kernel if discrete as a sampled Sinc) in frequency which means this approach is sensitive to non-white noise at other frequency locations (a strong tone very relatively far away would show up in your power measurement); the main lobes of the Sinc only roll-off at 6 dB per octave in frequency.  So for white noise, this approach is optimum, for other applications with colored noise, applying a properly designed filter to the signal to attenuate "out of band" noise before measuring the power would be recommended.
For measuring a single tone with the equivalent of a single DFT bin, the Goertzel is a well known choice that is efficient at doing this very thing.


A few questions about the output value of dsp.SpectrumAnalyzer

I don't really know what Matlab's spectrumanalyzer does internally, but I suspect it mimics the behaviour of an actual spectrum analyzer – which does not at all do a DFT and then square the output.

Instead, a SA shifts a filter across the spectrum. This filter's equivalent noise bandwidth is exactly how much white noise power gets let through and gets measured to derive a single point in that display. So, for white noise, there's always the factor "filter bandwidth" between your power spectral density and the power displayed in your plot.

Also, filters don't have to be unitary, i.e. there might be yet another factor, the filter gain, between your sine's power and the displayed power.

Also, in a SA you typically need to configure the "frequency advance step" separately from the filter bandwidth, hereby defining how much the windows "overlap". A quick skim of the matlab doc says the default overlap is 0 – so you might really just be missing your sine wave's "peak".

All in all, it's quite questionable why you'd want to use this visualization if you could also directly estimate the spectrum all at once, without losing any info, using any of the PSD estimators that Matlab can do – from the simple (windowed) FFT estimator, to the more complex parametric estimators.

This tool definitely only has the purpose of illustrating how measurements on a Scope can be modeled and how they differ from the mathematical idea of the PSD. Unless you need to model a SA, I really don't see any practical value in this tool.


Brightness equalization between different images

You could do Histogram equalization on the histogram of brightness values across all tiles.

Alternatively, you can loop over all tiles, accumulate the minimum and maximum brightness value of each tile, and then scale all brightness values in all tiles against the overall minimum/maximum.


You are correct that the total accumulated phase offset due to Doppler could be calculated by integrating the Doppler shift $\Delta f(t)$ over time. In addition to Doppler shift, however, you would also need to take into account any frequency difference between the transmitter and receiver's reference oscillators.

However, I would question the value of trying to do so. In most communications systems, there will be a random, unknown phase offset between the transmitter and receiver oscillators. This makes absolute phase measurements meaningless unless you have some reference. In real systems, phase synchronization is either achieved by exploiting some known characteristic of the signal (e.g. the known constellation of a PSK signal), or in some rare cases, via the explicit synchronization of the transmitter and receiver local oscillators.


Demonstrate the mean of the sample variance

The key to prove
$$
E\left[\widehat{\sigma_x^2}\right] = \sigma_x^2 \left(1 - \dfrac{1}{N} \right)
$$
is to add and subtract $\mu_x$ in the expression for the sample variance before working out the square. A detailed proof is available on the Wikipedia page for Bessel's correction.


2D Multi Scale Dot Enhancement Filter based on Gaussian Filter and Hessian Matrix

Questioner's answer...

sigma have the same units as x and y i.e. number of pixels.
In multi-scale filtering, the size of the filter must change when the sigma changes.
Obtain the number of pixels per one millimeter or the vice-versa. (I did this using the property of pixel spacing included in the DICOM metadata in Matlab you can do this as info=dicominfo('image.dcm'); and spacing=info.PixelSpacing;)

What I was doing wrong is not changing the filter size when I changed sigma in the different scales.

So, the solution to my problem is like this:
the size of the gaussian kernel should be 2 times or preferably 3 times the value of sigma on eather side of the origin as Mr. Cris suggested.
Gaussian filtering with Matlab's Image Processing Toolbox

cutoff = ceil(3*sigma);

In matlab there is two options, using fspecial() or imgaussfilt(), fspecial() is not recommended anymore in the newer versions of Matlab and latter is the recommended one.

h = fspecial('gaussian',2*cutoff+1,sigma);
B = conv2(A,h,'same');

or...

B = imgaussfilt(A, sigma, 'FilterSize', 2*cutoff+1);

Another thing that I did that corrected my results is converting my input image into type double using im2double() function.

Here is the result:


Why do we need to conjugate complex signals in autocorrelation and cross correlation

$\underline{Prologue :}$

Let me ask you another question, How will you compare two complex numbers $U$ (a+jb) and $V$ (c+jd)?
By comparing magnitude? Subtract them and take real part? Multiply them and compare?

Since any complex number involves two entities ( one for magnitude |$z$| and other for argument $\theta$ ) any comparison involves comparison of those two entities.
Suppose we represent $U$ and $V$ in a complex argument plane as $M e^{j \theta}$ and $N e^{j \phi  }$  if we multiply both of them together, we get

$U$ x $V =  M e^{j \theta}   *  N e^{j \phi } = MN e^{j ( \theta +  \phi )}  $

If we seek maximum value out of the above multiplication then

$|U$x$V|_{max} =  |MN|_{max} * |e^{j ( \theta +  \phi )}|_{max} $

Suppose our comparison is only about the relative location of $V$ with respect to $U$ in the complex argument plane, for the time being, forget about the magnitudes $M$ and $N$.

maximum value of the second term $\;e ^{j ( \theta +  \phi )}$  occurs at 1.

for that to happen
$\theta+\phi = 2\pi k \;\;\;\;$ where $\:k\:\epsilon\;( 0,1,2...$

if $k = 0$,

$\theta + \phi = 0$

$\phi =  - \theta $

ie, when two complex numbers with bounded magnitudes get multiplied, their maximum value occurs, when the second number ($V$) is in the negative argument direction of the first number ($U$), ie in the direction of its conjugate ($U^*$).

ie
maximum of $U \times V$ occurs when   $V$ points in the  direction of $U^*$ { the orange dashed line in the picture }

$\textit{-or conversely-}$

If $V$ is getting multiplied by $U^*$ then their maximum occurs when  $V \textit{points in the direction of } U$

In another words $U^* \times V$   gives nearness of $V$ w.r.t $U$ in terms of angle ( complex argument ) and it will decrease as a function of $cos \:\alpha$, where $ \alpha$ is the angle between $U$ and $V$.

$\underline{Epilogue:}$

Considering two signals $f$ and $g$, correlation is given as

$\hspace{1cm} (f \otimes g)(\tau) = \int_{-\infty}^\infty f^*(t) g(t + \tau) dt$

At any instant '$t$', a signal is just a point in the complex plane. So at any moment $f$ and $g$ are just two points in the complex plane. Then the job of comparison of two signal becomes a mere comparison of two complex points.
With this logic, to compare two signals with a given lag '$\tau$' between them, just multiply point by point in '$t$' by taking conjugate of the other signal. If you get a big number that means both reference and compared signals are just looking in the same direction in the complex plane.
When integrated on the whole '$t$', the value of the integration shows the how much similarity is between signal $f$ and signal $g$.
that is which is known as the cross-correlation between two signals.

{ I didn't talk about magnitudes $M$ and $N$, right? The method in which everything is put in a single benchmark by avoiding magnitudes $M$ and $N$ is known as Normalized cross-correlation }


Multiple cheap IMUs vs one expensive IMUs

That really will depend on the variance (standard deviation) of the measurements you get from the five vs the one.

Suppose your expensive IMU gives a measurement distributed as:
$$
x_\mbox{expensive}  \sim  N\left(x_\mbox{truth}, \sigma^2_\mbox{expensive}\right)
$$
and suppose your cheaper ones give measurements distributed as:
$$
x_\mbox{cheap}^i  \sim  N\left(x_\mbox{truth}, \sigma^2_\mbox{cheap}\right)
$$
for $i=1\ldots 5$ and that you form the average of these cheaper measurements to get the "real" measurement:
$$
\hat{x} = \frac{1}{5} \sum_{i=1}^5 x_\mbox{cheap}^i
$$
then $\hat{x}$ is distributed as:
$$
\hat{x} \sim  N\left(x_\mbox{truth}, \frac{\sigma^2_\mbox{cheap}}{5}\right)
$$
So if
$$
\frac{\sigma^2_\mbox{cheap}}{5} \le \sigma^2_\mbox{expensive}
$$
you should be OK with the cheaper sensors.

Of course, there may be other problems with your cheaper sensors: they may be biassed (so that the mean of their distributions will all be different from $x_\mbox{truth}$ ).

And I'm not quite sure what you mean by "voting"? Do you vote off the one furthest from the mean of the 5?  What if the mean changes for the remaining four?  I'll do a bit more digging to see if there's a viable algorithm there, but the direct mean of all five might be the easiest.


Does negative frequency actually exist or it is just theoretical?

Consider a wheel rotating counter-clockwise at one revolution per second. Its frequency of rotation is 1 Hz. If it rotates clockwise, its frequency of rotation is -1 Hz. It's as simple as that.


Minimum size for positive samples for HOG (Histogram of Gradient) detector training

You should use elements which are a power of 2, for easy and faster computation.
I'll say, start with 8x8, so you can only scale up.
And if you have to look only at some portions of this image, then you can crop the original one.


Autocorrelation of a uniform random process

As suggested I am adding my comments as an answer.
The sample autocorrelation function (ACF) for $n$ observations is given by

$\hat{p}_x(m) = \frac{\sum_{n=0}^{N-m-1}(x_n - \bar{x})(x_{n+m} - \bar{x})}{\sum_{n=0}^{N-m-1}(x_n - \bar{x})^2} $

To understand why you have a difference between figure 1 and figure 2, lets assume we are looking at observations $x_n$ being zero mean. Let be $y_n = x_n +c$ with c being a constant. Looking now at matlabs help for xcorr for the input $y_n$ we get

$r(m)=\sum_{n=0}^{N-m-1}y_n y_{n+m}^* = \sum_{n=0}^{N-m-1}(x_n + c)(x_{n+m} + c)$

which can be summarized to

$\sum_{n=0}^{N-m-1}(x_nx_{n+m}) + (N-m)c^2$

In the case of xcorr $m$ equals $m=-N,...,N$. Which means you are adding a triangle to your result.

To calculate the ACF matlab has the autocorr function, which takes care of normalization and demeaning of the input.


Are there any filters that look (vaguely) like this?

I would approximate the filter using an finite impulse design where the wrights of the FIR filter are computed based on the spectrum that the filter passes.  You can do this by choosing samples of the desired frequency response (from the graphs above) and taking the inverse discrete Fourier transform of the samples.  The result will be the impulse response of the desired filter.  The filter is then realized as a weighted sum using the impulse response sample values as the weights.  So your filter will be $A_0 + A_1Z^{-1} + A_2Z^{-2}...$  Where $A_0$ is the first value of the impulse response, $A_1$ is the second and so on.  $Z^{-n}$ represents a sample delay of n (coming from a Z domain representation of the filter structure).  This approach is sometimes referred to as the "frequency sampling" approach to FIR filter design.  I have used this approach to create time sequences that contain specific spectral content.  Google FIR filter design, frequency sampling.