Obtaining the system equations

Using Kirchhoff's voltage law:
$$
\begin{align}
  V(t) &= V_{R_1}(t) + V_{R_2}(t) + V_{R_3}(t) + V_I(t) \tag{1}\\
  V_{R_1}(t) &= V_{C_1}(t) \tag{2} \\
  V_{R_2}(t) &= V_{C_2}(t) \tag{3}
\end{align}
$$
And similarly using Kirchhof's current law:
\begin{align}
  I_{R_3}(t) &= I_{R_1}(t) + I_{C_1}(t) \tag{4}\\
             &= I_{R_2}(t) + I_{C_2}(t) \tag{5}
\end{align}

Now using the Voltage-Current relationships of the resistances and capacitors (assuming ideal components):

\begin{align}
  V_{R_1}(t) &= R_1 I_{R_1}(t) \tag{6}\\
  V_{R_2}(t) &= R_2 I_{R_2}(t) \tag{7}\\
  V_{R_3}(t) &= R_3 I_{R_3}(t) \tag{8}\\
  I_{C_1}(t) &= C_1 \frac{d\,V_{C_1}(t)}{dt} \tag{9}\\
  I_{C_2}(t) &= C_2 \frac{d\,V_{C_2}(t)}{dt} \tag{10}
\end{align}

Based on your comments, a current sink with a resistive internal resistance can be assumed. For analysis, it is common to replace such a load with its Norton equivalent circuit (see Norton theorem) which consists of putting the internal resistance (say $R_4$) accross an ideal current sink $I(t)$.
\begin{align}
  V_I(t) &= R_4 I_{R_4}(t) \tag{11} \\
         &= R_4 (I_{R_3}(t) - I_I(t)) \tag{12}
\end{align}

Numerical solution

To discretize the problem to discrete time step $t_n = n \Delta t = n/f_s$, we can make approximations for $(9)$ and $(10)$:
\begin{align*}
  I_{C_1}[n]
    &\approx C_1 \frac{V_{C_1}[n+1] - V_{C_1}[n]}{\Delta t} \\
    &\approx C_1 f_s \left(V_{R_1}[n+1] - V_{R_1}[n]\right) \\
  I_{C_2}[n] &\approx C_2 \frac{V_{C_2}[n+1] - V_{C_2}[n]}{\Delta t} \\
    &\approx C_1 f_s \left(V_{R_1}[n+1] - V_{R_1}[n]\right)
\end{align*}
This gives us the recurrence equations:
\begin{align*}
  V_{R_1}[n+1]
    &\approx V_{R_1}[n] + \frac{1}{C_1 f_s} I_{C_1}[n] \\
    &\approx V_{R_1}[n] + \frac{1}{C_1 f_s} \left(I_{R_3}[n] - I_{R_1}\right) \\
    &\approx V_{R_1}[n] + \frac{1}{C_1 f_s} \left(\frac{V_{R_3}[n]}{R_3} - \frac{V_{R_1}[n]}{R_1}\right) \\
  V_{R_2}[n+1]
    &\approx V_{R_2}[n] + \frac{1}{C_2 f_s} I_{C_2}[n] \\
    &\approx V_{R_1}[n] + \frac{1}{C_2 f_s} \left(I_{R_3}[n] - I_{R_2}\right) \\
    &\approx V_{R_1}[n] + \frac{1}{C_2 f_s} \left(\frac{V_{R_3}[n]}{R_3} - \frac{V_{R_2}[n]}{R_2}\right)
\end{align*}
And $V_{R_3}[n]$ begin obtained by substituting $(8)$ and $(12)$ into $(1)$:
\begin{align*}
  V_{R_3}[n]
    &= \frac{V[n] - V_{R_1}[n] - V_{R_2}[n] + R_4 I[n]}{1+R_4/R_3}
\end{align*}

which can readily be implemented in python (with $V(t)$ and $I_I(t)$ being constant for illustration, but the extension to arbitrary function is trivial) with:

# numerical solution
T    = 1.0/fs
Tmax = 20
t = np.arange(0,Tmax,T)

vr3 = np.zeros(len(t))
vr1 = 0
vr2 = 0
for i in np.arange(len(t)-1):
    vr3[i] = (V-vr1-vr2+R4*I)/(1+R4/R3)
    i3  = vr3[i]/R3
    ic1 = i3 - vr1/R1
    ic2 = i3 - vr2/R2
    vr1 = vr1 + ic1/(C1*fs)
    vr2 = vr2 + ic2/(C2*fs)
vr3[-1] = (V-vr1-vr2+R4*I)/(1+R4/R3)

Steady-state solution

As a validation of the above numerical solution, we can check whether it converges to a steady state which can be computed analytically fairly easily (assuming a steady state can be reached, which typically requires $V(t)$ and $I_I(t)$ to be near constant). To do this we note that in the steady state derivative terms go to zero. As a result, $I_{C_1}(t) = I_{C_2}(t) = 0$ and consequently, $I_{R_1}(t)=I_{R_2}(t)=I_{R_3}(t)$. Thus $(1)$ can be written as:
\begin{align*}
  V(t) &= R_1 I_{R_1}(t) + R_2 I_{R_1}(t) + R_3 I_{R_1}(t) + R_4\left(I_{R_1}(t) - I_I(t)\right) \\
     &= \left(R_1+R_2+R_3+R_4\right) I_{R_1}(t) - R_4 I_I(t)
\end{align*}
From which we can solve for $V_{R_3}(t)$:
\begin{align*}
  V_{R_3}(t) &= \frac{R_3}{R_1+R_2+R_3+R_4} \left(V(t) + R_4 I_I(t)\right)
\end{align*}

For illustration, the graph below shows an example for the following parameters:

# parameters
R1 = 2200.0
C1 = 0.001
R2 = 1000.0
C2 = 0.0047
R3 = 5600.0
R4 = 100000.0
V  = 5
I  = 1e-5
fs = 10

As a final note, for voltage source $V(t)$ and current sink $I_I(t)$ which are constant or are from a set of fairly simple functions, the solution could also be obtained analytically with the use of Laplace transforms.


How can spectrum of DTFT of sampled band-limited analog signal be in $[-\pi,\pi]$

The DTFT $X_d(\omega)$ of a discrete-time signal is periodic with period $2\pi$, so it's sufficient to consider the interval $[-\pi,\pi]$. "Restricted" to this interval means that due to the inherent periodicity there can be no new frequency components outside the interval $[-\pi,\pi]$. Of course, $\pi$ corresponds to half the sampling frequency.

If $X_a(\Omega)$ is the spectrum of a band-limited continuous-time signal, and if the sampling frequency is greater than twice the band-width, then the equality

$$X_d(\omega)=\frac{1}{T}X_a\left(\frac{\omega}{T}\right)$$

holds in the interval $[-\pi,\pi]$. Outside that interval, $X_a(\Omega)$ is zero, and $X_d(\omega)$ is continued periodically.


Understanding of wireless communication system from MATLAB point of view

Below is a sample MATLAB/pseudocode to generate a 4x4 MIMO system, for 1000 Monte Carlo trials:

M=4; MC=1000;
H= randn(M,M,MC);
x=randn(M,1,MC);
w=randn(M,1,MC);
Y=zeros(M,1,MC);
for j = 1:MC
   Y(:,:,j)= H(:,:,j)*x(:,:,j) + w(:,:,j);
end

This would generate y according to the equation in the question. You can view the MC values as different time instances. Meaning different channel, transmit vector and noise matrix for each instance.


Gradients enhances differences. Hence, highly non-smooth padding (like zero-padding) yields higher gradients, as they don't ensure continuity. Symmetric extensions preserve continuity, and anti-symmetric ones are likely to better preserve derivatives. There is a lot of literature on that, for instance in the filter bank domain. Suppose that the edge sample is $x[0]$, you can have whole-sample anti-symmetry or half-sample anti-symmetry:

$x[-1] = -x[1]+2x[0] $
$x[-1] = x[0] $ and $x[-2] = -x[1]+2x[0] $

This is illustrated in the picture, from the first paper mentioned below. One can find whole-sample (WS) and half-sample (HS) symmetries and anti-symmetries. A lot of image processing operators combine symmetry or anti-symmetry, so those can be useful.

A couple of pointers:

Lifted linear phase filter banks and the polyphase-with-advance representation (JPEG 2000 image coding applications)
Classification of Nonexpansive Symmetric Extension Transforms for Multirate Filter Banks
Unitary FIR filter banks and symmetry

Recently, unitary FIR filter banks with linear-phase have been
completely parameterized by exploiting the eigenstructure of the
exchange matrix. This correspondence gives new characterizations of
polyphase matrices of filter banks with various types of symmetry on
the filters. Using matrix extensions of the well-known hyperbolic and
orthogonal lattices, we give a alternative proof for the
parameterization of linear-phase unitary filter banks. A complete
parameterization of FIR unitary filter banks with each of the
different types of symmetries considered (not just linear-phase) is
also given. These results can also be used to generate non-unitary
filter banks with symmetries, though no completeness results can be
obtained. In some cases implicit, and in others explicit
parameterization of wavelet tight frames associated with these filter
banks are also given. This paper only considers filter banks with an
even number of channels. A similar theory can be developed if the
number of channels is an odd integer.

Aside, there are many other possible extensions, using constant continuation (or pixel replication), or higher-order (eg linear, polynomial) extrapolations:

Image padding, REFLECTION padding


From fft result one can affirm that the signal is affected by white noise?

These are just some hints on how to proceed which are too long for a comment. We cannot give a definite answer from just looking at the plots.

First, you should plot the magnitude of the FFT in dB, i.e. plot(20*log10(fft(signal))). Then, you would not need to plot different zoomed versions.

To refer to question 1) Is this the whole signal or just a part of it? Try analyzing a longer measurement cycle to become more certain. Certainly, there is some addition of white noise to the signal. But, additionally, it seems to me there is an additional sine-wave of period length 4 contained (referring to your first plot).

To do some more analysis:

Assuming you measure a constant acceleration, remove the mean from the signal. Then, the remaining signal just contains the measurement error.
Have a look at the (circular) autocorrelation of this error. In case of white noise, it should be ideally a single peak at delay 0. If the autocorrelation exhibits extra (significant) peaks, there probably is something extra on top of the white noise. (Looking at your FFT-plots, I actually expect some extra peaks).

To answer questions 2), 3) (I dont really understand your question 2)):

If your sampling frequency is 56Hz, your FFT output will range from 0 to 56Hz, where the distance between the bins is 56Hz/FFT-length. Furthermore, as your input signal is real, the frequencies above 56/2Hz will be a mirror of the first half of the frequency domain output, they do not contain additional information. As a general rule for Fourier Analysis: The frequency resolution $\Delta_f$ after FFT is only determined by the time duration $T$ of the original signal, given by $\Delta_f=1/T$. You can have a look at this article I wrote on my website about the FFT and frequency domain resolution as an introduction to the topic.


Real-time oscillation-like noise filtering

Marcus and Ben have very clearly identified the problem: the signal is simply too small relative to the least significant bits (LSBs) of the "virtual machine"'s ADC. The sensible solution is to increase the signal by use of gain. From the provided plots, a signal gain around 100 would be a good starting value. Alternatively, the ADC resolution would have to be increased (> 4 more bits), so the LSBs would be smaller. In a real system, the first option would be the only feasible option, since amplifiers are relatively inexpensive and high resolution ADCs are most certainly not. Since the goal is real-time noise filtering, rather than Lazarus filtering, i.e., post-mortem reversal of fatal mis-quantization, dithering and signal averaging is ruled out.


Summation or integration of the spectral components yields the total power (for a physical process)

yes, that's how "power density" is defined ...

  or variance (in a statistical process)

For any bound density, that's true.

  But would the total power correspond to the uncentralized second moment, so variance + mean²?

Yes, but if you have a process with a non-zero mean, that would imply that the spectral power density at $f=0$ diverges (a "Dirac delta $\delta(f)$" with "infinite height"), and that would imply the PSD isn't actually a function (it's still a distribution / test function, but seriously, this is just getting into math).

If you will, just imagine your "mean" is something at very low frequencies ($f\in[-\epsilon;\epsilon]$), but not concentrated to exactly $f=0$.

Still works out the same, for anything technical (because nothing will ever experience whatever process they're observing for eternity, so who's to say whether something happens at $0\,\text{Hz}$ or $10^{-24}\,\text{Hz}$?).

Honestly, the wording of the article isn't the greatest.


Synthetic sounds to describe motion in 3d or 100d

You say your points are strictly ordered:

  a sequence of points $X_1 \ X_2 \ X_3\  \ldots\ X_N$

Well, in that case, you can

either define "speed" as "euclidean distance (or any other metric) per in between two consecutive points":
$$s_n = ||X_{n+1} - X_n||,\quad n=\{1,\ldots,N-1\}$$
and interpolate a steady function $s(t)$ from that sequence, or
you can define speed as the first derivate of an interpolation function through all $X_i$. However, finding an interpolator for $N$ points in 100 dimensions isn't easy, and also, its utility would be pretty questionable, unless you have some goal in mind that clearly says why a specific interpolator is appropriate for your problem.

After you have a scalar value along a single axis, well, map to tone or loudness.

I really doubt the "illustrating" effect is overly great here: the choice of your interpolation will probably have the most important effect on what listeners will take as an impression from this, and that choice is not backed by the data itself – it's a model choice you make. Which would be ok, if you know something about the data (i.e. all the individual dimensions just contain harmonic oscillations would suggest classical low pass filtering would be an appropriate method to interpolate "speed" points), but if you don't and want to use this "audiolisation" to understand your data better, you'd end up fooling yourself. You will listen to your own assumptions, not the data.

However, you mentioned PCA, which indicates the points aren't actually to be understood sequentially! That of course makes things a lot harder. You'd essentially be looking for a $100-1$ hypersurface that describes how your points are mainly distributed in your $100$ dimensional space, and then look for the main gradient of that – and that leads to all the nice problems that you get with higher dimensional analysis (lack of "one" derivate of function, mapping values from something that isn't a vector space to a vector space, which introduces the need for maps on manifolds... meh).


Python bandpass filter - singular matrix error

I'd suggest checking again your parameter choices for

b,a = fd.iirdesign(Wp,Ws,Rp,As,ftype='butter')

Looking at the resulting filter coefficients, you have pole/zero cancellation at DC, i.e. at $z=1$.


Applying the CUSUM algorithm to a correlated random process

If you look at the derivation of CUSUM from first principles, the idea is to construct a likelihood ratio test statistic and wait till it exceeds a threshold. In case of CUSUM, the assumption is that samples are uncorrelated and so the test statistic simplifies to a sum. Since you are dealing with a different statistical change model, you can try writing a likelihood ratio statistic and see if you can simplify it to anything useful. There are also other algorithms out there for testing "shifts in variance" or "shifts in spectral content" for ARMA processes which you may want to read about.

See also: Ch 6 and 7 Detection of abrupt changes: theory and application
M Basseville, IV Nikiforov - 1993.


Segmenting a frequency series: do we need band-pass filtering?

I haven't heard of the second spectrum before, but yes, you can slice out a region in the frequency domain just like that.  This is equivalent to a maximal-length bandpass filter in the time domain and is in fact the best filter you can implement for a given FFT length.


Combining compressed measurements from the same source

If you have $k$ individual measured vectors $\mathbf y$, each taken with an individual measurement matrix $\mathbf \Psi$, you actually have an interesting, third option:
$$\mathbf y = \begin{bmatrix} \mathbf y_1 \newline \vdots \newline \mathbf y_k \newline \vdots \newline \mathbf{y}_K \end{bmatrix} = \begin{bmatrix} \mathbf \Psi_1 \newline \vdots \newline \mathbf \Psi_k \newline \vdots \newline \mathbf \Psi_K \end{bmatrix} \mathbf x$$
So now you effectively have a (much) taller measurement matrix, so you are effectively "under-sampling less" and that will improve your estimate $\mathbf x$ based on the above equation. How much, I cannot recall an expression for off the top of my head.

Note also that you may have a tall measurement matrix now due to the stacking of the $\mathbf \Psi_k$'s. That is, you do not necessarily have an under-determined system anymore. This gives you the option of for example using least-squares estimation instead of $\ell_1$-norm optimisation. If you know that $\mathbf x$ is sparse, then it still makes sense to use sparse estimation ($\ell_1$-norm optimisation etc.), but if $\mathbf x$ is only approximately sparse, compressible then I would try least squares estimation as well.


What function to apply to a signal to expose its burstiness (i.e., make the burst part more evident than the non-bursty part)

You can apply a simple moving average filter and based on the result of the filter make a decision if you are in the sparse area or dense area, as your average will certainly be larger in your dense area. There is a resolution vs false alarm trade space that you will need to make; your result will be more accurate the longer you can make your moving average filter but your precision in determining where the transition occurred wil be reduced.

So in this sense, to compare to what you tried, instead of determining how far each signal is from the mean, take the mean over shorter intervals (as you move through the sequence), and based on the result of the mean compared to a threshold, decide if you are in the sparse or dense area, and if your decision is sparse, you can clip the result as you have done.

An alternate and simple algorithm would be a window counter, where you allow a certain number of counts (window) between samples above a certain threshold: if the sample comes in above threshold before the count is up you declare a positive, and if not you declare a negative and the count in reset. Further, a certain number of consecutive positives or negatives is required to change state between dense and sparse (this adds hysteresis and avoids a chatter). I don't know the detailed characteristic of your dense to suggest actual count intervals, but this is something you can consider playing with.


Matlab IIR cheby2 bandpass, problems with Magnitude

Actually we are also facing the same problem, instead of (Z,P,k), even if you use (A,B,C,D), as soon as sos will come in syntax with g, your magnitude response will scaled by the attenuation factor. If you are concerned only with plotting magnitude response you can plot as freqz(sos) you will get exact magnitude response without scaling.


How I did it was: Use the original data word and fill in the new check bits. :

001111001101

Then again get the positions where the 1s are and convert that number to binary:

1 : 0 0 0 1

3 : 0 0 1 1

4 : 0 1 0 0

7 : 0 1 1 1

8 : 1 0 0 0

9 : 1 0 0 1

10: 1 0 1 0

1 0 1 0

So this means that the 10th bit is wrong

Since I don't have the actual  solution, i do not know if this is correct, bit it fits the way I learned how to do this.


Overview of regression methods between the data

There are multiple ways to recover the signal from noise. As you have suggested one such method, the low pass filtering which is a classical signal processing method.

The linear signal model that you have considered is $y[k] = f[k] + b + n[k]$, where
$y[k] =$ output signal, $f[k] = $ true signal, $b=$ noise offset and $n[k] = $  zero mean additive noise

If you want to look from a data point of view the methods of interest would be

Nearest Neighbor Approximation (or local average) : This method is similar to low pass filter but here you can choose a suitable window, $$\widehat{f}_{NN}[k] = \frac{1}{2L+1}\sum_{i=k-L}^{k+L} y[i] - \hat b$$ Here $\hat b$ can be obtained by an initial pilot sequence of the zero signal.
Basis Representation followed by Ridge Regression : With suitable signal class assumptions you can represent $f[k] = w_1 \phi_1[k] + \ldots + w_p \phi_p[k]$ - known as the basis expansion (In the example you have described a triangular wavelet basis would work). In this setting the goal is to learn the weight parameters $w_i, \ldots, w_p$ subject to the prior knowledge on noise, that is the mean and variance of $n[k]$ is known apriori. In case the prior information are unknown you can use the method used to estimate $b$, explained in method 1.
Alternatively, There is an equivalent formulation in convex optimization known as ridge regression which allows for regularization cost and no assumption on the noise parameters.

These methods work out mostly with a training dataset and hence mostly does non-linear (filtering in) signal processing.

Remark and Disclaimer : Please feel free to correct any mistakes in the arguments presented above, as I am an amateur in the area of modern methods of signal recovery.


Your block diagram is correct. However, it seems to me that your understanding of how to implement it is not. There are no 'bits' in the modulated signal! The input signal is a noisy analog bandpass signal, and the I and Q components are derived from this signal as shown in the block diagram: by demodulating and filtering the incoming signal. Then this (complex) signal is sampled at the symbol rate (you may need timing recovery to know when to sample) and fed to a slicer. The slicer decides which of the four possible symbols was (probably) transmitted in each symbol interval. After the detected symbols are known you simply decode each symbol according to the coding in the transmitter (each symbol carries two bits).

EDIT: I didn't mention matched filters even though they were emphasized in Dilip's comment. For an optimum receiver they are of course necessary, but given the level of the question it seems to me that the block diagram is OK, and that the idea is probably that it should be implemented as is. Even though it will result in a sub-optimum receiver it can perform well under mild channel conditions. The only real difference is the type of low-pass filter used (matched to the transmitted pulses or just eliminating out-of-band noise).


C/D is simply a "continuous-to-discrete" time conversion. In other words, sampling! :-) You can use $T$ to sample your continuous time signal $x(t)$ into a discrete time signal $x[n]$. Just replace $t=nT$ in your $x(t)$ equation, to turn it into a discrete time signal $x[n]$. D/C is the opposite operation. I think you can move on from here. Let us know if you can work it out now or if you need further help.


FIR Filter implementation vs Linear convolution implementation

I believe that the code is written in such a way that it gives a correct output only when the coefficients are passed in the reverse fashion. Comparatively, the code will be simpler if we pass the coefficients in reverse fashion.


Your Y axis likely should be scaled to frequency (Hz), not FFT bin index [0..127].  The frequency scale would be something like the bin_index * sample_rate / FFT_length .


Design Hilbert filter at matlab (pure 90 phase and magnitude 0db)

Hilbert Transformers are non-causal, i.e. they need to be delayed to be implementable. So you get the 90 degree phase-shift plus a bulk delay of 20 samples (half the filter length).

You see the 90 degree phase shift if you delay the original signal by 20 samples as well.

EDIT for Part 2:

Your lower bandpass cutoff is too high. It currently sits at 50 Hz. It needs to go down to 40 Hz or potentially a little lower. You may also have to crank up the number of points to get less amplitude ripple in the pass band. Something like

b = firpm(60,[38 950]*2/2000,[1 1],'h');

gets you about 0.1 dB ripple. Cranking the number of tabs up to 96 will reduce the pass band ripple to 0.01dB.


IIR filter with circular buffer - Arduino

First of all :

have you initialized all the previous data to 0 ? The x and y buffer should be initialized to 0. Otherwise, you could get some unexplained transients at the beginning. Maybe this is the offset that you see at the beginning. Maybe the compiler will set the initial data in the data at 0, but you should not rely on the compiler, it should be explicitely reset to 0.

Secondly : I think you have mixed up the c and d coefficients. I think the d coefficients should be used with the output (y) and the c coefficients for the input (x). I say that because d[0] = 1 which makes me think that this the a0 coefficient.  Usually the a0 coefficient is 1 in an IIR. It can be different than 1, but usually we try to normalize the coefficients by setting a0 to 1. https://en.wikipedia.org/wiki/Infinite_impulse_response

Third : Have you tried the coefficients in Octave/Matlab/Python ? Use freqz or Bode in Octave/Matlab. Do you get the expected frequency response?

Finally : This is not related to the matter at hand, but I notice that you use the modulo function. The modulo function is slow unless you're using powers of 2. You should consider having a circular buffer with a size of 4 instead of 10. You don't need more since you only go 4 samples in the past.


Fine-grained evaluation of speech recognition results

You have misconception that errors are caused by either language model or by acoustic model. There is no clear distinction because models are combined together in during search for best hypothesis. Interaction is complex and another major player is beam search pruning.

Most errors are originated from the acoustic model when it fails to give proper estimate for a certain frame, but the problem is that they propagate over next few words. For example, if you have an unknown word, it will not be recognized properly but the next several words will be wrong too because search state will not match the actual content. And phone-by-phone comparison will not be really meaningful.

So it's a bit more complex.


On a modern Pentium computer, it is possible to do many algorithms in real-time at 44.1 kHz. I suspect that what you want to do is possible, but what isn't possible is to tell you how to do it without a better explanation of what you are trying to do, which should include a better definition of what you mean by "peaks at amplitude".


I'm not sure if you ask about symmetry of the impulse response or symmetry of the frequency response. All filters with a real impulse response have conjugate symmetric frequency responses and symmetrical magnitude frequency responses. A symmetrical magnitude frequency response would affect the magnitudes of both frequencies $\omega$ and $-\omega$ identically so it could not possibly remove or attenuate a negative frequency while keeping the positive frequency intact. So the impulse response has to be complex.

The discrete impulse response $b_k$ of an ideal do-nothing filter is:

$$b_k = \frac{\int_{-\pi}^{\pi}e^{i\omega k}d\omega}{2\pi} = \left\{\begin{array}{ll}1 &\text{ if }k = 0\\\frac{\sin(\pi k)}{\pi k}&\text{ if } k \ne 0\end{array}\right.\\ = [\dots, 0, 0, 0, 1, 0, 0, 0, \dots]$$

The middle value $1$ is at $b_0$. Changing where the integral starts gives the ideal negative frequency removal filter:

$$b_k = \frac{\int_{0}^{\pi}e^{i\omega k}d\omega}{2\pi} = \left\{\begin{array}{ll}1/2&\text{ if } k = 0\\\frac{\sin(\pi k)}{2 \pi k} + i\left(\frac{1 - \cos(\pi k)}{2\pi k}\right)&\text{ if }k \ne 0\end{array}\right.= [\dots, 0, -\frac{i}{5\pi}, 0, -\frac{i}{3\pi}, 0, -\frac{i}{\pi}, \frac{1}{2}, \frac{i}{\pi}, 0, \frac{i}{3\pi}, 0, \frac{i}{5\pi}, 0, \dots]$$

The middle value $1/2$ is at $b_0$. The real part of the impulse response is zero at other values of $k$. The imaginary part has odd symmetry with respect to $\operatorname{Im}(b_0) = 0$. It is plain to see from the integral that the frequency response is zero for frequency $\omega < 0$ and constant ($1$) for $\omega > 0$.


The best values for these parameters depend on your application. You should probably try different values on a validation set to find out which ones work best for you.


You typically do a spectral analysis on a time series that has a single value per time instance. For example: what's the distance a car has driven after 1 hour. The answer can be 10 km or 50 km but it can't be both at the same time.

We first we need to understand WHY is there more then one value at the same time

Sometimes data can have multiple dimensions. In our example above you can express the position of the car in x, y, & z co-ordinate in which cases you have 3 values at the same time. In this case you would can analyze each dimension individually or you can derive a single quantity that you are actually interested in (like the total distance).

It could also be statistic, you have multiple estimates for the same data point. You could reduce this by taking the mean or doing a Monte-Carlo type of analysis.

Can you elaborate on the type of data and the reason for the multiplicity?


I've heard this referred to as demeaning a signal, or applying a DC notch (or, in some more hardware-related domains, a DC block). More generally, it might be called just highpass filtering the signal to remove the DC component.


There is an implant in me, is there a way to identify the frequency it is receiving signals on?

It is very hard to detect a device that is passively receiving a signal, and even more so when it is designed to operate undetected. Some random comments and ideas:

In the UK one has to pay taxes on TV ownership. I read somewhere that they employed a technique to detect unlicensed TVs, which are otherwise passive receivers. I don't know how the technique works, but you may want to research it.

Scanning the spectrum around you is useless. There are hundreds of different signals around us all the time; it will be impossible to detect if any one of them is intended for your implant.

Any receiver must use metal components, so you could use a metal detector on yourself. Get a very small receiver (such as one of the tiny WiFi dongles) and verify that your metal detector can detect it on the other side of your head; it should then be able to detect any device inside your body.


Lower-rate means that the rate of information bits you get through is lower, so $(128,32)$ is lower-rate than $(128,64)$.

  Also, would using a lower-rate ldpc code combat well against deletions in the channel or just misread bits?

Classically, codes like these do not help against deletion; they are designed to find flipped bits; which is much more usual in digital communication channels, where the receiver has a clock that keeps running. If you had to model a deletion channel in this framework, you'd just go ahead and assume a 50% BER for all bits following (and including) the deleted bit. You could probably derive some (probably Bernoulli) distribution of how many bits are deleted per word, and when they happen, and thus calculate how many additional bit errors your code needs to be able to correct.

There's no functionality in a LDPC to guess which bits were deleted; idealized, mathematically, you're just trying to stuff $\mathbb F^{128-N_{deleted}}$ vectors into a vector space with dimensionality ${128}$. That's a bit strange.

I'm sure there's awesome math extending certain classes of channel codes to do that well, but I'm not aware of those, and my guess is that they are very complex to understand.

That's why it's usually desirable to have channel coding come directly after the physical recovery of symbols, and not after some component tried to make sense of those (aside from soft decoding gains etc).


How to transmit voice signal and receive on radio?

The first thing you have to know is that transmitting on the FM band is something heavily regulated. Your emitter must not exceed a specific power or range. Check the legislation of the country in which the device operates.

For this you will need an FM transmitter. You can buy them pre-made, or build them with very few parts. It might even be possible to get it to work without additional hardware and some creative hacking (some smartphones or tablets use the BCM4329 chip for example, which is capable of FM transmission).

  And what happen if i ll broadcast on same frequency on which already data is broadcasting ?

The strongest signal wins.


How to choose an interpolation method for seismic ground motion signals

If there are differences then some likely sources might be  that the window on the Sinc wasn't wide enough, that the data wasn't really sufficiently bandlimited, or too much quantization noise in the samples.

The double integration of anything with system offsets or quantization noise is unlikely to be zero. (See random walk) One solution would be to compute the double integral of the first pass interpolated result and then bias the interpolator (noise filtered if needed) until the resulting double integral goes to zero.


What's a good classification algorithm for a large amount of classes?

I would suggest using Deep Learning and especially convolutional neural networks. This might be a generic approach to large scale classification tasks. One such example is identifying faces. As an example, 10.000 classes have been tackled by Sun et. al. in the following : http://ieeexplore.ieee.org/document/6909640/

And iterative learning is of course supported.


Does the Wiener–Khinchin-Einstein theorem hold for non-Gaussian processes? If so are there any assumptions?

first sentence of Wikipedia article:

  In applied mathematics, the Wiener–Khinchin theorem, also known as the Wiener–Khintchine theorem and sometimes as the Wiener–Khinchin–Einstein theorem or the Khinchin–Kolmogorov theorem, states that the autocorrelation function of a wide-sense-stationary random process has a spectral decomposition given by the power spectrum of that process.

(emphasis by me)

So this has nothing to do with Gaussian or not, this has to do with stationarity of your process – a process that is Gaussian for every single point in time isn't necessarily wide-sense stationary.

For example, if you have a random signal with an amplitude that is normal distributed for every instant, but whose variance decreases exponentially with time, you can't apply the Wiener-Khinchin theorem.

If you, on the other hand, have a process that actually only takes discrete values $-1$ and $1$, but does so with constant probability, then you can apply the Wiener-Khinchin theorem.

So:

  Does this relationship hold for non-Gaussian processes?

Yes.

  And if so are there any assumptions or restrictions in place regarding its use?

The actual restriction is wide-sense stationarity, as you've noticed:

  For if it holds for non-Gaussian processes, is it suffient for the process Weak or wide-sense stationarity, i.e. only the first and second order moment need to be stationary and higher order moments can be non-stationary, or are there more strict conditions?

Um, that's not exactly the definition of WSS that I have in mind:

Weak sense stationary means that the autocovariance only depends on the time difference, not on the absolute times. This also means that the autocovariance of zero shift, ie. $E[x(t)x(t)]-\mu(t)\mu(t)$ is constant, which can, by considering that these are two expressions for energy, only happen if the first moment $\mu$ is constant.


Can LMS filter converge to sparse solution?

...can we predict $h$ using LMS algorithm under all conditions?

LMS converges under most "reasonable" conditions.

The weights of the filter are adapted by a "signal" that is derived by the mean square error of the difference between two signals, which is convex. That is, it has a single minimum (when it does) and it will trickle towards it. If you notice the formula by which each weight is adapted, it does not include a cumulative term. It is just adapting the weights by the magnitude of the error. Therefore, it is supposed to "track" the situation smoothly until it eliminates the difference that is driving it. But if you now crank the $\mu$ up, it will be overshooting that single "minimum" in either side and it will go into "oscillations". Not in the "system is unstable and it is now whistling" sense but in the way that the logistic map goes into "oscillations"

  I was under impression that both s and h needs to be stationary to estimate h properly.

They would have to be stationary for the weight adjustment process to "converge". Otherwise, it will keep "hunting".

  In case, h is sparse (few non zero only) LMS cannot converge to sparse h. Is it true or false?

It is false...most of the times.

Given a reasonable choice for the learning parameter $\mu$ and enough time, LMS will converge to something. But, there are two "problems" that will affect your $h$ estimate: The fact that there is noise and the weights will be trying to minimise even that and the fact that you estimate $h$ from a small window in time. Therefore, the signal(s) better have exhausted the full repertoir of their variance within that window, otherwise, LMS will be "hunting" (smoothly, or wildely) around the minimum again.

Here is an example: $h_{channel}$ is innocent, maybe $h_{channel} = [1,0,0,0,0,0,0,0,0]$. You start with some random $h_{predicted}$ and feed through a sinusoid. The filter will now start adapting. Through the first few iterations, you get totally random "guesses". But there will be one guess, purely because of chance, that a coefficient will land higher than the rest. LMS will pick this up immediately and start improving on it. In the end, it will converge to something that looks like $h_{channel}$, if only the noise perturbations were not there. So, your $h_{predicted}$ is not going to look exactly like $h_{channel}$ with a few strong returns to show reflections and zeros in between, it will still have some values interdispersed with the strong reflections.

Try to think "vertically" into the coefficients of $h_{predicted}$. Because the LMS is like an Automatic Gain Control "Bank". LMS doesn't optimise the "filter" as a whole (i.e. looking at it "horizontally", as one $h$). It uses the error signal to tune each filter coefficient separately and through convolution it applies the effect of each one on the signal. Try to work out LMS, by hand, with a two coefficient $h$ to see what I am trying to say with this.

Hope this helps.

EDIT:

I am adding this small "post script" here, given the additional information that was provided in the comments.

It seems that what is sought after is the impulse response of a "corridor" over the medium of Radio Frequencies.

Depending on the room size, its construction and frequency of operation, the system will not even have the chance to "resolve" between returns.

A nanosecond of propagation at the speed of light is about 30cm of "distance". The distinct pulses that are expected here, with large intervals of "silence" between them, will pile up into "averages", creating a "complex" $h$. That $h$, by the way, is not going to be stationary, depending on the building construction, because you cannot control what the RF is bouncing off of OUTSIDE the limits of the corridor.

You might be better off with a slow chirp to obtain the frequency response of the channel (assuming that this is what you are after).

Hope this helps


Reasoning Behind Additional Phase Shift in DPSK

A QPSK constellation could in theory have any rotation. For example, $$\mathcal{S}=\lbrace 1, j, -1, -j \rbrace $$ is a QPSK constellation, and $e^{j\phi}\cdot\mathcal{S}$ is also a QPSK constellation. Matlab's phaserot parameter in the dpskmod command is $\phi$ in the equation above.

What is to be gained depends on the particular application I guess. In particular, note that $\phi=\pi/4$ is nice because then the ML decoding rule is equivalent to finding what quadrant the received point is in.


Seeking compressive sensing imaging demo in MATLAB

You can start with the Compressive Imaging Code code by J. Romberg, illustrating the paper  "Imaging via Compressive Sampling".

Another great source of information and codes is on Nuit Blanche.


What is the cause of slow roll off of windowed filters

A (non-full-width-rectangular) window attenuates data at its edges, thus narrowing the width (decreasing the 1st moment) of the bulk hump of the input.  Decreasing width in the time domain increases uncertainty in the frequency domain, which is represented by a filter response which has a shallower transition band stretching farther from the center frequency.

One way to roughly compensate is to start with a wider data window (longer FIR), which when "narrowed" by a non-rectangular window function will constitute roughly the same bulk "width" as before the widening.


Hint:

So you have proved that
$$\exp\left({j\frac \pi4 n}\right)+ \exp\left({j\frac{7\pi} 4 n}\right) =\exp\left({-j\frac \pi4 n}\right)+ \exp\left({j\frac \pi4 n}\right)$$

The second derivation does results in the same answer as well. It is not prohibited but one would want to simplify things here. It does not help factorizing with a factor $\exp({j\pi n})$ because $\exp({j\pi n}) = \left(-1\right)^n$ and not $-1$. But one could factor by $\exp({j2\pi n})$ for instance since $\exp({j2\pi n}) = 1$. With this in mind:

$$\exp\left({j\frac\pi4n}\right)+ \exp\left({j\frac{7\pi}{4}n}\right)={\exp\left({j2\pi n}\right)}\cdot\left[\underbrace{\exp\left(-j\frac{7\pi}{4}n\right)}_{\exp\left(j\frac{\pi}{4}n\right)}+ \exp\left(-j\frac{1\pi}{4}n\right)\right]$$
And the two forms are identical.


For a wide-sense-stationary random process, all the random variables
comprising the process have the same mean $\mu$ and variance $\sigma^2$, and
the variance is the integral of the power spectral density $S(f)$ less the square
of the mean:

$$\sigma^2 = \int_{-\infty}^{\infty} S(f)\,\,\mathrm df - \mu^2.$$

Generally, the mean is $0$ and $\sigma^2$ is just the integral of
$S(f)$, but when $\mu \neq 0$, then $S(f)$ includes
an impulse of $\mu^2\delta(f)$ at $f=0$ which contributes $\mu^2$ to
the integral, and this gets subtracted off by the $-\mu^2$ term in the
above formula. For discrete-time signals, the integrals should
be replaced by the appropriate sums.  But when there is a
(deterministic) signal component and a (random, zero-mean)
noise component, calculation of the signal-to-noise ratio SNR
(which has many definitions including one which says it is
the ratio of the signal energy to the noise variance times
the signal duration), requires separating out
the two numbers needed from knowledge of their sum (which is
what the periodogram will give you). This is a lot trickier
and the answer will typically depends on details of the signal that you
have not shared with us.


How do I add AWGN to an I and Q representation of a signal?

Yes, you can add AWGN of variance $\sigma^2$ separately to each of the two terms, because the sum of two Gaussians is also a Gaussian and their variances add up. This will have the same effect as adding an AWGN of variance $2\sigma^2$ to the original signal. Here's some more explanation if you're interested.

An analytic signal $x(t)=a(t)\sin\left(2\pi f t + \varphi(t)\right)$ can be written in its in-phase and quadrature components as

$$x(t)=I(t)\sin(2\pi ft) + Q(t)\cos(2\pi ft)$$

where $I(t)=a(t)\cos(\varphi(t))$ and $Q(t)=a(t)\sin(\varphi(t))$. If you wish to add AWGN to your original signal as $x(t)+u(t)$, where $u(t)\sim\mathcal{N}(\mu,\sigma^2)$, then you can add AWGN to each of the terms as

$$y_1(t)=\left[I(t)\sin(2\pi ft) + v(t)\right] + \left[Q(t)\cos(2\pi ft) + w(t)\right]$$

where $v(t), w(t)\sim\mathcal{N}(\mu/2,\sigma^2/2)$

Also note that because the in-phase and quadrature terms are additive, the AWGN can also be simply added to either of the two terms in the $IQ$ representation of $x(t)$ above. In otherwords,

$$y_2=I(t)\sin(2\pi ft) + \left[Q(t)\cos(2\pi ft) + u(t)\right]$$
$$y_3=\left[I(t)\sin(2\pi ft) +u(t)\right]+ Q(t)\cos(2\pi ft)$$

are statistically equivalent to $y_1$, although I prefer using $y_1$ because I don't have to keep track of which component has noise added to it.


What scale should a detected SIFT point be extracted at?

If I understand correctly what you are asking --
In general, the feature is found at the same scale as the SIFT detector says, but
in David Lowe's SIFT, the image is pre-smoothed with sigma: 0.5, so, you need to "subtract" this amount of smoothing from the sigma, so the "real" scale could be:

sqrt(sigma^2 - 0.5^2)

where sigma is the scale the feature was found at.
you have an excellent implementation of SIFT detector and descriptor at: www.vlfeat.org


Welch’s method is a scheme to reduce the fluctuations (noise) in spectral power estimations when you have only a single time sequence with which to work. Because you have multiple time sequences, I’m not convinced that Welch’s method is what you want. If you have time sequences $x1(n)$, $x2(n)$, and $x3(n)$ of different lengths, I suggest you zero-pad the shorter sequences so they’re equal in length to the longest of the original $x1(n)$, $x2(n)$, and $x3(n)$ sequences. (It’s my guess that no time-domain Hanning or Hamming windowing is necessary.) Perform three FFTs and compute the three corresponding spectral power sequences $P1(m)$, $P2(m)$, and $P3(m)$. Finally, average the three spectral power sequences, element-by-element, as:

Spec_Pwr_Ave = P1 + P2 + P3;

Then plot your ‘Spec_Pwr_Ave’ power-spectral sequence.


Column operations in LDPC generator matrix

Everything in me screams that's wrong, but it turns out if I just do the same thing on the decoding side it works out.

Why? That just means you reorder some bits. On the other end, when you revert that reordering, that's an operation that changes nothing about the code itself.

In fact, your system almost certainly has an interleaver, anyway, so you such reordering operations are absolutely integral to how it works!

  the conundrum that when I go to decode something encoded with this parity check matrix I have to know how the matrix columns were reordered on the transmit side.

What's surprising about that? Your receiver needs to have a decoder matching the encoder.

Don't forget that swapping columns is just a permutation matrix that you apply your data. Done!

  It seems bizarre to me that a designed system would put in a parity check matrix that requires column reordering to work.

Why? You say that, but you don't argue based on anything about that. In fact, many coding-theory proofs are based on random codes...

Again, imagine your reordering operations as permutation matrices that you just use after / before applying the parity check / generator matrix.


Determine when a constant signal changes to another constant value

As @AlexR mentioned, the digital 1st order derivative is the most straight forward and probably the most trivial solution to the problem, depending on what your signal looks like and the rate of changes involved between samples. The derivative can then be processed in a sliding window to find where it transitions and where it settles (derivative of constant will be close to zero). With these techniques you can reliably detect constants ONLY if the derivative of the signal is greater than the noise amplitude before the signal settles to its constant value. If the rate of change of the signal is slow (i.e. the sample to sample transitions are on the order of magnitude of your noise) then it will be difficult to detect the constant signal from this method. I scale the signals below to be on the same order of magnitude and also take the absolute value of the derivative.

close all
clear all
%signal
x = [ones(1,5)*255 253 251 ones(1,6)*240 240 230 220 200 193 187 180 170 164 140 130 130 120 110 80 ones(1,10)*70 80 95 112 125 137 140 170 180];

%repeat signal
y = repmat(x, [1 5]);

%max value 2^8
MAX = 255;

%noise amplitude
Namp = 5;

%add noise to signal
y = y+randi(Namp,1,length(y));

%1st order derivative
c = diff(y,1);

%plots
figure
stem(abs((c)/MAX),'color','b')
hold on
plot(y/max(y))
line([1,length(y)],[0,0],'color','r','linewidth',2)
line([1,length(y)],[Namp/MAX,Namp/MAX],'color','r','linewidth',2)
ww = 4
reg = zeros(1,ww)
for i = 1:length(c)-1
  reg = [reg(2:end) abs(c(i)/MAX)];
  if(reg(1) >Namp/256 & sum(reg<Namp/256)>= 3 )
      line([i-ww+2, i-ww+2],[0 1],'color','g','linewidth',2)
  end
end
title({'Red: Signal', 'Blue: Derivative', 'Green: start of constant signal'})
hold off

The signals generated from the derivative algorithm can be used to make a clean piecewise continuous function as seen below. The delay of which, will be a function of your window size-1 sample.

Updated code:

close all
clear all
%signal
x = [ones(1,20)*255 253 251 ones(1,4)*240 255 255 255 255 ones(1,5)*240 240 230 220 200 193 200 220 ones(1,10)*230 210 190 187 180 170 164 140 130 130 120 110 80 ones(1,10)*70 80 95 112 125 137 140 170 180];
%x = [ones(1,5)*130 ones(1,10)*255 ones(1,10)*20];
%repeat signal
y = repmat(x, [1 5]);

%max value 2^8
MAX = 256;

%noise amplitude
Namp =5;

%add noise to signal
y = y+randi(Namp,1,length(y));

%1st order derivative
c = diff(y,1);
c2 = diff(y,4);

hold off
%plots
figure
stem(abs((c)/MAX),'color','b')
plot(abs(c2)/MAX,'color','r')
hold on
plot(y/max(y))
line([1,length(y)],[0,0],'color','r','linewidth',2)
line([1,length(y)],[Namp/MAX,Namp/MAX],'color','r','linewidth',2)

ww = 4;
reg = zeros(1,ww);
filtered = []
piecewise =y(1);
for i = 1:length(c)-1
  reg = [reg(2:end) abs(c(i)/MAX)];

  %if we detect a transition to 0 derivative
  if(reg(1) >Namp/256 & sum(reg<Namp/256)>= 3 )

      line([i-ww+2, i-ww+2],[0 1],'color','g','linewidth',2)
      piecewise = y(i);
  end
  %if we detect a transition out of 0 derivative
  if( sum(reg<Namp/256)>= 3 & reg(end) >Namp/256 )
      line([i, i],[0 1],'color','y','linewidth',2)
      %piecewise = y(i)
  end

  filtered = [filtered piecewise];

end
title({'Red: Signal', 'Stem Plot: Derivative', 'Green: start of constant signal', 'Yellow: end of constant signal'})
hold off

figure
plot(filtered(3:end))
hold on
plot(y)
title({'Red: Signal', 'Blue: Cleaned Signal'})

If we take $x[n]$ to be the input and L to be the window length. First we take the 1st order digital derivative of $x[n]$ to get $\hat{x}[n]$

$$\hat{x}[n] = x[n+1]-x[n]$$

and then we examine $\hat{x}[n]$ in an L length window, noting that when the oldest sample in the window is large (signal is transitioning) and the rest of the samples in the window have settled close to zero (signal is constant) we are at an area of interest (we have found a constant signal after a transition). The longer the window L, the more samples you are considering after the transition before you trigger your algorithm. That is, if you know the signal will settle to a constant for a long period of samples, then choose a big L. If you want the output to transition quickly though, choose a smaller L. We know the smallest transition is zero but will be on the order of $N_{amp}$, the amplitude of the noise. If we scale everthing to the maximum amplitude the signal can be, $A_{max}$, then the newest L-1 samples are on the order of $\frac{N_{amp}}{A_{max}}$ and the oldest sample $x[n-L]$ is some value above that by $\epsilon$

We need to capture this mathematically to say: when the first $L-1$ samples are less than or equal to the noise amplitude and the sample L before is larger than the noise amplitude we have gone from a transition to a constant and the output should be the sample at this instance. Intuitively it seems that:

$$\frac{1}{(L-1)} \sum_{i=0}^{L-1}\frac{|\hat{x}[n]|}{A_{max}} \leq \frac{N_{amp}}{A_{max}}$$

should capture some essence of part of this, but doesn't quite capture that ALL samples in this $L-1$ window should be less than or equal to the noise amplitude.

it's tricky to capture this if-else statement mathematically. One way to get close to a single equation is to say the areas of interest (green lines in my plot) occur when

$$r[n] = \frac{\hat{x}[n-L]}{\frac{1}{(L-1)} \sum_{i=0}^{L-1}|\hat{x}[n]|} > 1 + \frac{\epsilon}{N_{amp}}$$

since the division above expresses:

$$\frac{\frac{N_{amp}+\epsilon}{A_{max}}}{\frac{N_{amp}}{A_{max}}}$$

This works ok, for example if you say you are only interested in $\epsilon$ greater than 5x your noise amplitude, then you get close, but it misses some transitions. It also seems to behave oddly when the window increases.

That is to say

$$y[n] = \left\{
\begin{array}{ll}
  x[n], \quad & r[n] > 1 + \frac{\epsilon}{N_{amp}} \\
  y[n-1], \quad & otherwise \\
\end{array}
\right. $$

code to reflect this math:

    close all
clear all
%signal
xx = [ones(1,50)*255 253 251 ones(1,50)*240 ones(1,4)*255 ones(1,10)*240 240 230 220 200 193 200 220 ones(1,20)*230 210 190 187 180 170 164 140 130 120 110 80 ones(1,25)*70 ones(1,2)*256 ones(1,25)*70 80 95 112 125 137 140 170 180];
%x = [ones(1,5)*130 ones(1,10)*255 ones(1,10)*20];
%repeat signal
x = repmat(xx, [1 5]);

%max value 2^8
MAX = 256;

%noise amplitude
Namp =5;

%add noise to signal
x = x+randi(Namp,1,length(x));

%1st order derivative
xhat = diff(x,1);

L = 5;
reg = zeros(1,L);
filtered = [];
y = [];
y(1:L) =x(1:L)

for n = L+1:length(xhat)-1
  LminusOneSum = 0;
  for i= 1:L-1
    LminusOneSum = LminusOneSum + abs(xhat(n-i))/MAX;
  end
  r(n) =  abs((xhat(n-L))/MAX)/(LminusOneSum*(1/(L)));

  if(r(n) > (1+6))
      line([n-L+2, n-L+2],[0 1],'color','g','linewidth',2)
      y(n) = x(n);
  else
      y(n) = y(n-1);
  end

end
title({'Red: Signal', 'Stem Plot: Derivative', 'Green: start of constant signal', 'Yellow: end of constant signal'})
hold off

figure
plot(y(L:end))
hold on
plot(x)
%plot(r)
title({'Red: Signal', 'Blue: Cleaned Signal'})


How to correct the phase offset for QPSK I-Q data

Assuming we have timing recovery resolved (meaning our system knows the correct time locations for our symbol decisions), then we can use the decided symbols along with our pre-decision values in a Decision-Directed Phase Detector. This will give us an error value that we can then feed into a phase lock loop to correct for phase error using traditional PLL methods.
Decision Directed Phase Detector
See my figure below for the implementation of a decision directed phase detector. The two vectors shown are the complex sample before making a decision ($V_2$), and the complex sample after the decision, which represents the closest location for each symbol in our locked constellation ($V_1). For QPSK, the normalized decisions would be 0.707 ± j0.707.

This describes a very simple approach to get a phase error term for QPSK (and higher order QAM systems) as an implementation of the phase discriminator. This is in fact a compact representation of the phase detector in the Costas Loop for a QPSK demodulator. The rest of the work for tracking assumes that you have knowledge in Digital PLL implementations to "close the loop" and is not covered here.
To explain the Decision Directed Phase Detector further for those extremely interested, I detailed below the more generic "cross product phase detector". Understanding this will give us full insight into the operation summarized in the figure above.
Cross Product Phase Detector
The cross product phase detector works based on the property that the imaginary portion (cross product) of a complex conjugate multiplication of two vectors is directly proportional to the sine of the phase angle between the vectors.  This has wide utility in digital carrier and phase tracking algorithms.  The easy way to see this is knowing that when we multiply two vectors, the angles add; so if we conjugate one of the vectors, the angles will subtract giving us the difference of the two angles.

The math for the actual complex conjugate multiplication reveals how simple the final  implementation can be (two real multiplies and an add), which is what makes this approach so prevalent:
Consider two vectors $V_1$ and $V_2$
$V1= I_1+jQ_1$
$V2= I_2+jQ_2$
$V_1^*V_2 = (I_1-jQ_1)(I_2+jQ_2) = (I_1I_2+Q1Q2) + j(I_1Q_2-I_2Q_1) = I_P+jQ_P$
The angle of the result is $\phi = atan2(Q_P,I_P)$
However in most cases for this application the processing intensive atan2 is not necessary, as we can make use of certain approximations as follows:
First the imaginary term $Q_P$ is proportional to the sine of the angle:
$Q_P = |V_1||V_2|sin(\phi)$
The proportionality constant just becomes part of our overall loop gain; for a phase detector we just need something that is linearly proportional to phase.
Second for small angles, $sin(\phi) \approx \phi$
A loop would drive the angle (error) to 0, therefore when tracking the angle will be small, so this approximation is valid for assessing the tracking performance of the loop. For acquisition considerations, we need to be aware of the reduced loop gain due to the changing slope of a sine function (unless calculating an atan which has the benefit of being linear throughout its range).
So bottom line for a cross product phase detector we take advantage of the following relationship between the I and Q values of two complex samples:
Given two complex samples [$I_1,Q_1$] and [$I_2,Q_2$], the cross product of the two samples is proportional to the phase between them:
$I_1Q_2-I_2Q_1 \propto \phi$
Thus we have provided a very simple and elegant way to provide phase discrimination in an all-digital PLL.
Note: The following result is rotation invariant, meaning if the constellation rotates past 45 degree error, it will lock in another position. A known training sequence (in the header) would be required to establish the correct orientation for your packet of data.


How to find the best Q format representation?

You need to calculate the transfer function from your input to each individual state variable.

This depends A LOT on how you implement your filter: I strongly recommend splitting it in second order sections and using either Direct Form I or transposed Form II for each section. Section order and pole/zero pairing can also make a big differences. Unless you have a REALLY good reason not to, follow this recipe: https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.zpk2sos.html

Once you have your filter topology sorted out you need to calculate all the transfer functions from the input to the state variables and look at the  frequency responses. At this point you can also play around with the individual section gains to make this as smooth as possible and manage outliers and peak frequencies.

Then you can pick the Q based on the maximum gain over all frequencies in each transfer function. Add another 3 dB of headroom for transients (sine waves are NOT the worst case signal).


Normalization of a signal with respect to another signal

I have had the same problem. My signals are coming from the MRI machine. In my case, the scaling is one of my options; in a way, I divided the signal A and B to the maximum of signal A, but I could not use z-score approach because the information has been embedded in the amplitude. If your signal is a non-amplitude modulated like a speech signal, you can use a z-score.


FFT vs DFT Run Time Comparison (Complexity Analysis) in MATLAB

Abhinav Jain, Welcome to DSP Community.

I build for you a proper testing of the run time comparison.

Few tips about timing in MATLAB:

Never time in a script. Always call a function to do the heavy lifting. When you run something from script it runs in the global scope which means MATLAB can't optimize it as well as it could if it was in a function.
When you time a function, make few iterations of the measurement. I like taking the median of this iterations. Others like the minimum. Usually the mean isn't good as it is sensitive to outliers.
MATLAB is a JIT based Script language. It means the first run of anything takes more time. Hence (2) is crucial.
Pay attetion that the result of the FFT / DFT is complex. Hence when you allocate memory for a complex array you should use - vArrayName = complex(zeros(arrayLength, 1));.
When you need the imaginary number $ i = \sqrt{-1} $ in MATLAB you should use 1i or 1j (No need for multiplication). Of course you can use 5i as well. Yet no need to multiply as it assists the JIT engine to understand you mean the imaginary number and not a variable.

Here are the results:

As you can see, indeed the DFT Implementation has running time which behaves as O(n^2).

The code goes like:

%% Simulation Parameters

vNumSamples     = 2:2:1024;
numIterations   = 6;

%% Generate Data

mDftTime = zeros(numIterations, length(vNumSamples));
mFftTime = zeros(numIterations, length(vNumSamples));

for jj = 1:length(vNumSamples)
    numSamples = vNumSamples(jj);
    vX = randn(numSamples, 1);

    for ii = 1:numIterations
        hDftTimer           = tic();
        vXDft               = ApplyDft(vX, numSamples);
        mDftTime(ii, jj)    = toc(hDftTimer);

        hFftTimer           = tic();
        vXFft               = fft(vX);
        mFftTime(ii, jj)    = toc(hFftTimer);
    end

end

%% Run Time Analysis

vDftMedian = median(mDftTime).';
vFftMedian = median(mFftTime).';

vDftMean = mean(mDftTime).';
vFftMean = mean(mFftTime).';

vDftMax = max(mDftTime).';
vFftMax = max(mFftTime).';

vDftMin = min(mDftTime).';
vFftMin = min(mFftTime).';

%% Display Results

figureIdx = figureIdx + 1;

hFigure = figure('Position', figPosLarge);
hAxes   = subplot(4, 1, 1);
% set(hAxes, 'NextPlot', 'add');
hLineSeries = plot(vNumSamples, [vDftMedian, vFftMedian]);
set(hLineSeries, 'LineWidth', lineWidthNormal);
% set(hLineSeries(2), 'LineStyle', ':');
set(get(hAxes, 'Title'), 'String', {['DFT vs. FFT Run Time - Median']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', {['Input Size']}, ...
    'FontSize', fontSizeAxis);
set(get(hAxes, 'YLabel'), 'String', {['Run Time [Sec]']}, ...
    'FontSize', fontSizeAxis);
hLegend = ClickableLegend({['DFT'], ['FFT']});

hAxes   = subplot(4, 1, 2);
% set(hAxes, 'NextPlot', 'add');
hLineSeries = plot(vNumSamples, [vDftMean, vFftMean]);
set(hLineSeries, 'LineWidth', lineWidthNormal);
% set(hLineSeries(2), 'LineStyle', ':');
set(get(hAxes, 'Title'), 'String', {['DFT vs. FFT Run Time - Mean']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', {['Input Size']}, ...
    'FontSize', fontSizeAxis);
set(get(hAxes, 'YLabel'), 'String', {['Run Time [Sec]']}, ...
    'FontSize', fontSizeAxis);
hLegend = ClickableLegend({['DFT'], ['FFT']});

hAxes   = subplot(4, 1, 3);
% set(hAxes, 'NextPlot', 'add');
hLineSeries = plot(vNumSamples, [vDftMax, vFftMax]);
set(hLineSeries, 'LineWidth', lineWidthNormal);
% set(hLineSeries(2), 'LineStyle', ':');
set(get(hAxes, 'Title'), 'String', {['DFT vs. FFT Run Time - Max']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', {['Input Size']}, ...
    'FontSize', fontSizeAxis);
set(get(hAxes, 'YLabel'), 'String', {['Run Time [Sec]']}, ...
    'FontSize', fontSizeAxis);
hLegend = ClickableLegend({['DFT'], ['FFT']});

hAxes   = subplot(4, 1, 4);
% set(hAxes, 'NextPlot', 'add');
hLineSeries = plot(vNumSamples, [vDftMin, vFftMin]);
set(hLineSeries, 'LineWidth', lineWidthNormal);
% set(hLineSeries(2), 'LineStyle', ':');
set(get(hAxes, 'Title'), 'String', {['DFT vs. FFT Run Time - Min']}, ...
    'FontSize', fontSizeTitle);
set(get(hAxes, 'XLabel'), 'String', {['Input Size']}, ...
    'FontSize', fontSizeAxis);
set(get(hAxes, 'YLabel'), 'String', {['Run Time [Sec]']}, ...
    'FontSize', fontSizeAxis);
hLegend = ClickableLegend({['DFT'], ['FFT']});

if(generateFigures == ON)
    saveas(hFigure,['Figure', num2str(figureIdx, figureCounterSpec), '.png']);
end

As you can see in the code above I run numIterations = 6 each algorithm for each number of samples.
Then, on the time array mDftTime / mFftTime I can analyze the Median / Mean / Min / Max run time.
This is the proper way to analyze function run time in MATLAB.

The full code is available on my StackExchange Signal Processing Q51516 GitHub Repository.

Remark
You're doing your first steps in DSP world and big part of being modern engineer is about programming code.
You should invest time in learning to properly write code and style in a way it is easy to under stand and manage for you and your partners.
Try learning from other people examples.


If you consider the whole set of potential wavelet transforms, then you have a lot of flexibility.

For instance, should you use 1D continuous complex wavelet transforms, by analyzing the modulus and the phase of the scalogram, and provided you use well-chosen wavelets (potentially different for the analysis and the synthesis), and a proper discretization, you can:

find discontinuities and irregularities of a signal and its derivatives
find break point location by wavelet ridge extrapolation
denoise
perform matched filtering based on templates (with complex continuous or discrete dual-tree wavelet frames)
analyse (multi-)fractalty
analyse frequencies (with Gabor wavelets for instance)

Due to the redundancy, and the quantity of available wavelets (not the same is best for different purposes), they could  appear a little less efficient for the analysis of pure stationary and harmonics signals, for which Fourier is better suited.

The main drawbacks are:

for fine analysis, it becomes computationaly intensive
its discretization, the discrete wavelet transform (comp. efficient), is less efficient and natural
it  take some energy to invest in wavelets to become able to choose the proper ones for a specific purpose, and to implement it correcly.


Is there an OBJECT BASED Image Segmentation method, with pure non-segmented training images?

This is not a complete answer (I would comment directly your question but not allowed to due to rep), but you could look at "Learning and detecting shape models v1.3 (September 2009)" (http://groups.inf.ed.ac.uk/calvin/software.html). The paper's name is: "From images to shape models for object detection" by Ferrari, Jurie & Schmid

The idea is to learn a contour based model of object which is then possible to deform. Actually, you have to segment a little bit by providing simply bounding boxes (x-y coordinates of the top-left and bottom-right corner) of the object in the training images. You may also use training images that are already cropped around the object which you want to learn the model.

I tried their Matlab implementation and creating the bounding boxes for my training images was about as fast as manually cropping them in Gimp or else...

EDIT: BTW, if you wonder how to create groundtruth files with this program, I used a text editor such as "Notepad++" to create a text file. In the file, I wrote, separated by spaces, the coordinates of the bounding box, e.g.:" 10 33 572 458 " (no quote) then saved the file with an appropriate name (check how they do for their sample files) with the dummy extension ".groundtruth", their Matlab implementation will recognize the files.


Inverse $\mathcal Z$-transform of system with an 8th order pole

Applying the geometric series directly

$$
H(z) = \frac{1}{1-\alpha z^{-8}} = \sum_{n=0}^{\infty}(\alpha z^{-8})^n = \sum_{n=0}^{\infty}\alpha^n z^{-8n}
$$
If $h_n = \mathcal Z^{-1}\{H(z)\}$, then
$$
h_n = \begin{cases}
a^{n/8} \; , \;  \text{n is multiple of 8} \\
0 \; , \; \text{otherwise}
\end{cases}
$$


Your Fourier space image is purely real, so the inverse Fourier is going to be radially symmetric about the origin. Because the origin is at (0,0), it takes a shift to put the inverse image in the center of your matrix.


Translating Matlab scrambler function to C

I think the mistake would be this line

current[0] = xor ^ bin[bin_ndx] ^ 1;

Remove the xor with 1 at the end.

current[0] = xor ^ bin[bin_ndx];

The MATLAB figure shows 1 because that is the multiplier (scale factor) and may not mean xor with 1.


Detection milliseconds signal using a DVB-T dongle as SDR in C (need idea)

Your C program can link to librtlsdr (source available on github), set the frequency and sample rate, then filter the incoming IQ sample stream from the RTL USB dongle (demodulate if necessary), and (one possibile algorithm out of many) implement an envelope matched filter.


Direct Digital Synthesis - Super Sampling - Unknown Harmonics

I think that with a 250 MHz sampling clock, you should be able to generate a signal with a 125 MHz max frequency.
Also, what type of DDS do you use? Is it based on a table of sine/cosine or is it a CORDIC algorithm? It has an impact on precision.
Is the signal path from DDS to DAC ok? Do they work at the same sampling frequency?
If not, a resampling module should be used or correctly configured.


What is the formula for the frequency spectrum?

You are likely looking for the "Power Spectral Density" which is determined by taking the complex conjugate multiplication of the Fourier Transform. The complex conjugate process will give the expected real result and represents a power quantity (which has no phase). For the DFT this is often presented in normalized fashion such as dBFS (dB relative to full scale) or dBc (dB relative to the total power in the signal).  If the time domain signal is in volts, then scale the DFT by $1/N$ and the complex conjugate multiplication will represent the power in each bin in terms of Watts assuming a 1 ohm load (or adjust further for true power if the resistance is known). Since it is now a power quantity to convert to dB use $10Log_{10}(F(\omega)F^*(\omega))$. Where $F(\omega)$ represents the DFT or DTFT result with proper scaling, depending on which units are used.

Further as a "density" the result will be power per unit of frequency (such as dBc/Hz). The DFT and DTFT without any further windowing has an equivalent noise bandwidth that is the inverse of the time duration $T$ of the time domain signal, so therefore is equal to $1/T$ (in whatever units are used for $T$). Windowing complicates this further, but see this post for more details on that:

How to calculate resolution of DFT with Hamming/Hann window?


Why the energy is a good feature extraction for detect disturbances in signal processing?

In short, energy is not a good feature per se. Disturbance enhancement sometimes requires preprocessing, and might not be resistant to outliers.
The notion of disturbance in your  question could be made more  precise.

However, since the energy squares amplitudes of signal values, and is generally integrated over several samples, it has the ability to emphasize disturbances whose amplitude dynamics vary from the "normal" ones.

Aside, within the scope of automated detection, energies are quantities that are relatively easy to minimize: least-squares are common techniques in signal processing.

Yet I would not consider energy as a good feature in general. Depending on the need, other nonlinear operators can be used, as can be seen in the family of Teager Kaiser operators, see for instance: Signal processing using the Teager Energy Operator and other nonlinear operators.


Why not use synthetic datasets for training machine vision (deep learning) algorithms?

Well, when you're training something against a dataset, you're training it to work on that kind of data.
Since we have no idea of what the neural network "uses" from the input images to fulfills its purpose, these might be things where our computer graphics aren't quite like real-world images.
We're "trained" with real-world imaging, perceived through our eyes with our brains. The neural net has none of that. It might simply do the segmenting, and that's actually pretty likely, on the specific kind of anti-aliasing done by the 3D renderer. Great, now you have a neural net that can segment generated images very well, but can't deal with real-world photography. Not quite sure what you've gained there!
But: sometimes, you actually do that, if you can describe your real-world data well enough – but that's usually things like RF signals, abstracts datasets etc., and often you'd intentionally add noise, so that the noise equivocation is more prominent than the differences between synthesized data and the real world.
You'd always have to verify, either way.


Separate two sinusoids with very close frequencies (5 Hz difference) in a 44.1 Khz audio signal, with good time-domain resolution

I believe there are a couple of reasons why your answere is not coming out as expected.  What you have in your problem is the spectrum of the sum of three signals

$G(\omega) = S_1(\omega) + S_2(\omega) + N(\omega)$

where $G(\omega)$ is the spectrum you showed, $S_1(\omega)$, and $S_2(\omega)$ are the spectrum of the signals of interest and $N(\omega)$ is the noise spectrum.  You are attempting to recover signals $S_1(\omega)$, and $S_2(\omega)$ with high fidelity, which may not be possible in this case.  For simplicity, I won't worry about noise here.

Your solution assumes a couple of things.  First it assumes a very small bandwidth for each signal.  You mention that the recovered signals have a slow fade in, which would indicate a band-limited signal.  As an example, a infinite length sinusoid in the time domain would be delta function in frequency.

Also, the spectrums $S_1(\omega)$, and $S_2(\omega)$ overlap, which is problematic.  Think of it as a system of linear equation.  In your case, you have two unknowns, $S_1(\omega)$, and $S_2(\omega)$, but only one equation.  This leads to an underdetermined system that has 0 or infinite solutions.

To overcome this situation you may attempt to assume additional information about the signals $S_1(\omega)$, and $S_2(\omega)$.  Your previous assumptions of limited bandwidth and non overlapping spectrums lead to your previous answer.  In looking toward a better solution, I would suggest new assumptions.  For example, you could attempt to model your signals $s_1(n)$, and $s_2(n)$ as windowed sinusoids and attempt to estimate window function from the data.  Additionally, you may attempt to use a windowing function on the time domain signal to limit the side lobes of the frequency spectrum to improve isolation as well.


Part of your misunderstanding comes from the fact that there are many ways in which the radar signal processing chain is implemented. Depending on the type of radar, targets of interest, hardware, etc., some methods are more appropriate than others. We will consider pulsed-Doppler radar here.

In the chain you describe:

In modern pulse-Doppler systems using pulse compression, the act of performing pulse compression is virtually synonymous with matched-filtering. You can combine these two terms to mean the same thing and simplify your thinking. You are requiring that range detection happen before processing Doppler, but this not be the case. You can very well detect a target using its Doppler information only while being denied range. It is mainly during acquisition of the target that being denied either of the measurement dimensions is truly detrimental. In the tracking phase, you can "track through" and continue to update target velocity information using Doppler only to then later be able to receive good range measurements.

In addition, the system you describe is a more classical approach of using a Doppler filter bank to separate clutter and to determine the target's radial velocity. Many radar systems require that a signal make it past the Doppler filter bank to even consider performing a measurement. Hence why the block diagram you show has the "detection" block coming after the "Doppler filtering" block.

What the block diagram suggests:

The block diagram can represent two main schemes:

Single echoes from targets that make it past the filter bank are considered to perform a measurement. If a target appears at the output of one of these filters (i.e. adequate SNR and mitigated clutter), you can declare a detection where range and Doppler are measured simultaneously. Multiple pulses may be integrated to improve SNR.
The radar system performs Doppler processing by forming a range-Doppler map. This map can be seen as a 2D array formed by collecting multiple target returns after matched filtering for multiple pulse-repetition intervals (PRIs). Once a certain number of matched-filtered returns are collected (usually a power of 2), a DFT is performed in what is called the "slow-time" dimension. In other words, a DFT is performed across every range gate (a range gate is equivalent to each sample of the matched-filter output). After this step, you now have a 2D map where one axis is the delay (via the matched-filter) and the other axis is the Doppler frequency. Detection now finds the cell(s) according to a certain threshold. A cell is "detected" if it passes a certain threshold, and you have performed range and Doppler measurements simultaneously by virtue of simply picking the cell. Hence, you use the range-Doppler map to discriminate using Doppler and perform detection at the same time! You can go further and use both measurement dimensions to discriminate possible detections. One application being to keep track of an airborne radar's altitude return.

These are just the basics! You can implement all sorts of filtering and signal processing algorithms on the range-Doppler map alone to get even more information about your environment and your targets of interest.


If you read the link which describes the second mapping
$$ w = \frac{z+1}{z-1}$$

  Thus inside of the unit circle in z-plane maps into the left half of w-plane and outside of the unit circle in z-plane maps into the right half of w-plane. Although w-plane seems to be similar to s-plane, quantitatively it is not same

It states that $w$ is not the $s$ plane.

Therefore, if your aim is in studying the bilinear transformation for digital filter design then you can safely use the first, well known, bilinear transformation relation, and its inverse:

$$\boxed{s \leftarrow \frac{2}{T}\frac{z-1}{z+1}}$$

Otherwise, the paper you linked is dealing with stability analysis of discrete-time systems using bilinear transform of its own selection...


Visualizing the "frequency shift" theorem in FFT

In the line

Inverse_GShifted=real(ifft(GShifted));

I have changed to

Inverse_GShifted=real(fftshift(ifft(GShifted)));

This is because the ifft output is from $0 \le n \le N-1$, but your actual time output is symmetric about $n=0$. So what you have to do is to center it around $n=0$ by doing fftshift() command. It will give the output of $-N/2 \le n \le N/2-1$ whose plot is below

Regarding the zoomed output, the IFFT of gaussian function $G(f)$ is again Gaussian $g(t)$. It is getting modulated or point-by-point multiplied by $e^{j\omega_0 n}$ where $\omega_0 = 1400/4500\pi$ (considering $-\pi \le \omega_0 \le \pi$), which is the shift in frequency that you applied. ie, $g[n].*e^{j\frac{14}{45}\pi n}$. The period of this waveform turns out to be $\approx 6.43$ so you will see the peaks repeated after every 6 or 7 samples which is correct as per the zoomed plot. $14/45$ being irrational may be the reason why period is not an integer.


How can I generate offset QPSK symbols from QPSK symbols

The idea of OQPSK is very simple:

Consider your QPSK as two BPSK, one orthogonal to the other. It helps remembering that you're free to rotate the constellation, so that the constellation points are on the I and Q axis. Instead of choosing one point of a four-points constellation based on two bits, you now choose two separate constellations points from two two-point constellations.
pulse shape I and Q independently with a pulse shaper that has more than 2 Samples per Symbol. Typically, use 5 or more.
offset (in time) either component so that it lags half a symbol time behind the other

If you already have QPSK-modulated baseband, you can still do step 3.


Finding the Amplitude spectrum of a FM signal

The difficulty with respect to this question as well as a related one is that the authors of the text that the OP is using have chosen to adopt a definition of an FM signal that is not
the one used by most other writers.  This unfortunate choice muddies
the study of what is already a difficult-to-understand modulation method
(as compared to amplitude modulation and even phase modulation). Most people
use take an FM signal to be of the form given in
Wikipedia where,
for a modulating signal $x(t)$ (assumed to be of magnitude at most $1$), the FM signal is
of the form
$$A \cos\left(2\pi f_c t + 2\pi f_\Delta \int_0^t x(\tau)\,\mathrm d\tau \right)$$
where $f_c$ is the carrier frequency and $f_\Delta$ is the maximum frequency deviation from the carrier.
The instantaneous frequency is the derivative of the phase (argument of the $\cos$) with respect to time $t$, and is thus
$$\begin{align*}
\omega_{\scriptstyle{\text{inst}}}
&= \frac{\mathrm d}{\mathrm dt}\left[2\pi f_c t
+ 2\pi f_\Delta \int_0^t x(\tau)\,\mathrm d\tau \right]
= 2\pi f_c + 2\pi f_\Delta x(t) ~ \text{rad/sec},\\
f_{\scriptstyle{\text{inst}}}
&= \frac{1}{2\pi}\frac{\mathrm d}{\mathrm dt}\left[2\pi f_c t
+ 2\pi f_\Delta \int_0^t x(\tau)\,\mathrm d\tau \right]
= f_c + f_\Delta x(t) ~ \text{Hz}.
\end{align*}$$
When $x(t)$ is a sinusoidal signal, say $\cos(2\pi f_s t)$,
its integral is
$$\int_0^t \cos(2\pi f_s \tau)\,\mathrm d\tau
= \left.\frac{1}{2\pi f_s}\sin(2\pi f_s\tau)\right\vert_0^t
= \frac{1}{2\pi f_s}\sin(2\pi f_s t)$$
is also a sinusoid at the same frequency, and so for
this special case, the FM signal can be written as

  $$A \cos\left( 2\pi f_c t + \frac{f_\Delta}{f_s} \sin(2\pi f_s t)\right)
= A \cos\left(2\pi f_c t + k\sin(2\pi f_s t)\right).$$

Note that the instantaneous frequency of the FM signal
varies (in sinusoidal fashion) between $f_c-k$ and $f_c+k$ Hz.
Note also the difference between this standard result
and what the OP claims is the formula in his textbook
(I have not read the text myself and have no idea what
it actually says):

  $$v(t) = 8 cos (2\pi [10^5 + 4 sin(5 \times 10^3 t)]\,\mathbf t\,)$$

The extra $\mathbf t$ at the end changes the instantaneous frequency
to
$$2\pi10^5+8\pi sin(5 \times 10^3t) + 8 \pi \, \mathbf{t}\,  (5 \times 10^3) cos (5 \times 10^3 t)$$
(as correctly calculated by the OP).  But,
instead of the instantaneous frequency fluctuating between
fixed limits, the maximum frequency deviation increases with time.
The claim that this is an FM signal (in the
usual sense of the word) is certainly not tenable.

If OP Samuel has correctly copied the problem from
his textbook, then the problem, as stated, has no solution.
Perhaps it is just a typo in this particular problem.
On the other hand,
if the same "additional" $t$ occurs in
occurs in the text in the definition of an FM signal,
and this formulation is used throughout the book, then
the book certainly has a major blunder in it.


How can I filter a noise with pattern type in audio?

Simple idea: let some pros do the work. Send the audio through a cellphone to another one via a WhatsApp, Duo or whatever voice call app, even classic network call works, dealer's choice. The noise reduction algorithms in cellphones are designed to do exactly what you want to do. Reduce noise while conserving human speech. The outcome will depend to some degree on the market, the phone's firmware was designed for. US customers like it, when the background noise is cut out completely, accepting some degradation of speech. European customers prefer more naturally sounding outcomes with some noise left.

I think this would be worth a shot.


Best way to implement variable phase shift on FPGA?

The allpass you're describing is a multiplication with multiple numbers. Hence, if you want to shift your phase, the only mathematically viable and easiest way is indeed a multiplication with $e^{i\phi}$.

I don't see why that would be overly slow – in a modern FPGA, you'd probably utilize DSP slices to do the (real-valued) two multiplications and additions that make up a complex $\cdot$ complex multiplication.

Now, if you need to increment your phase, so you need to *calculate the value of $e^{i\phi}$ first, then there might be room for optimization. In fact, digital tuning happens in FPGAs, in the DSP chains of software defined radio devices like the Ettus USRP series.

There, a pretty standard CORDIC algorithm was employed to approximate a complex sinusoid $e^{i2\pi f}$ of arbitrary relative frequency $f$. Have a look at the  USRP X3xx's GPL-licensed RX DSP DDC chain:

cordic_z24 #(.bitwidth(cwidth))
  cordic(.clock(clk), .reset(rst), .enable(run),
     .xi(to_cordic_i),. yi(to_cordic_q), .zi(phase[31:32-zwidth]),
    .xo(i_cordic),.yo(q_cordic),.zo() );

The cordic_z24 is really just a multi-stage CORDIC implementation with pre-computed atan coefficients.

It runs at complex 64MS/s even on 2005 Spartan FPGAs, and doesn't really blink much of an eye about 200MS/s on modern Kintex7 FPGAs.

It avoids explicit multiplication alltogether.


I think I can come up with two suggestion to address this problem.

First (simple)

You can simply use HSV color space, and try to find right thresholds. It seems in HSV colour space, the car is pretty distinguishable from the background but the rail is the problem now.

I=imread ( your image );
J=rgb2hsv(I);
imshow(0.1<J(:,:,2) & J(:,:,2)<0.2 & 0.2<J(:,:,1) & J(:,:,1)<0.54,[]);

Second (a bit more complex)

Since your background is presumably static, you can remove background and the remaining would be your car (and of course the rod). I suggest to use optical flow for this purpose.


Depends on the amount of clipping and what other information about the signal you have. It's pretty easy to detect clipping (wave form constant at +xmax or -xmax). Once a region has been found you can try to replace the clipped region with an interpolation from the non-clipped regions around this. The best interpolation scheme will depend on what else you know about the signal. You can try polynomial expansions, Taylor series, Frequency domain interpolation, etc.


To elaborate the point of my comment:

What you did there is a kind of circular convolution, where you thought of the impulse response $h(t)$ as periodic, i.e. the system triggered by $\delta(0)$ would oscillate until infinity. The periodicity is represented by your "wrapping around" the impulse response. The output you calculated is actually just one period of the output signal.

This is different from the linear convolution that the task formulation obviously expects to be applied here. With the linear convolution, the impulse response is thought to be the graph representing $h(n)$ padded by zeros from $-\infty$ to $0$ and from $3$ to $\infty$. Here, if the system is triggered by $\delta(0)$, you would get $h(n)$ and then all zeros to infinity.


Why does a system with no poles have a finite support?

If a system doesn't have finite poles, then its transfer function is of the form:

$$H(z) = \frac{Y(z)}{X(z)} = a_Nz^N+a_{N-1}z^{N-1}+...+a_1z+a_0$$

So if you go back to the time domain:

$$y[n]=a_Nx[n+N]+a_{N-1}x[n+N-1]+...+a_1x[n+1]+a_0x[n]$$

Note that if $x[n]=\delta[n]$, then $y[n]=h[n]$, the impulse response of the filter.

$$h[n]=a_N\delta[n+N]+a_{N-1}\delta[n+N-1]+...+a_1\delta[n+1]+a_0\delta[n]$$

It's easy to see that the equation above will return a $h[n]$ that is non-zero for a limited range of values of $n$.

This means that the filter is FIR. Therefore, it has a finite support and it is BIBO stable, leading to your second question.

  Why if the system has a finite support it means that it is BIBO stable?

For BIBO stability of an LTI discrete-time system, one must check that the following condition is satisfied:

$$\sum_{n=-\infty}^{\infty}|h[n]|<\infty$$

If $h[n]$ is the impulse response of a FIR filter, then it has some length $L$:

$$\sum_{n=-\infty}^{\infty}|h[n]|=\sum_{n=K}^{L+K-1}|h[n]|$$

Due to this summation being finite, it converges. Ergo, any FIR filter is BIBO stable.


Shifting a fourier spectrum by subpixel amount in python

You'd not do it in frequency domain at all. Use the shift property of the discrete Fourier transform!

Simply multiply the rows of your original image with an $e^{j2\pi \Delta f_x x/W}$ pointwise ($x$: pixel index in that row, $W$: width of image) before transforming to frequency domain to shift by $\Delta f_x$ in row direction. Same for shifts in column direction; multiply columns with $e^{j2\pi \Delta f_y y/H}$ pointwise ($y$: pixel index in that column, $H$: height of image). There's no restrictions on the "fineness" of $\Delta f$ in either direction, and any 2D shift can be understood as a shift in row and one in column direction.

If necessary, transform your image from frequency to spatial domain, do the multiplications, and transform back.


Combining two filters (audio processing)

In order to combine two filters, you need to convolve their impulse responses, and you haven't done that in this case.

In fact, according to this code snippet, the notch filter is incorrectly built.

I'm sure once you review it you'll figure out your errors. Good luck!


7 is a very small number of taps and you're expecting very sharp cutoff. Also, the cutoff frequency needs to be normalized such that 1.0 is half of the sample rate. Try a much larger order. e.g.

[b,a] = fir1(300,4/fs/2);

Also beware that the beginning of the filtered output is going to have a big transient so you might try a longer waveform.


How to calculate the PSD from the complex calculated STFT?

In general, if you have complex spectrum and need PSD in dB the mathematical equation is

$$P_{xx} = 20\cdot\log_{10}|X_{x}|,$$
where $P_{xx}$ is your PSD in dB and $X_{x}$ is your complex STFT spectrum (variable Spectro in your case).

In SciPy documentation for scipy.signal.spectrogram is mentioned, that you can compute spectrogram with different modes (‘psd’, ‘complex’, ‘magnitude’, ‘angle’, ‘phase’). You have chosen mode='complex', which returns complex STFT, so you should get PSD from equation above. If you have chosen mode='psd' you would get PSD e.g. $S_{xx}$ (Spectro in your case) in the first place and to convert it to dB your equation is $$P_{xx} = 10\cdot\log_{10}S_{xx}.$$

However, if you want to plot same spectrogram as the MATLAB one I would choose equivalent matplotlib.pyplot.specgram, which is very similar and you can plot PSD already in dB with this tool.


Amplitude calculation of the original signal from its FFT

Whoa, clarification is needed. The peak amplitude ($A_p$) of an $N$-length $f_1$ input sinusoid can only be found using

$A_p = 2M/N$

(where $M$ is the the magnitude of the $f_1$ spectral component) if $N$ samples of that $f_1$ input sinusoid represent an exact integer number of cycles.

If that condition holds then you can measure $A_p$ for the $f_1$ sinusoid as well as the peak amplitudes of the harmonics. But those peak amplitude values apply only to the individual sinusoids that comprise your FFT input signal. You CANNOT use the individual measured peak amplitudes to measure the peak amplitude of your nonsinusoidal FFT input.  That was the point made in jojek's comment.


For ergodic processes, time averages (defined by integrals over time) and ensemble averages (defined by expectations with respect to probability distributions) are identical. This means that the autocovariance is the same, no matter if defined by a time integral or by an expectation:

$$C_X(\tau)=\lim_{T\rightarrow\infty}\frac{1}{2T}\int_{-T}^{T}X(t+\tau)X(t)dt=
E\{X(t+\tau)X(t)\}$$


How to shift a signal to the center of the frequency domain?

I've found that shifting to (513, 513) works best, because in MATLAB, the zero frequency is at the upper left corner of the unshifted Fourier transform.  The operation fftshift flips the 1st and 3rd quadrants, meaning that the zero frequency is now at the top-left of the bottom-right quadrant.


Time domain representation of a modulated signal when symbol rate >> carrier frequency?

tl;dr: you can't. That's basics!
longer story:

You ask about symbols.
When the symbol rate is much larger than the carrier frequency, you've basically lost, and cannot transport that many symbols per second; you don't have the bandwidth; you get at most twice the carrier frequency in bandwidth, namely from 0 Hz to twice the carrier frequency.

However, you're asking how to modulate 10¹⁵ bits on a 1 Hz carrier, and that's a different question: theoretically, you can take forever to transport these 10¹⁵ bits, so there's no problem here. Realistically, we'll all be dead by the time you've finished your transfer, but whether or not that is  a problem is of mere philosophical interest.

If you, however, mean you want to transport 10¹⁵ bits per second on a 1 Hz carrier wave, you have less than two channel access per second, so you'll need enough SNR that you get more than  $2^{5·10^{14}}$ different constellation points.
If you reasonably try this, I hope someone stops you in time: as physics dictates the difference between two detectable constellation points needs to be at least the Planck constant, in fairly good approximation, the action needed to get from the all-0 bitword to the all-1 bitword would be at the very least,
$$\frac{2^{\left(5\cdot 10^{14}\right)}\cdot h}{\pi} =\frac{2^{\left(5\cdot 10^{14}\right)}\cdot 2^{-110.2}\frac{\text{kg·m²}}{\text{s}}}{\pi} \approx 2^{\left(4\cdot 10^{12}\right)} \frac{\text{kg·m²}}{\text{s}}\approx 10^{4\cdot10^{11}} \text{J·s}.$$
That is incredibly much action. This is nuts. Your photon of 1 Hz would have to have an energy of $10^{4\cdot10^{11}}\,\text{J}$ to be able to walk a path of that action; that's as much energy in a single photon as the sun emits in 487 years. This photon would end the universe as we know it. There'd be only heat and chaos. And then, who the hell will care about your 10¹⁵ bits?


Compute the distance transform. You can do that very efficiently. This will give you the Euclidean distance to the closest boundary pixel, as well as the index of that pixel. Then, normalize all the distances in the mask so that the maximum becomes 1. This is required for the next stage: At each point initiate a Gaussian kernel. The kernel size $h$ (and similarly the standard deviation $\sigma$) should be a function of the weight, obtained from distance normalization. You can either use a coefficient - linear proportionality - or a complex kernel, depending on the blurring requirements.

Finally, blur the signal at each pixel using the generated mask. Each pixel gets a unique mask and thus a different blurring. Note that, some of the pixels could get the same blur.

You only need to adjust the kernel size $h$, as for typical Gaussian kernels, $\sigma$ is a function of $h$. The function is simply a lookup table:

   $h$  $\sigma$
    3    (0.65)
    5    (0.87)
    7    (1.43)
    9    (1.88)
    11   (2.31)


Filter response is changed according to its order

An order of 30 is way too high, especially for a band reject filter with such a narrow band as you want.

Your filter will have 15 poles in a very small semicircle between frequencies 58 and 62 Hz (and another 15 poles in the negative frequencies).

Your a coefficients correspond to an order 30 polynomial with the 30 roots corresponding to the poles of the filter.

Note that any minuscule numerical approximation (which you are not free from) will alter the position of all roots. You will very probably get roots ending up outside of the unit circle (which makes your filter unstable).

So the ugly frequency response for the high order filter is expected.

As a rule of thumb, do not try to implement IIR filters of order above 10.


Which one is the correct graph for $u[2-n]$?

HINT: For which value of $n$ does the argument of $u[2-n]$ become zero? That's where the step occurs. For which values of $n$ is the argument non-negative? That's where your unit step equals $1$. If you think about it for a minute, it should become really easy.


Implementation of IIR Filter in dsPIC 33EP - Scaling

The tf2sos function takes an input filter of order
$N$, given by $H(z)=\frac{\sum_i^{N} b_iz^{-i}}{\sum_i^{N} a_iz^{-i}}$ and returns coefficients for $N/2$ second-order filters $H_k(z)=\frac{b_{k0}+b_{k1}z^{-1}+b_{k2}z^{-2}}{a_{k0}+a_{k1}z^{-1}+a_{k2}z^{-2}}$ and a gain $g$ such that

$$
H(z) = g \prod_kH_k(z)
$$

So, regardless of what normalization you want from the function (none, inf or two), the above always holds. So to answer your questions:

the gain of each filter section is different (i.e. not necessarily one). (But, you dont care about that)
you do not need to adapt the gain between the filter stages. What comes out of the first stage is directly sent to the second stage. Then, after all stages were processed, you multiply the outcome of the last stage with $g$, what you receive from the tf2sos function.


Are Total Harmonic Distortion and Noise (THD+N) and THD correlated in practice?

Q1 Are %THD+N and %THD correlated in practice in the case of common
  audio devices (speakers, headphones, amplifiers)?

This really depends. THD is just a contributor to THD+N. If the noise dominates the numbers can be vastly different and if the THD dominates the numbers are basically the same. Passive components (passive speakers or headphones) have no inherent electrical noise so THD tends to dominate unless you hit a point where you trigger mechanical or acoustic noise as well. A good example of non-harmonic noise is port noise on a vented speaker.

  Q2 How important is %THD in the reproduction of the input signal in
  speakers or headphones?

Depends on your definition of important on what you care about, but specifically in speaker design we keep a very close eye on THD since speakers are fairly non-linear.


My example isn't in fact below $f_s/2$; I blundered with $f_s=N$. Adjusting this eliminates aliasing, for linear, exponential, and hyperbolic chirps - and, in fact, for any signal with $\phi '(t) < f_s/2$.
This can be verified two ways:

For $\omega_\text{max} = \text{max}(\phi'(t))$ (assume $\phi(t)$ includes $2\pi$), the worst-case scenario is a straight line with slope $\omega_\text{max}$ (a curve can't go any steeper else it's no longer "max", and flatter is "more favorable"). But we know if $f_\text{max} < f_s/2$, a pure tone won't alias, thus the worst case is accounted for.

Aliasing occurs if the finite difference $\phi[n] - \phi[n-1] = \phi(n\Delta t) - \phi((n-1)\Delta t) \geq \pi$ for any $n$, where $\Delta t=$ sampling period, or increment of the t vector.  The greatest such difference is $\omega_\text{max} \cdot \Delta t$ - setting the inequality: $\omega_\text{max} \Delta t < \pi \Rightarrow \omega_\text{max} < \pi / \Delta t \Rightarrow f_\text{max} < 1/(2\Delta t)= f_s/2$.

Thus if I ensured $f_s/2$ was met, phi could never exceed $\pi$.

Example in question, corrected: if we seek (tmin, tmax) = (0, log(511)), our $N$ must be large enough so that $\Delta t$ is small enough for $\text{max}(\phi'(t)) < \pi$. The max is 2*pi*exp(log(511)) = 2*pi*511 = wmax, so max finite difference is wmax*dt < pi -> dt < pi/wmax, and dt = t[1] - t[0] = (tmax - tmin)/N, thus (tmax - tmin)/N < pi/wmax -> N > wmax/pi*(tmax - tmin) = 2*pi*511 / pi * log(511) = 2*511*log(511) = 6373.6.
N = 6374
t = np.linspace(0, np.log(511), N, endpoint=False)
phi = 2*np.pi * (np.exp(t) - 1)
x   = np.cos(phi)


MATLAB - remove the frequency at zero in FFT

Frequency bin at zero is simply mean value of your signal. Just take a look on definition of DFT, for zero frequency $k$ we get:

$$\left. X[k]=\sum_{n=0}^{N-1}x[n]e^{-i 2\pi  n\cdot k} \right |_{k=0} \Rightarrow X[0]=\sum_{n=0}^{N-1}x[n]e^{-i 2\pi  n\cdot 0} = \sum_{n=0}^{N-1}x[n] $$

So after you calculate the FFT and divide by number of samples, indeed you get the mean value of your signal.
In order to remove that, please apply the following before the FFT calculation:

Lift = Lift - mean(Lift);