Custom nearly-linear-phase IIR filter computation gone wrong

I believe that there's something going wrong with the implementation of the filter. The design algorithm lets you prescribe the maximum pole radius, so if you prescribe e.g. $r_{max}=0.98$ the filter is guaranteed to be stable with some reasonable stability margin. So using floating point arithmetic it seems to me that nothing can go wrong. If you use Matlab's filter() function, be careful to call it filter(b,a,x) with b being the numerator coefficients. I've seen people use this function as filter(a,b,x) (which seems more logical to some folks). The latter function call will give you trouble, because filter() will interpret the numerator coefficients b as denominator coefficients, so any zero outside the unit circle will become an unstable pole. And due to the phase approximation of the design algorithm there will almost always be some zeros outside the unit circle.

EDIT:

This is example #7 from this site, with the desired group delay value set to 40 samples (instead of 37). This is what I understood from your comments you're trying to do. This is the code:

tau=40;M=40;N=6;r=0.98;
om=pi*[linspace(0,0.2,20),linspace(0.23,1,77)];
D=[exp(-1i*om(1:20)*tau),zeros(1,77)];
W=[ones(1,20),1000*ones(1,77)];
[b,a]=mpiir_l2(M,N,om,D,W,r);

[H,w] = freqz(b,a,1024);
gd = grpdelay(b',a',1024);
subplot(2,1,1), plot(w/2/pi,20*log10(abs(H)));
grid, title('|H(f)|')
subplot(2,1,2), plot(w/2/pi,gd);
grid, title('passband group delay'), axis([0,.1,35,45])

I ran the design program in Octave. I believe that you have a problem with running the software in the current Matlab version. Note that the software was written about 17 years ago, and I'm not sure if and how you (or someone else) adapted it to the current version. You can find the filter coefficients below and I'm convinced that a filter implemented with these coefficients will work fine.

These are the coefficients:

b =

  -5.3754e-04
   1.4765e-03
  -1.7103e-03
   1.1605e-03
  -3.7906e-04
   6.8635e-05
  -2.9013e-05
  -2.0759e-05
  -4.4175e-05
  -1.0035e-05
  -8.8312e-06
   3.8369e-05
   3.4096e-05
   5.6969e-05
   1.6945e-05
   5.8498e-06
  -5.1954e-05
  -5.4722e-05
  -7.7327e-05
  -2.8043e-05
  -1.5210e-07
   7.7270e-05
   9.5478e-05
   1.2200e-04
   5.9885e-05
   1.6209e-06
  -1.1856e-04
  -1.8065e-04
  -2.3618e-04
  -1.7031e-04
  -5.9434e-05
   1.7262e-04
   4.0741e-04
   6.8756e-04
   8.7908e-04
   9.9248e-04
   1.3572e-03
  -1.4580e-04
   3.0311e-03
  -1.8688e-03
   1.9026e-03

a =

    1.00000
   -4.65499
    9.58904
  -11.07801
    7.53566
   -2.85503
    0.47033


Digital to analog aliasing or mirror query... DAC can output negative frequency?

I think you are confused by negative frequencies and what they mean so let me add this explanation.

When you see a spectrum that contains "positive" and "negative" frequencies, each of the frequencies are of the form:

$$e^{j\omega t}$$

Where $\omega$ is the frequency (in this case angular frequency as $2\pi f$ with f being the frequency in Hz.

The general form of $e^{j\phi}$ is a phasor on the complex plane with magnitude one and angle $\phi$, so $e^{j\omega t}$ is simply a phasor with unit magnitude rotating with rate $\omega$. If $\omega$ is positive, the rotation is counter-clockwise: a positive frequency, and if the rotation is negative, the rotation is clockwise and represents a negative frequency. So the concept of positive and negative frequencies only applies to frequencies of the form of the complex exponential and is the direction of rotation.

Consider Euler's identity for a cosine that demonstrates the concept of negative frequencies well:

$$cos(\omega t) = \frac{1}{2}e^{j\omega t} + \frac{1}{2}e^{-j\omega t}$$

Here we get two phasors on the complex plane, one rotating clockwise, one counter-clockwise, and adding the two results in a signal that always stays on the real axis going back and forth sinusoidally.

So therefore a simple cosine wave with a constant frequency would contain one positive (exponential) frequency and one negative (exponential) frequency. We also see from this that ANY real signal must have positive and negative frequencies, AND the negative frequencies MUST be of the opposite phase. Like the cosine example above, each phasor component spinning in a counter-clockwise direction (the positive frequencies) MUST have an equal magnitude phasor spinning in the opposite direction (negative frequency component) in order for the sum to stay on the real axis (which would then be a real signal with no imaginary components. This is said to be "complex conjugate symmetric" when the negative frequency components have equal magnitude and opposite phase compared to the positive frequency components.

In order to IMPLEMENT a signal that does not have a spectrum with this property (positive and negative frequencies not complex conjugate symmetric), you would need to have TWO values in implementation, one to represent the real number and one to represent the imaginary number (or you could do magnitude and phase, either way you need two signals). You see this in the baseband part of your diagram where the real is represented as "I" for "in-phase" and the imaginary is represented as "Q" for quadrature. So in that part of your diagram we are indeed representing a complex signal, and the spectrum at that point can be asymmetric.

I notated your diagram to demonstrate this further. I will call the IQ signal the "baseband signal" and because it is indeed complex, it's spectrum can be assymetric as I have drawn. This complex baseband signal is then digitally upconverted (frequency translated) by multiplying it by $e^{-j\omega t}$ since  $e^{-j\omega t} = cos(\omega t)-jsin(\omega t)$ and then taking the real part:

$$(I_1+jQ_1)(I_2-jQ_2)$$
$$= (I_1I_2+Q_1Q_2)+j(I_1Q_2-I_2Q_1)$$

Taking the real part of this results in the form in your block diagram prior to the DAC:

$$(I_1I_2+Q_1Q_2)$$

Given as:

$$I cos(2\pi f_{IF}t) + Q sin(2\pi f_{IF}t)$$

With $I_1$ and $Q_1$ as the I and Q inputs to the digital multiplier and $I_2$ is the cosine waveform into the upper digital multiplier and $Q_2$ is the sine waveform input into the lower digital multiplier in the diagram)

So the signal where we maintain two datapaths (truly a complex signal) will have an asymmetric spectrum while the signal where this is only one datapath must be real and therefore must have a spectrum that is complex conjugate symmetric. I like this diagram as shown as it clarifies that the first signal I circled is indeed complex and has an asymmetric spectrum, while the second one I circled, while it looks like the spectrum is asymmetric if you just considered the positive frequencies alone, is not as we are referring to symmetry between the positive and negative frequency axis.

So we see how negative frequencies exist throughout the entire block diagram, and it is not that the DAC is creating them.  The baseband I and Q signal coping out of the encoder is complex and has an assymetric spectrum, this gets digitally upconverted (and is still complex until summed), the output of the summer is now one signal path so must be real and therefore now has a symmetric (complex conjugate) spectrum.  The DAC will convert this to the analog domain and the bandpass filter will filter out the images of the DAC (positive and negative--it is a real filter so it too is complex conjugate symmetric in its passband). So assuming the BPF is centered on $f_{IF}$, the analog spectrum would look identical to the spectrum prior to the DAC, extend to $\pm \infty$ while the unique DAC input spectrum extends only to $\pm f_s/2$.

The "image" shown  in the diagram are NOT the images due to digital sampling, but the image due to the analog upconversion. This is a real multiplier and as such its output will contain the sum and the difference of the input (review the expansion of $cos(\alpha)cos(\beta)$ to see this), so we get the terms RF+IF and RF-IF at both postiive and negative frequency locations. (However you can also see that we are just frequency translating the same spectrum at the output of the BPF by convolving it with the two frequencies of the RF sinusoid.)

There is more not detailed about the DAC upconversion that I didn't cover but you should be aware of including the option to pass a higher Nyquist zone by setting the bandpass filter on one of the DAC images, and that the DAC itself will introduce a Sinc roll--off in frequency due to the staircase reconstruction.


Relationship between time average and cross spectrum

I will infer from the "discrete Fourier coefficients" in your question that we are considering periodic signals, otherwise we should be considering a continuous spectrum (the Continuous-Time Fourier Transform). I will also be using the definition
$$\left< x(t), y(t) \right>_t = \frac{1}{T}\int_{\alpha}^{\alpha +T} x^*(t) y(t) dt$$ where the superscript $*$ denotes complex conjugation, $T$ is the shortest period that both $x(t)$ and $y(t)$ share, and $\alpha$ is an arbitrary real number. If these definitions are not correct, please alter your question and provide the correct ones.

We know that we can decompose $x(t)$ and $y(t)$ using a Fourier series expansion
$$ x(t) = \sum_{n=-\infty}^{\infty} c_{x,n} \exp\left(\frac{2 \pi n t}{T}\right)$$
and
$$ y(t) = \sum_{n=-\infty}^{\infty} c_{y,n} \exp\left(\frac{2 \pi n t}{T}\right).$$

We know from Parseval's Theorem that
$$ \sum_{n=-\infty}^{\infty} c^*_{x, n} c_{y, n} = \int_{-\frac{T}{2}}^{\frac{T}{2}} x^*(t) y(t) dt.$$
Notice the equality. Also notice that the right hand side is merely the product of $T$ and $\left< x(t), y(t) \right>_t$. This means that whenever the left hand side is zero, so is $\left< x(t), y(t) \right>_t$ and whenever $\left< x(t), y(t) \right>_t$ is zero so is the left hand side. I believe this resolves your question. Let me know if you have objections, and I will try to address them.


Intuition behind impulse response terms in convolution

The coefficient $h[n]$ is the value of the system's response at time $n$ when the input signal was an impulse at $n=0$. Obviously, that's why we call $h[n]$ the impulse response. From this you can see that for a system to be causal, $h[n]$ must be zero for $n<0$, otherwise the system would "know" in advance that an impulse will come at $n=0$.

Note that any discrete-time signal $x[n]$ can be written as a sum of unit impulses:

$$x[n]=\sum_{k=-\infty}^{\infty}x[k]\delta[n-k]\tag{1}$$

Since the response to a shifted impulse $\delta[n-k]$ is $h[n-k]$, and since the system is linear and time-invariant, from $(1)$ the output signal must be given by

$$y[n]=\sum_{k=-\infty}^{\infty}x[k]h[n-k]\tag{2}$$

which is of course the discrete-time convolution. Consequently, positive indices of $h[n]$ correspond to the memory of the system, that's why $h[1]$ is multiplied with the past input sample $x[n-1]$, etc.


If $g(x)= f(N-1-x)$ then
$$\begin{aligned}
  (\operatorname{DCT}\:g)(u)
    =& c(u)\cdot\sum_{x=0}^{N-1}g(x)\cdot\cos\bigl((2x+1)u\cdot\tfrac\pi{2N}\bigr)
  \\=& c(u)\cdot\sum_{x=0}^{N-1}f(N-1-x)
        \cdot\cos\bigl((2N-2N+1+2x)u\cdot\tfrac\pi{2N}\bigr)
  \\=& c(u)\cdot\sum_{(N-1-x)=(N-1)}^{0}f(N-1-x)
        \cdot\cos\Bigl(\bigl(2N-1-2(N-1-x)\bigr)u\cdot\tfrac\pi{2N}\Bigr)
\end{aligned}$$
substitute now $N-1-x=:\xi$,
$$\begin{aligned}
  (\operatorname{DCT}\:g)(u)
    =& c(u)\cdot\sum_{\xi=0}^{N-1}f(\xi)
        \cdot\cos\Bigl(\bigl(2N-1-2\xi\bigr)u\cdot\tfrac\pi{2N}\Bigr)
  \\=& c(u)\cdot\sum_{\xi=0}^{N-1}f(\xi)
        \cdot\cos\Bigl(\pi\bigl(\tfrac{2N}{2N}-\tfrac{2\xi+1}{2N}\bigr)u\Bigr).
\end{aligned}$$
In the cosine, we have now a shift by an integer multiple ($u$) of $\pi$. Such a shift in a sine or cosine function is equivalent to multiplying the function with $(-1)^u$, So the result is
$$\begin{aligned}
  (\operatorname{DCT}\:g)(u)
    =& c(u)\cdot\sum_{\xi=0}^{N-1}f(\xi)
        \cdot\cos\bigl(\pi\tfrac{2\xi+1}{2N}u\bigr)\cdot(-1)^u
  \\=& (-1)^u\cdot(\operatorname{DCT}\:f)(u).
\end{aligned}$$
Not quite the expression you wanted to show, but this should be the correct one.


Inverse wavelet transform on wavelet coefficients

Because, as far as I know, there is no discrete wavelet scheme related to the morlet wavelet, admitting neither orthogonal nor biorthogonal implementation. See for instance the discussion in Complex Morlet function and DWT.

Using waveinfo('morl'), you indeed get that: DWT: no

 Information on Morlet wavelet.
    Morlet Wavelet
    Definition:
    morl(x) = exp(-x^2/2) * cos(5x)
    Family                  Morlet
    Short name              morl

    Orthogonal              no
    Biorthogonal            no
    Compact support         no
    DWT                     no
    CWT                     possible

    Support width           infinite
    Effective support       [-4 4]
    Symmetry                yes

So you might obtain something only (and I am not sure it is valid for your data) using some form of inverse continuous wavelet transform, for instance in Matlab  icwt.


You may have problems related to discontinuity and this can give small glitches in your signal, if you hear your signal you will notice pops/clicks artifacts at every 100 epoches edge, one way to try avoid this problem is apply one window function in a overlapped segments between your Epoches.

When you apply a window like a Hann window this reduces the amplitude of the components of the outlet ends of the window and its suppress the edge artifacts when reconstructing your signal to be continuos, for it you need use technics like overlap and add to works like a crossfading between every overlapped segments.

OK you asked for example, i can show you the basic concept, imagine two sinusoids at 100hz, if you try join you will have problems, take a look:

In this plot is clear that the end of the sine wave does not match with the beginning of another, and this cause the problem described above.

One way to do try avoid it is apply a crossfade, here an example:

Fs=22050;
f=100
winsize=2048;

%creating signal
signal1= 0.9*sin(2*pi*f/Fs*(1:winsize));
signal2= 0.9*sin(2*pi*f/Fs*(1:winsize));

%use plot(in) to see the first image plot in this answer
in=[signal1 signal2];

%create a hann window
window1=hann(winsize)';
window2=hann(winsize)';

%i'm using just the last half values of the hann window to skip fade in
window1(1:winsize/2)=ones;

%i'm using just the first half values of the hann window to skip fade out
window2(winsize-(winsize/2):winsize)=ones;

%apply window in signal1 and signal2
teste1=(signal1.*window1);
teste2=(signal2.*window2);

%apply 50% crossfade between signal1 and signal2
OverlapCrossfade=teste1(length(teste1)/2:length(teste1)-1) +  teste2(1:length(teste2)/2);

%rebuild the signal after crossfade
join = [teste1(1:length(teste1)/2) OverlapCrossfade teste2(length(teste2)/2:length(teste2)-1)];

Now if you do plot(join) you can see it:

OW this is better i cant see glitches (but we can see amplitude variation), now all that you need to do is apply some type of envelope adjust.


orthogonal family and pulse shaping filter

In general for a pulse shaping filter and any orthogonal code family: No.

UPDATE: See @Engineer's answer which I confirmed by excluding the first row in the set of Walsh codes (the all 1 row) the remaining rows which have a balance of 1's and 0's are indeed all orthogonal.

*SECOND UPDATE: @MBaz has pointed out here that with Root-Raised Cosine filtering the entire code set will remain orthogonal: Orthogonal Codes for Band Limited Channel

A detailed mathematical proof along with the criteria where orthogonality can be maintained would be a much better answer, but since no one has submitted one I am at least providing here a simulation result that demonstrates a specific example of an orthogonal family where after pulse shaping with a raised cosine filter the orthogonality was not maintained. (The original signal can of course be extracted at the sampling locations with no distortion given the zero-ISI property of a raised-cosine filter, so that extracted signal would still be orthogonal after sampling, but the question was specifically on the orthogonality of the pulse shaped waveform itself ($g(n)*s_n(n)$)

The formal definition for orthogonal signals is that their inner product (also called the dot product) is zero.

Consider a simple case of the Hadarmard $H_4$ matrix, with each row forming a set of four orthogonal Walsh codes:

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

Before pulse shaping it is easy to confirm orthogonality, for example between row 1 and row 2 using the dot product of index by index product and sum:

$$(1\times1)+(1\times-1)+(1\times1)+(1\times-1) = 0$$

And likewise betweeen row 1 and row 4:
$$(1\times1)+(1\times-1)+(1\times-1)+(1\times1) = 0$$

The result of a simulation of this set with a raised cosine filter suggests that orthogonality is maintained in the case of the first two rows, but is not between row 1 and row 4. The resulting waveforms from the first two rows after upsampling by 10 and pulse shaping with a raised cosine filter are shown in the figure below (the sample locations for decoding the original signal are at samples 50, 60, 70 and 80):

The computed dot product for these two waveforms was $4.4E-16$. (Close enough to zero so orthogonal).

Row 1 vs Row 3 is similarly quite low at $9.1E-16$ (again orthogonal).

However
Row 1 vs Row 4 was $-0.465$ which is a very strong correlation and clearly no longer orthogonal.

The plot showing the pulse-shaped Row 1 and Row 4 is shown below.

The dot product is the sample by sample product and then the sum (or integration) of the result. Plot of the product of each of the waveform pairs gives further insight into the loss of orthogonality. Note how the product of the pulse-shaped row 1 and row 2 is very symmetric with similar positive and negative areas under the curve, while the product or the pulse shaped row 1 and row 4 has a significantly larger negative area, thus destroying orthogonality in that case:


Binomial Approximation of Gaussian Distribution

Let $X \sim B(n, p)$. At the third layer $\begin{bmatrix}1 & 3 & 3 & 1\end{bmatrix}$ of the Pascal's triangle, we have $X \sim B(3, \frac{1}{2})$. Which means $P(X = k) = {3 \choose k} (\frac{1}{2})^k (\frac{1}{2})^{3-k} = \frac{{3 \choose k}}{2^3}$. The variance is given by $\sigma^2 = n p q = 3 \cdot \frac{1}{2} \frac{1}{2} = \frac{3}{4}$ and the std by $\sigma = \sqrt{\frac{3}{4}}$. You can derive the variance by looking at $n$ independent Bernoulli trials.

In general, a 1d binomial filter/kernel of size $n$ with $p = \frac{1}{2}$ is given by $\begin{bmatrix}P(X = 0) & P(X = 1) & \cdots & P(X = n-1)\end{bmatrix}$. As $n \to \infty$, this kernel will approximate the normal distribution with mean $\frac{n}{2}$ and variance $\sigma^2 = \frac{n}{4}$. If you want a different variance/std, you should change $p \neq 0.5$.

As somebody already wrote in the comments, you should look at the De Moivre-Laplace theorem to prove this. However, this is only a special case of the more general central limit theorem. By showing that the binomial distribution is asymptotically equivalent to the normal distribution as $n\to\infty$, I would argue that we have shown already both directions.


Question About $ Q $ Matrix (Model Process Covariance) in Kalman Filter

I'm the author of the textbook linked to above. This is a new account and thus I am not allowed to reply to that answer.

Anyway, the Gaussian chapter covers the definition of a covariance matrix. In general though your Q matrix will be full, not a diagonal, because there is correlation between the state variables. For example, suppose you are tracking position and velocity in x-axis for a car. Process noise (which Q is modelling) might be things like wind, bumps in the road, and so on. If there is a change in velocity due to this noise, there will also be a change in position. They are correlated, and so the off-diagonal elements will be non-zero for those terms.

The Kalman math chapter tells you how to compute the Q matrix for various scenarios. Admittedly there is more work to be done on that section. It's one of the more opaque subjects as most sources either gloss over it or just throw an integral at you and expect you to figure out the rest.

Finally, the companion FilterPy project contains code to compute the full Q matrix for you given common scenarios. The code is in common/discretization. For analytically difficult situations there is a routine to find the result numerically following an idea from Golub and Van Loan; the textbook is lacking on this topic, however.


`SSIM(YUV)` vs. `SSIM(Y)`: Values Are Very Different

I figured it out: the ssim() function in MATLAB, does not handle RGB images, but instead threats it as a 3D-image (with the third dimension time for example, to measure the quality of a video). So it interpreted my 3-color-component image as a video with 3 frames.


Exact average error probability of QPSK in Rayleigh fading channel

From Tse and Viswanath, "Fundamentals of Wireless Communication" (which I highly recommend): $$p_e = \frac{1}{2}\left(1-\sqrt{\frac{\text{SNR}}{2+\text{SNR}}}\right)$$ Note that this is the probability of symbol error.


Is what I am doing correct?

Yes

If not, how would I reproduce the SINAD function?

By replicating the exact same algorithm that the Matlab function uses.
Your answer is more correct, but you also "cheated". You have access to the original sine and noise signals. sinad() ONLY has access to the sum of both. Hence it needs to do an estimation using windowing and Periodogramm analysis. That's less precise than calculating it directly from the original signals. If you really want to replicate what sinad() is doing, try an algorithm that only uses the combined noise+sine signal.


Filtering can be done in the frequency domain which is actually a very efficient technique (and it can very well be, that Matlab does this internally). However, for very long signals it's not as straight-forward as "taking the FFT and applying a frequency response". You can read up on overlap-add filtering which is a possibility to filter in frequency domain. However, this is mainly done because it can be faster.

However, this question is independent from filter design which is a completely different thing. Every FIR filter has a corresponding frequency response and it does not matter whether you convolve in time domain or apply the frequency response in frequency domain. So your approach becomes a question of filter design and the consequences it has.

Just setting unwanted frequencies to zero might completely eliminate those frequencies but it usually comes at the price of a significant ringing in the time domain which is usually unwanted (Wikipedia on ringing). So in fact, what you are proposing is just one way to design a filter and frankly not a good one in most cases.

By allowing a transition band in which the frequency response can gradually go from passband to stopband, degrees of freedom are gained that can be used to improve other properties of the resulting filter (for example eliminate the ringing or obtain a shorter impulse response). That's why Matlab implements so many different filter types, they all have different properties and selecting the most suited one is actually part of designing your signal processing system.

This topic is actually quite complicated and I suggest to read up on filter design. Whole books can be filled with this.


Rayleigh fading channel as a scalar quantity

Rayleigh channel is defined between transmit and receive antennas, if there is one transmit and one receive antennas then there is one Rayleigh fading coefficient.

If there is one transmit and multiple receive antennas then there is a vector or fading channel.

If there is one receive  and multiple antennas then there is a vector or fading channel.

If there are multiple trasmit and receive antennas then you have a Rayleigh fading matrix represting the wireless channel


how to sample and modulate a gps gold code

Have a look to the definition of the C/A signal here: where it says t you need to put your time, where it says i, you need to put your index i  from the for loop, so that concept is correct.

I think your confusion comes from using the frequency in the equation instead of the sampling time which is your time vector.


Band-edge timing error detector, time accumulator

Thanks for your hints and answers. That helped me.
Ive got my code to work now. Its not fine tuned but here is my solution:
-As loop filter i use a PI Filter with a relatively low Ki.
The loop filter is just updated at a 1 Sps/Sym rate when the time accumulator gives an overflow.
-The time accumulator is a modulo interpolation controller
-The interpolator is a farrow structure for cubic interpolation
The IQ plot isn't stationary because the phase isn't locked. But i wont add a following pll anyways because it is a pi/4-DQPSK signal. A bit of phase shifting shouldn't bother me


I have seen people use MFCC coefficients for speech related analysis.It is a  done using filter banks and logarithmic scales for representing speech.


Pulse shape filtering is used to constrain the signal bandwidth (with what would otherwise be a Sinc function in frequency given the rectangular pulse) and as a Nyquist filter done in such a way as to not introduce inter-symbol interference.
For CDMA and specifically direct-sequence spread spectrum in many cases we have no concern for spectral containment (GPS for example since by the time the signal reaches any receivers it is well below the noise floor and therefore not an interference to other receivers in adjacent bands). But IF we did, then it would definitely be the individual chips that should be pulse shaped, since it is their base pattern that would define the spectral occupancy.


Signal Acquisition in Compressed Sensing

Just to restate everything for clarity, the Compressed Sensing problem is defined as the following: given a signal $x$ of length $N$, we measure the projection of $x$ by some projection operator, $\Phi$ of size $M \times N$,

$$ y = \Phi x,$$

where $y$ is are our $M$ measurements of the original signal. We could also state this in terms of the inner products of $x$ with the rows of $\Phi$,

$$ y_i = \left< \phi_i, x\right>.$$

In the CS sampling procedure we have an inherent dimensionality reduction since we are taking a projection of $x$. The degree of this dimensionality reduction is usually referred to in terms of the ratio $\frac{M}{N}$ (aka the subsampling ratio or "Subrate"). By obtaining $y$ during acquisition we are simultaneously performing acquisition and dimensionality-reduction rather than full-resolution sampling followed by compression (DCT, DWT, etc.). The measurements $y$ are read off of the sensor, quantized (still an open problem, but decent results can be obtained through simple scalar quantization. See Laska's dissertation for more novel approaches), entropy coded (pick your flavor), and then transmitted or stored.

You touch on a good key point when it comes to CS signal acquisition, namely, if the calculation of this projection is costly then the advantage of CS for acquisition systems would appear small. However, depending on the context and type of signal, such projections can be accomplished in the analog domain (requiring no computation). The single-pixel camera (SPC) is a great example of this.

Researchers have also been able to greatly reduce sampling times for MRI using CS techniques. Specifically, conventional MRI techniques sample along different radial lines within the frequency domain. Each radial line is a measurement of the MRI device requiring some sampling time. Traditionally, for higher resolution MRI, more radial lines must be acquired, incurring longer acquisition time (problematic for the MRI of small children). However, this process can also be characterized as linear projections within the frequency domain. Because of this, CS techniques can be used to recover MRI images from many fewer measurements, allowing for high resolution MRI with much shorter MRI appointments. This was, in fact, the context in which CS was first applied. For more info, this paper by Lustig et al might be a good starting point.

In summation, the usefulness of CS to a particular context really depends on how you implement the projection. Getting it right can require some out-of-the-box thinking to come up with a novel sensing strategy. Thankfully, frameworks such as the SPC are generalizeable to a wide range of different signal contexts, so we don't have to re-invent the wheel every time.


Need to calculate the VARIANCE of GLCM matrix of an image in MATLAB

This link has details code of how to calculate GLCM variance and other statistical features in matlab.
Reference link https://www.imageeprocessing.com/2019/04/texture-measures-from-glcm-matlab-code.html
For meanX and meanY calculation


Welch’s Algorithm of Averaging Modified Periodograms

Understanding the Windowing Method in PSD Calculation

Please go through the link .. Someone has posted their code snippet which matches with MATLAB output. Hope that might give you better understanding


I suspect you're not doing symbol synchronization. You need to synchronize and then sample the output of your receiver filter at the correct times. Only then will you see the tradiational constellation diagram. Look for "Gardner algorithm" on this website; another well-known algorithm is Mueller and Muller.

One thing you can do to get confidence that your signal is correct is to plot its eye diagram. If the eye is open, then all you need to do is synchronize and sample. If it's closed, you have a different problem.


channel impulse and frequency response: H?

No, $H_{rx,tx}$ is the channel matrix with $(r,t)^{th}$ element representing channels attenuation and phase rotation on transmit symbols $\vec{x}_{tx}$, from $t^{th}$ transmit antenna to $r^{th}$ receive antenna. The equation is neither in time-domain nor in frequency domain. It is on transmit symbols. For example, in OFDM systems, you get $y_{rx}$ (received symbols) only after FFT operation at the receiver, so you can understand the equation as in frequency domain. But for a Non-OFDM based system, like CDMA, you get received symbols $y_{rx}$ after demodulation with matched-filter and that would be in time-domain.

You can understand this equation as a relationship between $rx$ number of received symbols $y_{rx}$ and $tx$ number of transmitted symbols $x_{tx}$ in MIMO scenario where $n_{rx}$ is the AWGN for each of the $rx$ receive chain and $H_{rx,tx}$ is the channel matrix affecting the symbol transmitted from $t^{th}$ antenna and received at $r^{th}$ antenna.
Since symbols are represented as vector points complex plane, hence each channel can only affect each symbol by changing its magnitude and giving a random phase rotation. It doesnot matter time-domain or frequency-domain. That depends on your Receive chain implementation of getting $\vec{y_{rx}}$.


Good features / algorithms for recognizing car-models in images

I would have a look at the so called "bag of words" or "visual words" approach. It is increasingly used for image categorization and identification. This algorithm usually starts by detecting robust points, such as SIFT points, in an image. The region around these found points (the 128 bit SIFT descriptor in your case) is used.

In the most simple form, one can collect all data from all descriptors from all images and cluster them, for example using k-means. Every original image then has descriptors that contribute to a number of clusters. The centroids of these clusters, i.e. the visual words, can be used as a new descriptor for the image.  Basically you hope that the clusters an image its descriptors contribute to, is indicative of the image category.

Again, in the most simple case, you have a list of clusters, and per image, you count which of these clusters contained descriptors from that image and how many. This is similar to the Term Frequency/ Inverse Document Frequency (TD/IFD) method used in text retrieval. See this quick and dirty Matlab script.

This approach is actively researched and there are many much more advanced algorithms around.

The VLfeat website contains a nice more advanced demo of this approach, classifying the caltech 101 dataset.  Also noteworthy, are results and software from Caltech itself.


On what it means for a system to be minimum phase

The minimum phase system will have the minimum group delay for a given magnitude response. The phase response will have the least excursion over the frequency domain (due to all zeros being inside the unit circle). Since group delay is the negative derivative of phase with respect to frequency, a minimum phase system also has the minimum group delay.
The minimum phase excursion for a minimum phase system is well demonstrated with the graphic below, knowing that the complex frequency response (magnitude and phase) of an FIR filter can be determined graphically through phasors rotating at a rate according to the coefficient number (the first coefficient rotates zero times, the second coefficient rotates once, etc as the frequency axis is swept from $0$ to $2\pi$). Both filters shown below will have an identical magnitude response, but the one with the larger coefficient appearing first is minimum phase while the one with the coefficients reversed (the reverse filter) is maximum phase. We can see with the minimum phase system how the total phase deviation is minimized, while in the maximum phase case the phase deviation will be maximized, but in both cases having the same magnitude.
The transfer function of the first filter is $1-0.5z^{-1}$ or $(z-0.5)/z$ with all zeroes inside the unit circle as one zero at $z=0.5$ (minimum phase!). While the transfer function of the second filter is $0.5-z^{-1}$ or $(z-2)/z$ with all zeroes outside the unit circle with one zero at $z=2$ (maximum phase!). Notice too how in the minimum phase case, with the dominant tap placed earlier, how most of the signal will be at the output of the filter sooner (least delay).
When all the zeros are inside the unit circle, the frequency response as equivalently mapped out as in these graphics by rotating phasors will never enclose the origin, thus the phase excursion never exceeds $2\pi$.

For further details and how this relates to Cauchy's Argument Principle please also see: Can I set a constraint on the first tap of an FIR filter such that its inverse is stable?


A system with simple distinct poles on the imaginary axis (and note that the origin is on the imaginary axis) and no poles in the right half-plane is called marginally stable. If you have poles with multiplicity greater than $1$ on the imaginary axis, or if there are poles in the right half-plane, then the system is unstable.

For discrete-time systems, the same is true if you replace "imaginary axis" by "unit circle".

Transients in marginally stable systems do not decay, but neither do they grow without bounds. In practice, we usually want to avoid marginally stable systems. An ideal integrator is an example of a marginally stable system.


Frequency shift property of Fourier transform

Your Eq. $(2)$ is the frequency shift property; Eq. $(1)$ is wrong. This is easy to show:

$$\begin{align}\mathcal{F}\left\{e^{j\omega_0t}x(t)\right\}=&\int_{-\infty}^{\infty}x(t)e^{j\omega_0t}e^{-j\omega t}dt\\&=\int_{-\infty}^{\infty}x(t)e^{-j(\omega-\omega_0)t}dt\\&=X(j(\omega-\omega_0))
\end{align}$$

Your 'Case II' has nothing to do with the frequency shift property, because it is no frequency shift. For a frequency shift you need to multiply by a complex exponential $e^{j\omega_0t}$, and not by a real-valued exponential $e^{at}$.

You can come up with a rule for the Fourier transform of $x(t)e^{-at}$, $a\in\mathbb{R}$, but this is more tricky than the frequency shift property. For the frequency shift property, if you know that $X(j\omega)$ exists, then you know for sure that also the Fourier transform of $x(t)e^{j\omega_0t}$ exists. On the other hand, if $X(j\omega)$ exists, it is not certain that also the Fourier transform of $x(t)e^{-at}$ exists. This has everything to do with the Laplace transform and its region of convergence:

$$\begin{align}\mathcal{F}\left\{x(t)e^{-at}\right\}&=\int_{-\infty}^{\infty}x(t)e^{-at}e^{-j\omega t}dt\\&=\int_{-\infty}^{\infty}x(t)e^{-(a+j\omega)t}dt\stackrel{?}{=}X(a+j\omega)\end{align}$$

The last equality only holds if the integral converges, i.e. if $s=a$ is inside the region of convergence of the Laplace transform of $x(t)$, and if the Laplace transform of $x(t)$ has no singularities on the imaginary axis, i.e. if its Fourier transform has no Dirac delta impulses. The latter explains why the above 'rule' doesn't work for $x(t)=u(t)$ (as in your example).


Resources for an introduction to computer vision

If (from the comments) you are interested in an introduction to computer vision (for the difference, see Question on difference between CV and IP), you could take a look at the introductory course of Coursera on Computer Vision..

I am also currently enrolled, but have not had much time to look into it, but it seems to me as a decent course (with the knowledge of computer vision I already have) to get an overview of the different subjects in computer vision, and the techniques that are used to overcome problems.

Next to this, there is a free download of the book Computer Vision by Richard Szeliski, which is also a great source for the basic principles of computer vision, accompanied with a lot of examples.

With this basic understanding and feeling with the subjects of computer vision, you will most likely understand more complicated things easier, because you can place things in perspective.


The QRS complex is a "spike" that lasts about 100 ms. There is also the P wave that comes before Q which indicates atrial contraction, and the ST segment after. The QRS complex is associated with the ventricles contracting. Blood is being pushed into the major arteries. This part of the cycle is called systole. QRS is like a spark plug if you would permit a simplistic analogy. It represents timing not blood flow. The waveform that represents blood flow is the instantaneous pressures in each  artery and vein and those waveforms are much smoother. The pulse on your wrist is a pressure waveform.

See here for a good tutorial.

While an average heart rate is about one a second, there is a lot of beat to beat variability even when an individual is resting. It's not really periodic in the Fourier decomposition sense.   Signals that are not periodic have a continuous frequency spectrum and the majority of the energy for QRS is in the 40 to 100 Hz band. Blood pressure waveforms have most of their energy at the lower frequencies.


In general, the symbol error probability is given by

$$P_E=\sum_{i=1}^MP[s_i]P[E|s_i]\tag{1}$$

where $P[s_i]$ is the probability that symbol $s_i$ was sent, $P[E|s_i]$ is the error probability given that $s_i$ was sent, and $M$ is the number of symbols. If all symbols are equally likely, $(1)$ simplifies to

$$P_E=\frac{1}{M}\sum_{i=1}^MP[E|s_i]\tag{2}$$

For $16$-QAM it is possible to derive an exact expression for the error probabilities $P[E|s_i]$ (and hence for the symbol error probability $P_E$). You can find the derivation in most digital communication text books. A good approximation (better than the union bound) is given in this answer.

The union bound approximation gives an upper bound on $P[E|s_i]$:

$$P[E|s_i]\le\sum_{j\neq i}P[s_j\text{ chosen}|s_i\text{ sent}]\tag{3}$$

Under the given assumptions we have

$$P[s_j\text{ chosen}|s_i\text{ sent}]=Q\left(\sqrt{\frac{d^2_{ij}}{2N_0}}\right)\tag{4}$$

where $d_{ij}$ is the distance between symbols $s_i$ and $s_j$. With $(3)$ and $(4)$ the union bound approximation of $(2)$ becomes

$$P_E\le\frac{1}{M}\sum_{i=1}^M\sum_{j\neq i}Q\left(\sqrt{\frac{d^2_{ij}}{2N_0}}\right)\tag{5}$$

Note that generally you will get a sum of $Q$ functions with different arguments because the distances $d_{ij}$ are usually not all the same. This is also the case for $16$-QAM.

The solution you provided is still more conservative (hence, worse) than the union bound. That solution uses the fact that the $Q$ function is monotonically decreasing, which implies

$$Q\left(\sqrt{\frac{d^2_{ij}}{2N_0}}\right)\le Q\left(\sqrt{\frac{d^2_{min}}{2N_0}}\right)\tag{6}$$

with $d_{min}=\min_{i\neq j}\{d_{ij}\}$ being the smallest occurring distance between symbols. Plugging $(6)$ into $(5)$ we finally get the bound

$$P_E\le (M-1)Q\left(\sqrt{\frac{d^2_{min}}{2N_0}}\right)\tag{7}$$

which, as mentioned above, is even less tight than the union bound.


General Questions MIMO Communication theory

As it says in the documentation, SampleRate is the signal's sampling rate, not the symbol rate.
The path delays together with your symbol rate determine whether the channel is flat or frequency-selective.
You should choose Rayleigh if there is no line-of-sight component in the received signal, and Rice otherwise. Regarding the path delays: the model is that the gain you observe at each sampling instant is the result of many paths, whose arrival times are so close that you can't resolve them. So, in your model, you always have a finite number of paths, each of which has random Gaussian gain.
Each path has a different length and arrival angle, so it can have a different Doppler.

As a side note, and expressing only my opinion, it is more educational to try to create your own models of the fading channel than using a Matlab toolbox. Toolboxes are created to increase an expert's productivity, not for learning.


problems calculating phase propagation (in C)

you don't wanna do this on a bin-by-bin basis.  here's a paper by Miller Puckette that is seminal regarding this issue.  what you want to do is identify specific frequency components (perhaps each lobe in the spectrum would be a specific frequency component) and multiply that entire $n$th lobe by $e^{j \phi_n}$  where $\phi_n = \omega_n \ T_n$  and  $\omega_n$ is the frequency of the $n$th component and $T_n$ is the time adjustment to that component to match its phase to the corresponding component in the previous frame.

each bin gets multiplied by $e^{j \phi_n}$, but for most bins that factor is the same as for an adjacent bin.


Why does the excitation signal appear, separated, at high quefrencies in the cepstrum?

When using an FFT, an evenly spaced sequence of events in one domain usually produces a strong component in the other domain at a location related to the spacing of the events in the first domain.

A voiced speech signal usually includes a lot of harmonics which are evenly spaced in the frequency domain.  These evenly spaced events in the frequency domain will thus produce a component in the quefrency domain whose location is related to the spacing of the harmonics.  And the spacing of the harmonics in the frequency domain is related to the spacing of the exciter events in the time domain (the period of the fundamental).


What is the cut-off frequency for this low pass filter?

Look for Sallen and Key filter.
Fc = 1/(2 * pi * sqrt(R1 * R2 * C1 * C2))
Or use a calculator, found one here: Sallen-Key Low-pass Filter Design Tool
(It appears to handle "k" and "n", as a bonus)
Fc = 232 Hz, if you want to check your math.


In order to block the DC component in a signal, the high-pass filter can be used.

In the figure above, a basic high-pass filter which is known as the passive RC high-pass filter can be seen. The working principle of it depends heavily on the capacitor’s behaviour at varying frequencies as the formula depicts:
$$X_{C}(j\omega) = \frac{1}{j\omega{C}} \Omega$$
The formula is named as the capacitive reactance which is actually a capacitor’s resistance (not the internal resistance) by which an alternating current is attenuated while it passes through it.
So, let’s analyse the given circuit. If the Kirchhoff’s voltage law is applied on that circuit;
$$-x(t) + {\frac{1}{C}\int_{t_{0}}^{t} i(t) \, dt} + V_{C}(t_{0}) + i(t)R = 0$$
$$-x(t) + {\frac{1}{C}\int_{0}^{t} i(t) \, dt} + V_{C}(0) + i(t)R = 0$$
Let’s assume that the capacitor’s inital voltage $V_{C}(0) = 0V$.
$$-x(t) + {\frac{1}{C}\int_{0}^{t} i(t) \, dt} + i(t)R = 0$$
$${\frac{1}{C}\int_{0}^{t} i(t) \, dt} + i(t)R = x(t)$$
If the last equation is differentiated with respect to time and the relation $y(t) = i(t)R$ is used according to the circuit;
$${\frac{1}{C}i(t)} + R\frac{di(t)}{dt} = \frac{dx(t)}{dt}$$
$${\frac{1}{RC}y(t)} + \frac{dy(t)}{dt} = \frac{dx(t)}{dt}$$
By rearranging the last equation, we obtain;
$${\frac{dy(t)}{dt} + \frac{1}{RC}y(t)}  = \frac{dx(t)}{dt}$$
Interpreting what is happening in a high-pass filter, low-pass filter, etc. by looking at their differential equations can be demanding. Since these circuits are used to make changes on signals with regards to changes in frequencies, analysing these circuits in the frequency domain can be quite useful. In order to do this, the Fourier transforms of the both sides of the last differential equation can be taken.
$${j\omega}Y(j\omega) + {\frac{1}{RC}}Y(j\omega) = {j\omega}X(j\omega)$$
$${Y(j\omega)}(j\omega + \frac{1}{RC}) = {j\omega}X(j\omega)$$
Consequently, the transfer function $H(j\omega)$ becomes;
$$H(j\omega) = \frac{Y(j\omega)}{X(j\omega)} = \frac{j\omega}{j\omega + \frac{1}{RC}}$$
As it can be seen from the formula, basically, the transfer function results in a complex number. In order to obtain its magnitude, the absolute value of it can be calculated.
$$|H(j\omega)| = |\frac{Y(j\omega)}{X(j\omega)}| = |\frac{j\omega}{j\omega + \frac{1}{RC}}|$$
$$|H(j\omega)| = \frac{|j\omega|}{|j\omega + \frac{1}{RC}|}$$
$$|H(j\omega)| = \frac{\omega} {\sqrt{ {\omega^{2}} + \frac{1}{R^{2} C^{2}} }}$$
This equation relates the input voltage of the passive RC high-pass filter to the output voltage by the input voltage’s frequency. Actually, this equation explains the main working principle of the passive RC high-pass filter by itself.
In order to make the transfer function’s behaviour i.e. the high-pass filter’s behaviour over frequency perceptible, let’s create an example circuit with defined values of $R$ and $C$ and illustrate the magnitude graph of the transfer function. I’ve chosen the values arbitrarily and assumed that $R = 20 k\Omega$ and $C = 1 \mu{F}$.
The MATLAB code for drawing the magnitude graph of the transfer function is given below:
f = 0: 0.01: 20;  % Frequency range in Hz

R = 20000;  % Resistance value in Ohms
C = 0.000001;  % Capacitance value in Farads

omega = 2 * pi * f;  % Angular frequency in rad/s

mtf = omega ./ sqrt((1 / power(R * C, 2)) + power(omega, 2));

plot(f, mtf, "r", "LineWidth", 2);
grid on;
xlabel("Frequency (Hz)");
ylabel("|H(f)|");
title("Frequency Response of the Passive RC High-Pass Filter");

I’ve plotted the graph with respect to the frequency as most of the time, we are dealing with classic frequency rather than angular frequency.
The resultant plot is given below:

As it can be seen from the plot above, the corner/cut-off frequency of the example high-pass filter is around 8 Hz which can be calculated by the formula;
$$f_{C} = \frac{1}{2\pi{RC}} Hz$$
Corner/cut-off frequency is the frequency at which the output signal’s power which is transferred to the load (in our case, it is the resistor $R$) is half of that of the input signal. Alternatively, it is the point at which the output voltage’s magnitude is approximately 0.707 times the input voltages’s magnitude.
So, for our passive RC high-pass filter; if the input signal’s frequency is chosen lower than the corner/cut-off frequency of the filter, the input signal will be attenuated dramatically.
Since a DC signal has the frequency of $f_{DC} = 0 Hz$ ideally, it can be removed from the input signal.
However, I don’t have any idea about why the manufacturer has chosen the cut-off frequency as 4 Hz but choosing it lower than the minimum frequency value of the human hearing range which is 20 Hz is obvious as we don’t want to miss a frequency component while we are listening to music, conversations, etc. via a loudspeaker.
This specific cut-off frequency rating could be about the designer’s state of mind or the circuitry’s needs.


Assuming 0-based vectors, hop size $H$, and window size $M$, the frame number $N_f$ will contain the DFT of a window that starts on sample $n_1 = N_f \times H$ and ends on sample $n_2 = n_1 + M -1$.

That means that for each STFT frame, you shift the window by $H$ samples. Note that $H$ and $M$ not necessarily coincide.


What are good ways to detect signal clipping in a recording?

I was in the middle of typing an answer pretty much exactly like Yoda's.  He's is probably the most reliable but, I'll proposed a different solution so you have some options.

If you take a histogram of your signal, you will more than likely a bell or triangle like shape depending on the signal type. Clean signals will tend to follow this pattern. Many recording studios add a "loudness" effect that causes a little bump near the top, but it is still somewhat smooth looking. Here is an example from a real song from a major musician:

Here is the histogram of signal that Yoda gives in his answer:

And now the case of their being clipping:

This method can be fooled at times, but it is at least something to throw in your tool bag for situations that the FFT method doesn't seem to be working for you or is too many computations for your environment.


Shifting signal smaller than discrete step

I don't do much image processing but maybe the simpler 1D case will help:

A delay in the time domain is a phase shift in the frequency domain, this notion can be used to shift your FFT in the frequency domain by a rotation much smaller than a sample

$$x(n-n_o) = W_N^{kn_o}X(k)$$

$$W_N^{kn_o} = e^{\frac{-j 2 \pi k n_o}{N}}$$

You need to apply this principle row wise and column wise. I tried giving it a go by building a 2D matrix of the x y shift and then applying it to the image. For integer shift this works, when you get into fractional shift (0-.99) for $n_o$ it gets weird.

Non-fractional shift of a "cat eating fancy ice cream" to prove that it does what I think it does:

Matlab code:

close all
clear all

%load image
I = imread('cat_eating_fancy_ice_cream.jpg');

I = im2double(I);

%make phase shift matrix for frequency domain shifting
sz = size(I)
y = e.^((-j*2*pi*[0:sz(1)-1]*100)/sz(1)).';
x = e.^((-j*2*pi*[0:sz(2)-1]*200)/sz(2));
shift = y*x;

%convert image to grayscale for easier manipulation
pF1 =(fft2(rgb2gray(I)));

%phase shift frequency domain info
pF = pF1.*shift;

%plots
figure
subplot(121)
imshow(rgb2gray(I))
subplot(122)
imshow(ifft2(pF))


Can anyone suggests some basic readings for video processing?

I recommend Video Tracking: Theory and Practice  as comprehensive overview of what have been done in this area for last 15 years. Also, 3D Computer Vision: Efficient Methods and Applications for quick understanding useful algorithms.

This books will be a good start, but when you have enough knowledge, you should start looking into research papers, because information there is more updated than in books.


motion compensation for range migration algorithm

I have been trying to incorporate MoCo to the extended Omega-k as well following Reigbar's paper (https://ieeexplore.ieee.org/document/1648550). I understand that for the first order MoCo, the range deviation needs to be calculated from a reference range (i.e. image center) and for the second order, the differences between the deviation for all ranges and the deviation for the reference range need to be used. However, calculating the deviation in the above recommended method didn't work for me. Especially, I'm not sure if I calculated the deviation for the second order MoCo right. I used r = (c*fast_time)/2 for calculating all the range positions and then subtracted those from the non-ideal platform positions in order to know their range deviations. Is that correct?

Thank you.


Sequences that are "time-invariant" must be constants (not varying with time at all).
What the OP is actually asking about is a (discrete-time) system whose response $y$ to the input $x$ is defined to be
$$y[n] = \begin{cases}\sum_{k=n_0}^n x[k], & n \geq n_0,\\??? & n < n_0\end{cases}$$
where the $???$ is there because according to some traditionalists, if $n < n_0$, then that sum is an empty sum which has value $0$ but more liberal-minded might treat that sum as representing $$x[n] + x[n+1] + x[n+2] + \cdots + x[n_0-1] + x[n_0]$$ (with more serious-minded folks even negating that sum in analogy with the standard result that for $a<b$, $\int_b^a = -\int_a^b$). So, for the traditional meaning of things and $n_0$ denoting a (finite) integer, the system is not a time-invariant system. It is time-invariant (with unit pulse response $h[n]$ being the unit step function $u[n]$) if $n_0 = -\infty$ as noted in a comment by Robert Bristow-Johnson on the main question.


FIR filters: direct form / transposed FIR

A finite impulse response (FIR) digital filter implements the following convolution sum

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

for each output sample $y(n)$, where $x(n)$ is the discrete-time input signal, $h(n)$ is the filter's impulse response, and $N$ is the filter length. The values $h(n)$ are also called filter taps, and $N$ is then referred to as the number of taps. The filter described by Equation (1) is also called an $N$-tap filter.

Direct-form and transposed direct-form are just different implementations, i.e. different ways to compute the sum in (1). In theory they are identical, but when computed with finite precision, there can be differences between the different implementations. The direct-form FIR structure is also called tapped delay line or transversal filter.

The two realizations below are the direct-form structure (transversal filter, tapped delay-line) and the transposed structure (from Oppenheim and Schafer, Discrete-time Signal Processing):


Duality in the discrete-time Fourier series

Based on the references and the definitions given, here is a proof.
Given:
$$x[n] \leftrightarrow a[k], \ and$$
$$y[n] \leftrightarrow b[k]$$
Fourier Series expansion of $x[n]e^{j\frac{2\pi}{N}nm}$ would be derived as follows:
$$\frac{1}{N}\sum^{N-1}_{n=0}x[n]e^{j\frac{2\pi}{N}nm}e^{-j\frac{2\pi}{N}nk} = \frac{1}{N}\sum^{N-1}_{n=0}x[n]e^{-j\frac{2\pi}{N}(k-m)n} = a[k-m]$$
The factor $\frac{1}{N}$ gets consumed in the definition of Fourier Series coefficients. Therefore, the following relationship can be stated:
$$x[n]e^{j\frac{2\pi}{N}nm} \leftrightarrow a[k-m]$$
Similarly, for the second case of multiplication in time-domain, you need to figure out that the inverse Fourier Series relation or the synthesis expression will be the following:
$$x[n] = \sum^{N-1}_{k=0}a[k]e^{j\frac{2\pi}{N}kn}$$
I am not adding the proof of the above which can be easily derived by plugging the expression for $a[k] = \frac{1}{N}\sum^{N-1}_{n=0}x[n]e^{-j\frac{2\pi}{N}kn}$ on the RHS.
Then all we are required to do is proving that the inverse Fourier Transform of convolution of $a[k]$ and $b[k]$ would be product $x[n]y[n]$. Which is as follows:
$$\mathcal{FS}^{-1}(a*b)[k] = \sum^{N-1}_{k=0}(a*b)[k]e^{j\frac{2\pi}{N}nk} = \sum^{N-1}_{k=0}( \sum^{N-1}_{l=0}a[l]b[k-l]) e^{j\frac{2\pi}{N}nk}$$
$$=>\mathcal{FS}^{-1}(a*b)[k] =  (\sum^{N-1}_{l=0}a[l]e^{j\frac{2\pi}{N}ln})(\sum^{N-1}_{k=0}b[k-l]e^{j\frac{2\pi}{N}(k-l)n})$$
$$=> \mathcal{FS}^{-1}(a*b)[k]  = \mathcal{FS}^{-1}a[k] \mathcal{FS}^{-1}b[k] = x[n]y[n]$$
I have replaced $m = (k-l)$ in the proof above to get $y[n]$.
Hence both the equations (5.69) and (5.71) are true and correctly given in the book.
I think you are trying to use equation (5.67) to prove equation (5.69) and (5.71). That is not correct, because (5.67) is stating Duality of Fourier Series expression between Continuous Time and Discrete Time. Equation (5.69) and equation (5.71) are duality of Shift-Property and Convolution Property. The second reason you cannot use Duality Expression (5.67) is because (5.67) states that if $\mathcal{FS}g[n] = f[k]$, then $\mathcal{FS}f[n] = \frac{1}{N}g[−k]$. That is it takes Fourier Series of the Fourier Series itself assuming $f[n]$ in time-domain. But in (5.69) and (5.71), the time-domain sequences remain $x[n]$ and $y[n]$ only, so (5.67) does not apply here. There are 3 Duality Properties stated in equation (5.67), (5.69) and (5.71). They are not the consequence of one another. They are a consequence of similarity in definition of Continuous time Fourier Series and Discrete Time Fourier Series.


Differential Equations applications in Signal Processing

I don’t think it will hurt or will be a waste of time.

You can think of SP in a number of ways and those ways overlap. It’s a technolgy, a science,a branch of mathematics, an engineering discipline, and your view is largely tied to your self image.

A lot of SP is based on physical models.  A lot is based on hueristics.

SP overlaps. You might benefit from taking music or statistics. A lot of people take the longer path.

Differential equations are a way to describe a lot of phenomena that SP is used to sense and often control.

You should consider what motivates you.

To be honest a lot of what I took had to do with what was offered that semester. Getting that piece of paper on time was a huge practical consideration.


Estimate sampling frequency given an array of values

Given an array x of size 2000, is there any way to calculate (or estimate) the sampling frequency [...]? The values are the only information I have about the signal.

No.

  I know the signal has been sampled with Fs = 360 Hz but I want to know if there is any way to extract this information from these plots.

No.

The sampling frequency is an intrinsic characteristic of a discrete sampled signal. Any operation on the sampled signal uses Fs either implicitly or explicitly. The only way to infer the sampling frequency is to inject or detect a pilot tone on a properly sampled signal.

For example, the Discrete Fourier Transform does not make explicit reference of the sampling frequency. The DFT simply decomposes the signal into the sum of those sinusoids that can be supported by the "data frame" of length $N$ that it is presented with. But, when you then want to go from "wave number" or "bin", back to the physical frequency it represents, then you have to know the Fs.

If you don't use the "right" Fs, then you might run into the condition of aliasing where you basically cannot trust the numbers you are looking at because they are an uncontrolled overlap of the signal spectrum and itself. This is why, the anti-aliasing filter is a key component in the analog pathway of the signal BEFORE the analog to digital conversion. The anti-aliasing filter will make sure that no frequency very near or above $\frac{Fs}{2}$ will make it through the system.

If we were to remove this filter, the digital processing part's outputs would become unreliable.

Having said this and as is pointed out by @Fat32, if you are absolutely sure that the sampled signal adheres to this restriction and it was sampled without aliasing then you can infer the Fs. BUT, this cannot be done with a random Electrocardiogram. Or rather, it really depends on the sort of accuracy you are after.

A standard electrocardiogram waveform's harmonics can extend to 100Hz.

The question now is, what is the actual bandwidth of THIS specific ECG that you have in your hands? That one, may indeed have been obtained with a standard montage but it contains noise and baseline wandering. Also, depending on what the patient was doing, the signal would be different and its bandwidth would be different too.

Therefore, which DFT bin corresponds to that 100Hz is uncertain.

If you are interested in a ball park figure, you can look at a point in the signal's DFT where its power drops below an average background noise and take that point as your 100Hz nominal value and work backwards to what Fs was used. But that would certainly not be accurate, especially when "strange" Fs are used (like 128.5 Hz).

So, in the general case, it is impossible to "second guess" the sampling frequency from a discrete signal's values.

The only exception to this is to ensure that no aliasing will occur (i.e. sample high enough) and use a pilot tone. That is, a tone or other signal whose spectrum you know exactly or can control fully.

If there was a single pilot tone in the signal then you could say, I know that the pilot tone is at 1kHz, it shows up on the DFT output as a high output on bin k, then I can work backwards to estimate the correspondence of bins to physical frequencies.

Could there be such a pilot tone in your signals?

Yes, it is likely that your signals contain power line interference. This will show up as a consistent set of spikes in the DFT because it is not only that you would be looking at a sinusoid of 50Hz but you would get a few harmonics along too.

The frequency of the power line is not exactly 50Hz (in Europe at least) but it has to adhere to strict standards, so perhaps for the relatively short duration of the ECG it could be used as a pilot tone at an assumed value of 50Hz.

Of course, it is not necessary that the ECG contains power-line noise. It may have been suppressed by a filter or other shielding at the point of acquisition. In that case, you are back into "it cannot be done" response.

Hope this helps.


The problem is not sufficiently specified, because the range of admissible values of $n$ is missing. Here I make the assumption that we consider $n>0$. With this assumption we have

$$X(z)=\sum_{n=1}^{\infty}x[n]z^{-n}=\sum_{n=1}^{\infty}\frac{z^{-n}}{n^2}\tag{1}$$

And that's the point where we might get stuck, if we didn't have a list of mathematical series, or if we didn't know about the polylogarithm, which is defined by

$$\text{Li}_{s}(z)=\sum_{n=1}^{\infty}\frac{z^n}{n^s},\qquad |z|<1\tag{2}$$

where $s$ is an arbitrary complex number.

In your case, $s=2$ and the corresponding function is called the dilogarithm or Spence's function.

Comparing $(1)$ and $(2)$ we get for the $\mathcal{Z}$-transform of $1/n^2$ for $n>0$

$$X(z)=\text{Li}_2\left(\frac{1}{z}\right),\qquad |z|>1\tag{3}$$

Another way to arrive at the solution is to use the differentiation property of the $\mathcal{Z}$-transform:

$$nx[n]\Longleftrightarrow -z\frac{dX(z)}{dz}\tag{4}$$

Applying $(4)$ twice will give you the result. You need the correspondence

$$u[n-1]\Longleftrightarrow \frac{1}{z-1},\qquad |z|>1\tag{5}$$

In a first step you'll arrive at the transform of $1/n$, $n>0$, and in a second step you'll arrive at the transform of $1/n^2$, $n>0$.

In this case you will be using the integral representation of the dilogarithm:

$$\text{Li}_2(z)=-\int_0^z\frac{\ln(1-u)}{u}du\tag{6}$$


So you have
$$
\int_{-\infty}^{+\infty} \delta(u) \delta(t-u) du
$$

Try using the sifting property of $\delta(t)$:

$$
\int_{-\infty}^{+\infty} f(u) \delta(T-u) du = f(T)
$$


First of al it's a very exciting idea to capture a room's impulse response. Hence congratulations on your challenges. Assuming that you have correctly implemented the software he provided in his page, then the only cause of the differences are; (Note: I can say nothing about his method of computing the RIR. It can be a proper method or not... Only the possible cause of differences you observe I'll discuss)

1- Speakers: He uses monitor quality Yamaha speakers, you use crappy PC speakers (just like mine or yours; they are all crappy when compared to a Yamaha monitor speaker)

2- Microphones: He uses high fidelity, high linearity, low noise, high sensitivity microphones , you use web-cam crappy microphone. (again we all use those PC microphnes but they are crappy compared to any professional one)

3- Speaker and microphone placement: it's very vital that you place the speakers and the microphones to the correct places, otherwise you will be measure something that you don't want or need to.

4- Room: His impulse response is from a church. Your's is from a living room. When compared to such a structure, so you should be measuring different room impulse reponses of course.

However, even the most interesting part is, how and why do you judge your impulse response from its graphical plot appearance? You should judge it from its acoustical result point of view.

Have you tried your calculated room impulse response to add a reverb effect to dry audio, to see if it works? i.e.,to see if the dry audio becomes as if it was played and recorded in that room you refer to its impulse response?

You should do this test with a head-set ! otherwise you would be applying the room impulse response twice when played back in the same room...


Neither of the two equations match with the diagram. It's actually quite straightforward: at the input of the delay you have

$$x[n]+\frac12\big(y[n]+y[n-1]\big)\tag{1}$$

Then a delay of $N$ samples is applied to compute the output, hence

$$y[n]=x[n-N]+\frac12\big(y[n-N]+y[n-1-N]\big)\tag{2}$$

But here I've found a different diagram, which corresponds to the last equation in your question.


Removing lines of certain length from an image in MATLAB

This is a program which removes the line segments in a single column of length which is more than the threshold

function [ret_im] = remove_vertical(in_image)

[NZ_r, NZ_c] = find(in_image);

ori = in_image;

%scan through each row of location of non-zero pixels
% If three (threshold) or more pixels are placed in a single column, turn them to zero
thresh = 3;
for i=1:1:size(NZ_r,1)-1

    probe = NZ_c(i); % Set the initial search element
    len = length(find(NZ_c(:)==probe));
    if( len > thresh)
        for j=1:1:len
            if((NZ_r(i+1) - NZ_r(i)) == 1)
                in_image(NZ_r(i),NZ_c(i)) = 0;
            end
        end
    end
end

ret_im = in_image;

end


Energy and power of the product of an exponential and sinusoidal signal

Let's write the signal $x(t)$ in the general form

$$x(t)=e^{-\alpha t}\cos(\omega_0t)u(t)\tag{1}$$

Note that its energy is given by

$$E_x=\int_{-\infty}^{\infty}|x(t)|^2dt=\int_{0}^{\infty}e^{-2\alpha t}\cos^2(\omega_0)tdt\tag{2}$$

Since this is a homework type problem I won't solve it for you but I'll give you some hints:

Expand the integrand of $(2)$ as
$$e^{-2\alpha t}\frac12(1+\cos(2\omega_0t))=\frac12e^{-2\alpha t}+\frac12e^{-2\alpha t}\frac12(e^{j2\omega_0t}+e^{-j2\omega_0t})$$
Now you have two integrals over basic exponential functions, which should be very straightforward to solve.
The final result should be
$$E_x=\frac{1}{4\alpha}+\frac{\alpha}{4(\alpha^2+\omega_0^2)}$$
If the energy $E_x$ is finite, what can you conclude about the power of $x(t)$?


DDS generated chirp and antialiasing filter : what is the best suited filter for a chirp?

I recommend a digital FIR filter to both compensate for the DAC Sinc droop and added group delay distortion combined with a simpler analog filter focused on rejecting the higher images out of the DAC with the minimum number of components (this would eliminate using a Bessel filter, while a Butterworth or elliptic filter would be viable candidates).  The analog filter cost and complexity can be minimized by utilizing digital pre-distortion to equalize the analog filters group delay variation with a reasonable tolerance.
For further details on the design details for the anti-alias filter, please see this post: Where should I set my anti-aliasing filter corner frequency for this signal?
This post may also be of interest as an implementation approach where I detail all the code for an optimized chirp for a flat FFT response which would also have minimum aliasing effects under the same sampling conditions. Instead of a time domain filter, the response is scaled with a Tukey window which provides a constant envelope chirp over most of its duration, while scaling the very start and end of the chip minimize aliasing effects. This could be done with a weighting factor of the amplitude of the chirp at the output of the DDS. Please see the very bottom of this post for more details showing a flat FFT response over most of the chirp duration: How can I plot the frequency response on a bode diagram with Fast Fourier Transform?


Both are reasonable approaches and it is foreseeable that either one could outperform the other empirically. The Euclidean distance assumes the data to be isotropically Gaussian, i.e. it will treat each feature equally. On the other hand, the Mahalanobis distance seeks to measure the correlation between variables and relaxes the assumption of the Euclidean distance, assuming instead an anisotropic Gaussian distribution.

If you know a priori that there is some kind of correlation between your features, then I would suggest using a Mahalanobis distance over Euclidean. Also, note that Z-score feature scaling can mitigate the usefulness of choosing a Mahalanobis distance over Euclidean (less true of min-max normalization though). The major drawback of the Mahalanobis distance is that it requires the inversion of the covariance matrix which can be computationally restrictive depending on the problem.


How can I describe a sum of sinusoids by an equivalent lowpass signal?

The sum of two sinosoids with different frequencies $f_{c1}, f_{c2}$ does not add up to a single sinosoid. Only, if both sinosoids have the same frequency, you can combine them into a single sinosoid using the equation you state at the end of your post. See also for example this post.


How to determine SNR required for adding AWGN noise to clean IQ samples?

Note that the MATLAB result is floating point (or may be fixed point but with much higher precision) while the VSA hardware is fixed point with lower precision that also has dynamic range limitations in the analog portion of the hardware as well. To see this directly and more clearly, try these two experiments:

Create a single carrier tone and capture results at different signal levels of that carrier. Ensure the sampling rate is much higher than the tone frequency (such as > 12x higher to capture the 3rd and 5th harmonic).

Observe the output with a properly windowed FFT (I recommend the Kaiser window) and increase the capture duration until the higher harmonics are clearly visible. (increasing the duration will lower the FFT noise floor).

From such testing you can determine the SNR (based on the wideband noise floor) as well as spurious free dynamic range (SFDR) of your test set-up. The SNR is estimated as the difference between the tone and the integrated power that is not at the tone frequency (and outside the main lobe of the window used).  The SFDR is the difference in power between the tone and the highest spur which may or may not be at a harmonic of the signal (In some definitions of SFDR the harmonic components are not considered and that distortion is captured in Harmonic Distortion). In any event find the maximum tone power that maximizes SFDR, and this sets a mark for "Full Scale" from which the other metrics can be computed.

If you proceed with this, please post your results. What would be interesting too is a plot of Pin vs Pout for your tone, from which you could also determine where it starts to compress. The point where Pout is 1 dB lower than the otherwise linear relationship is known as the 1 dB compression point.

To be complete you could also do two tone testing, where you generate two sinusoids closely spaced (but wider than the width of the main lobe of your window in the frequency domain) and measure the two-tone third order intercept point. See Wikipedia for more details on that test:

https://en.wikipedia.org/wiki/Third-order_intercept_point

I would want to understand all these parameters in conjunction with the use of any test equipment such that I understand what is being introduced by the test equipment rather than in the waveform itself.

Finally when adding AWGN to the signal to simulate the SNR for a discrete system, the AWGN is a power spectral density that is flat across the entire sampling bandwidth. The SNR of interest would be the portion of the spectrum that is in-band---thus if you are sampling at a higher rate, the AWGN with the same variance would be spread across a wider bandwidth resulting in lower SNR given the portion of the noise that would occupy the same bandwidth of the signal. Any noise outside of the signal bandwidth can be filtered so is not part of the SNR computation.


Can I run a GA to optimize wavelet transform?

You should be able to use GA or any other derivative-free optimization methods to solve your problem [1]. The correct parameter values can be determined using your application. If I'm trying to optimize for compression then I will look into the entropy of the resultant coefficients. Deciding on the actual values for your parameters will depend on how you plan to use them. The correct parameter values are simply the ones that work better.


There isn't any obvious effect in the frequency domain. However, extending the duration of the cyclic prefix does make your signal more resistant to multipath channels with a large delay spread. As a rule of thumb, you need to size the cyclic prefix duration to cover the entire delay spread. If this condition is met, then the linear convolution that is imposed on your signal by the multipath channel will manifest itself as circular convolution on the cyclic-prefixed signal. This circular convolution can be equalized by simple pointwise multiplication in the frequency domain at the receiver.

This is why one of OFDM's advantages is often cited as the fact that it permits the use of a "one-tap equalizer." In single-carrier systems, equalization can be one of the harder things to get right, so anything that can simplify that portion of the receiver is a help.


Yes, there is a way to combine radar sensor data.  There are instances where you can correlate the two signals with each other, but it depends on what you are trying to do.  A couple of points about your examples:

Chin/cheek.  At 25 MHz, any reasonable aperture (antenna) will likely have a large enough beamwidth that both sensors would see the same thing: both cheek and chin.  That said, you may see a Doppler shift from both objects.

Top.  Conjugate one of them, provided that all other parameters are exactly the same.  That said, one sensor is already measuring the exact event that you describe.  If you need higher SNR, there are other things you can do.

And don't forget that each radar can pollute each other's data


Confussion on roll-off factor, sampling time and Raised Cosine pulses.

An easy way to see what's going on is by plotting the signal's eye diagram.

With low roll-off factors, the eye diagram will show you that deviating a small amount from the optimum sampling time will decrease the error margin. Note that the signal peaks are not the ideal sampling points. This is the eye diagram of a sinc pulse (zero roll-off):

With larger roll-off factors, the eye remains quite open for a larger period of time around the optimum sampling time. This is why it's said that these pulses are more tolerant of timing errors. This is the eye diagram when using pulses with roll-off equal to 1:

Mind you, this doesn't mean that you can be sloppy about timing synchronization. It's just that small errors in timing will result in larger variation in the sample values of low roll-off pulses compared to large roll-off pulses.


Adequate representation of frequency domain amplitude/magnitude of FFT of a signal

The correct way to group multiple bins together is to multiply each complex fft bin output by its complex conjugate (which gives the bin power) then add all the bin powers together and divide by the number of bins in the group. If you want to display in db (which is the conventional approach) then take 10*log10() of the result. Negative db values are normal and expected.

Note that if your audio signal is music, it is quite normal that the spectrum falls off with frequency.

Also from your description it sounds like you are not using a 2 ^n number of points. It's best to stick with 2^n data buffer sizes to maximize the efficiency of the fft computation.

Bob


If you have audio that will be played by TV beforehand, then this problem becomes almost similar to echo cancellation, i.e. adaptive filters could be used and cancellation will be great.

What we have is the audio that will be played through TV, this audio will be used as reference for adaptive filter and input to the adaptive filter would be the capture(noisy music in earphone/headphone).
so adaptive filter will try to estimate the TV's audio in the earphone capture, then this estimate could be subtracted from earphone capture to get noise free(TV noise) music in earphone/headphone.


"Is Windowed FFT called smoothing?" No.
"Do I need to do anything else apart from ifft(X_3term(m)) to perform inverse of Xthree−term(m)" No.
Here is what is going on:
Suppose
$$ X[k] = \sum_{n=0}^{N-1} x[n] e^{-i \frac{2\pi}{N} nk } $$
This is the conventional definition of the DFT which is unnormalized.
Now do the following with it:
$$
\begin{aligned}
Y[k] &=  \alpha X[k] - \frac{\beta}{2}X[k-1]-\frac{\beta}{2}X[k+1] \\
&= \alpha \sum_{n=0}^{N-1} x[n] e^{-i \frac{2\pi}{N} nk }
- \frac{\beta}{2} \sum_{n=0}^{N-1} x[n] e^{-i \frac{2\pi}{N} n(k-1) }
- \frac{\beta}{2} \sum_{n=0}^{N-1} x[n] e^{-i \frac{2\pi}{N} n(k+1) } \\
&= \alpha \sum_{n=0}^{N-1} x[n] e^{-i \frac{2\pi}{N} nk }
- \frac{\beta}{2}  \sum_{n=0}^{N-1} e^{ i \frac{2\pi}{N}n} x[n] e^{-i \frac{2\pi}{N} nk }
- \frac{\beta}{2}  \sum_{n=0}^{N-1} e^{-i \frac{2\pi}{N}n} x[n] e^{-i \frac{2\pi}{N} nk } \\
&=  \sum_{n=0}^{N-1} \left[ \alpha - \frac{\beta}{2} e^{ i \frac{2\pi}{N}n} - \frac{\beta}{2} e^{-i \frac{2\pi}{N}n} \right] x[n] e^{-i \frac{2\pi}{N} nk } \\
&=  \sum_{n=0}^{N-1} w[n] x[n] e^{-i \frac{2\pi}{N} nk } \\
\end{aligned}
$$
Where
$$
\begin{aligned}
w[n]  &=\alpha - \frac{\beta}{2} e^{ i \frac{2\pi}{N}n} - \frac{\beta}{2} e^{-i \frac{2\pi}{N}n}  \\
&= \alpha - \beta \cos \left( \frac{2\pi}{N}n \right) \\
\end{aligned}
$$
So
$$ y[n] = \text{DFT}^{-1} ( Y[k] ) =  w[n] \cdot x[n] $$
You might also enjoy:

The Zeroing Sine Family of Window Functions


Shannon formula, why square the channel gain?

The Shannon channel capacity equation, as I remember is

$$\begin{align}
R &= W \log_2 \left( 1 + \frac{S}{N} \right) \\
\\
  &= W \log_2 \left( \frac{S+N}{N} \right) \\
\\
  &= W \big( \log_2(S+N) \, - \, \log_2(N) \big) \\
\\
  &= \frac{W}{6.02} \big( (S+N)_\text{dB} \, - \, (N)_\text{dB} \big) \\
\end{align}$$

where $S$ is the signal power and $N$ is the noise power (both over the same bandwidth $W$).

Noise power for a given bandwidth (thin enough that noise power density $N_0$ is considered constant is

$$ N = N_0 W $$

and the received signal power is the transmitted power, $P_t$, times the power gain of the channel (assuming both are constant over the bandwidth $W$).  The power gain of the channel is the square of the gain in amplitude because instantaneous power is the square of instantaneous amplitude.  So the received signal power is

$$ S = |H|^2 P_t $$

assuming all of the transmitted power is in the channel bandwidth $W$ and the channel gain, $H$, is constant over the bandwidth $W$.


How to find quantization noise in my simulation?

The quantization error is the quantization noise. It is two different terms for the same thing.


I'm not sure if I understand your question correctly, but I'll try to make a few things clear. First of all, note that windowing and frequency sampling are not iterative. On the other hand, the Remez exchange algorithm (used in the ParksMcClellan program) is iterative.

All methods for the design of linear phase FIR filters impose the linear phase constraint in the formulation of the design problem, i.e., the filters are always guaranteed to have a linear phase. The difference between the different methods is either the way the optimal solution is computed, or the design criterion, i.e., the way the error to be minimized is defined. Some common design criteria are the least squared error criterion (used in the Matlab function firls.m), or the maximum (Chebyshev) error criterion, where the maximum error is minimized (used in the Matlab function firpm.m).


Opencv calcHist and calcBackProject in plain english

See What is Back Projection in the openCV tutorials

An image histogram measures the distribution of colour (and brightness) of the pixels in an image.

If you take an image and identify a region of interest eg. a hand, and calculate the histogram of the pixels in that object.

Then take that histogram and a second image and essentially reverse the process - you pick the pixels in the second image that match the histogram from the first. It's this reverse process that gives it the name back-projection.

You then make the assumption that areas of the image in the second image that have the same colour distribution as an object in the first image are an image of the same (or similar) object.


Scipy resample, "fourier method" explanation

Suppose you have initially a real-valued sequence x of length N. The function is basically doing this:

To upsample, it transforms to the frequency domain and adds N/2 zeros at the end. Then it transforms back to the time-domain.
To downsample, it transforms to the frequency domain and deletes the second and third groups of N/4 elements (which correspond to the half with the highest frequency components). Then it transforms back.

You can try it on your own running step by step the source code for the resample() function.

The theory regarding these methods can be found in this answer.


How to Implement the moving average filter in time domain in MATLAB?

HINT (really an answer but not providing code so its a hint):
In your setup, a $M$-point moving average filter uses the average of the previous $M-1$ samples and the current sample to compute the current output. So for $n=5$, the output will be $\frac{1}{5}(x[1]+x[2]+x[3]+x[4]+x[5])$. In general, for $n=k$ the output of this $M$-point MA filter is:
$$ \frac{1}{M}\sum_{n=k-M+1}^kx[k]$$


Counting the number of objects intersected by the image border in MATLAB

In my opinion, dilation erosion are very basic tools. You have a very strong information base and quite a decent input image to make such decisions.

Here is my view:

Given a reasonable success you have shown to go from Fig 1 to 3, you can identify and segment individual signs.
Assuming that you have known signs a prior, you can apply a decent fast algorithms to do pattern matching. In case if the exact patterns are not known- you can just identify the outer shape of the pattern.
Based on the classification, you can always define the centroid of each matched pattern and its respective width and height. If the centroid X,Y position is too close to the border - i.e. $centroid(x) < 0$ or $centroid(x) > imagewidth - shapewidth $ it is outside the edge, similarly you can apply for Y axis as well.
Given that you are only concerned about what falls on the edge - you should start with each edge only and start pattern matching there. Start matching the partial pattern/shape and if the partial pattern/shape does match that object IS being cut on the edge.

Here are some references that might help you formulate the problem well.

This paper is a very good to understand lot of basics about Signs/tokens you are dealing with.

Anil K. Jain and Aditya Vailaya Shape-Based Retrieval: A Case Study with Trademark Image Databases  Pattern recognition 1998, vol. 31, no9, pp. 1369-1390

There are many research elements which deals with partial or occluded shape/pattern matching.

Eli Saber,Yaowu Xu, A. Murat Tekalp Partial shape recognition by sub-matrix matching for partial matching guided image labeling Pattern Recognition 38 (2005) 1560 – 1573

Will expand this answer for more specific queries if you take this approach.


Most likely yes. FSK is already generated using many samples per symbol. Look at Matlab fskmod function.


In DSP - is it true that the only type of interfering wave occurs when frequencies match?

This is a question where theory and practice will diverge somewhat.

Theory

I assume that by "frequencies" you mean sinusoids of varying frequency. In that case, yes, you're correct, for any amplitudes $a_1$ and $a_2$ and frequencies $f_1$ and $f_2$ (where $f_1 \ne f_2$), the signals

$$
x_1(t) = a_1e^{j2\pi f_1 t} \\
x_2(t) = a_2e^{j2\pi f_2 t}
$$

are orthogonal to one another. The exact meaning of orthogonal varies from context to context, but in this sense, you can take it to mean that you could isolate the two signals from one another. If you had infinite time and infinite resources (and perhaps the ability to bend physics to your will if you want to realize it in hardware), you could build a filter that would pass $x_1(t)$ and eliminate $x_2(t)$, or vice versa. In theory.

In our theoretical utopia, the relative amplitudes of $x_1(t)$ and $x_2(t)$ don't matter either. Your theoretical processor can have unlimited dynamic range and/or numeric precision, so having one really large and one really small sinusoid isn't a problem.

Practice

Practical systems are hamstrung by several limitations that can invalidate the above claims. Among them:

If $f_1$ is really close to $f_2$, it can be difficult to resolve them as two separate sinusoidal components. If you were trying to do so digitally using short-time Fourier analysis, for example, you need a longer observation time in order to properly separate the two frequencies in your spectral estimate.

Another similar way of looking at it: if $f_1$ and $f_2$ are close to one another, and you wanted to build a filter that would pass one and eliminate the other, you would need a filter with a very sharp transition band. It can be shown that, qualitatively, the width of a filter's transition band is inversely proportional to its overall delay. So, as your transition band gets smaller and smaller, the delay required to implement the filter gets larger and larger. Furthermore, such a sharp transition band would require higher- and higher-order filters, possibly beyond the point of being practically realizable.
Practical systems do not have unlimited dynamic range. There will be some finite ratio between the maximum and minimum amplitudes that a system can reliably measure. This can limit your ability to resolve sinusoids that have widely-varying amplitudes, even if their frequencies are far apart. A very strong signal can desensitize your observer to the point where you can't see smaller signals.

Try shining a bright flashlight and a small LED at your face at the same time. Even if they are comfortably spread apart in space, you'll have a hard time seeing the LED due to the bright flashlight; your visual system has limited dynamic range.


Extracting positive frequencies of discrete-time signal

The mistake lies in the modulation of the low pass impulse response. The low pass response equals $1$ for $\omega\in(-\pi/2,\pi/2)$ and is zero otherwise, and you want a frequency response that is zero for $\omega\in(-\pi,0)$ and $1$ for $\omega\in(0,\pi)$. Consequently, you need to shift the spectrum by $\pi/2$, and not by $\pi$ as you did. The corresponding impulse response of the system becomes

$$h[n]=\begin{cases}\displaystyle\frac12,&n=0\\\displaystyle\frac12 e^{jn\pi/2}\frac{\sin(n\pi/2)}{n\pi/2},&n\neq 0\end{cases}\tag{1}$$

which can be rewritten as

$$h[n]=\begin{cases}\displaystyle\frac12,&n=0\\\displaystyle0,&n\textrm{ even},n\neq 0\\\displaystyle\frac{j}{n\pi},&n\textrm{ odd}\end{cases}\tag{2}$$

So only the sample at $n=0$ is real-valued, all samples at even indices (other than zero) vanish, and all samples at odd indices are purely imaginary.

Consequently, the real part of the impulse response is just a scaled unit impulse, and its imaginary part is a scaled discrete-time Hilbert transformer:

$$h[n]=\frac12\big(\delta[n]+jg[n]\big)\tag{3}$$

with $g[n]$ the impulse response of an ideal discrete-time Hilbert transformer.


Impulse response if initial conditions are given

The reason you're having trouble with this is because the impulse response of a system tells you something considerably different from a system's behavior when it is starting with non-zero initial values -- so no one bothers to treat this in the literature.

Problem But I think this method will give wrong answer because due to impulse input there must be some changes occur in initial conditions, so I wanted to know

A way to reconcile the ambiguity you see (and, typically, something you need to do if you're getting down and dirty with the details around the Dirac delta functional) is to make a distinction between $t = 0$, $t = 0^+$, and $t = 0^-$, where $0^+ = 0 + \epsilon$ and $0^- = 0 - \epsilon$ where $\epsilon$ is that famous infinitesimal from calculus that's greater than zero but arbitrarily small.
Then recall that the Dirac delta functional is zero for $t < 0^-$, and zero for $t > 0^+$, and -- strictly speaking -- undefined for $t = 0$

Due to change in initial conditions, which equation will change, zero input or zero state or both?
What is best approach to get solution of problem like this?

With a bit of hedging about with what you mean by initial conditions, your original idea -- that you can just add the impulse response and the system behavior with initial conditions -- is correct.
You lead yourself to this by recalling that if $\delta(t)$ is undefined at $t = 0$, then so are its effects.  So to be entirely picky and technical, the impulse response is only defined for $t > 0^+$, not before.
Then you hedge your initial conditions statement by specifying the initial conditions at $t = 0^-$.  Without that $\delta(t)$ in there, the state of the system at $t = 0^+$ would equal its state at $t = 0^-$, so (to an engineering level of rigor) this works.
With these details specified, you can just add the two responses.

Is by using Laplace transform this problem can be solved without any change in initial condition?

You can use the Laplace transform as an aid to solving this, but at least the way that I was taught to insert initial conditions into a Laplace-domain solution was to use Delta functionals and their derivatives.  You are, essentially, doing the same thing I did before, it's just that instead of taking the initial conditions at $t = 0^-$, you're finding the correct weights for $\delta(t)$, $\delta^2(t)$, etc., so that the states of your system match your initial conditions at $t = 0^+$.
Then to add in the impulse response, you just add in an extra one times $\delta(t)$.  The effect will be the same as your proposed solution.


Why does my plot of the Blackman window not match everyone else's in literature?

There is a missing $2\times$ factor inside the second cosine, and you get a centered $201$ point windows with the change of variable $x-100$. See at Wolfram:

$$ 0.42 - 0.5*\cos\left(2\pi\frac{x-100}{200}\right) + 0.08 \cos\left(4\pi\frac{x-100}{200}\right)\,.$$


I would expect that you simply cut the outer parts of the image away. Certainly, there is always a hassle to figure out the proper pixel in the center, but that is a detail of the FT-algorithm employed.

Remeber: Sampling denser than you need increases the FOV. Hence, the image you simulate should simply be larger with the original image (that you get by simulating 128 sampling points in the setting you described above) sitting at its center. Cutting off the outer void parts is therefore the technique you need, and not an interpolation technique.


There are many different ways to extract phase shift, the simplest one I think is using normalized cross-correlation using instruction 'xcorr' and then finding the index where the maximum correlation is placed. Other technique as you mentioned is using FFT and spectrum multiplication and ifft(equals to correlation in time) and/or calculation of  arctan(complex/real) of bins. Other methods like Hilbert transform is comprehensively detailed in below links:

https://stackoverflow.com/questions/27545171/identifying-phase-shift-between-signals

Best method to extract phase shift between 2 sinosoids, from data provided


This can be explained with the help of some knowledge of probability theory.
Since the output signal is a hybrid signal consists of many(say,Nfft=2048) single carrier signals, which are QAM-modulated by the complex symbols derived from the information bits. Thus leads to(according to the Central Limit Theorem in probability theory) a near Gaussian distributed signal.

As is known to us, Gaussian distribution is not bounded in its range, hence resulting in a high PAPR in OFDM signal. While, in SC_FDMA, things are different. The output signal is essentially a single carrier modulated by the complex symbol, no summation of Nfft signals here. Thus a low PAPR is obtained.

This can also explained intuitively in the time domain. Say, you have a single carrier(in SC-FDMA) modulated by A QAM modulated scheme. The peak of the signal waveform is almost bounded by the carrier amplitude it self. Whilst, in OFDM, where you have a large number of carriers, each with a different frequency. The peak of the output signal is not necessarily bounded by the highest amplitude, since at some time point, all the peaks of the carriers maybe coincide together and thus sum up to form a ultra high peak. While, at some time point, its negative counterpart can be appeared. Although both of the positive ultra high peaks and negative of it are rare in time, it indeed leads to a high PAPR.


Scipy fourier transform zero frequency spike (from DC offset) - de-meaning and hanning window have no effect

It doesn't happen with a random signal.  Your signal must have low frequency content around 0 Hz that shows up even after you've nulled out 0 Hz itself?
import numpy as np
from scipy import signal
from scipy import fft
import matplotlib.pyplot as plt

a = np.random.randn(174001)
a -= np.mean(a)

win = signal.hann(174001)
dt = 0.01
n = 174001

X = fft.fft(a*win, n=n)
freqs = fft.fftfreq(len(a), dt)
plt.plot(freqs, np.abs(X))
# plt.yscale("log")
plt.show()

If you zoom in on 0 Hz, I bet you can see that the 0 Hz spike itself goes away when you null out the mean, but the content around 0 Hz stays the same.


Java - Convert Image to Binary (Two Colours/Black & White)

Converting it to "negative" will give you single channel result (0-255). To convert this image to binary, you need thresholding. This threshold value will have a HIGH (0-255) and LOW (0-255). These values are dependent on your application (what you wish to segment). For instance, you could setup an empty image structure and set all values to '0'. Set '1' to channel values that fall between HIGH and LOW. However, remember that performing naive segmentation like this may not give you good results.


What is 'Normalized frequency in the range [0,1)', à la DTMF & Goertzel algorithm

The normalized frequency '1' corresponds to the Nyquist frequency, i.e. half the sampling frequency. So the normalized frequencies $w$ are obtained from the frequencies $f$ in Hz by

$$w=\frac{2f}{f_s}$$

where $f_s$ is the sampling frequency.

The frequency (in Hz) corresponding to a certain bin $i$ ($0\le i\le N-1$) of an FFT with length $N$ is

$$f_i=\frac{if_s}{N}\quad i=0,1,\ldots,N-1$$

For real-valued signals, the bins $i=N/2+1$ (assuming $N$ is even) to $i=N-1$ are redundant.


ESPRIT, like the French word.

For two reasons:

Out of vanity of myself, I learned it like that, I used it in presentations like that, and all I need for that to be correct is for a majority to recognize it to be correct. Hence, by stating this publicly, I make it correct! (just kidding, but when you look at language, there's no right or wrong – there's only conventions)
A. Paulraj himself pronounces it French in this presentation, and he should know :)

So, I disagree with @Envidia on how it's pronounced, but I agree with them on: It doesn't really matter.


Nothing is missing in your understanding. When we calculate power that way, it is assuming a $1 \Omega$ load and that amplitudes are given in volt, but almost nobody ever says that even in textbooks. Unless it is explicitly stated otherwise, I think it is safe to assume this.


How to efficiently calculate only the low coefficients of a zero-padded FFT

If you only a few bins then the following may be very efficient for you:
1. Simply do the DFT of at each frequency you need.
2. Use the Goertzel algorithm for each frequency in question.


40kHz would be the two sided bandwidth, in which case sampling requirement is 40khz minimum


Say you quantize a signal with $L=2^n$ levels. Then, every quantized sample is represented by $n$ bits. So, as you increase $L$, you need more and more bits to represent the samples. This means that transmitting the signal will require transmitting more bits.

Also, keep in mind that, in order to convey bits at a rate of $R_b$ bits per second, you need a signal with bandwidth at least $B=R_b/2$ hertz (assuming binary communications).

Now (this is important) assume that you want to transmit the signal in real time, like in a voice application, or in streaming. If you increase $L$, you will have to transmit more bits per second, which means a larger $R_b$, which in turn implies a larger $B$.

As an example, say you're transmitting telephone-quality voice with $L=256$, $n=8$, and with 8000 samples per second. Then, $R_b = 8 \times 8000 = 64000$ bits per second, which require a signal of bandwidth $B > 32000\text{ Hz}$.

Say you want to improve the signal quality and go with $L=1024$ ($n=10$). Now, you need $R_b = 80000$ and $B>40000\text{ Hz}$. An increase in $L$ produced an increase in $B$.

Note that, if there is no real-time requirement, you can instead keep a signal with the same bandwidth (so, same bit rate) and instead take longer to transmit. So, in general, increasing $L$ will cause an increase in bandwidth and/or an increase in transmission time.


Finding crosses and their orientation on a given image

Try to use cv2.HoughLinesP (see an example) to detect the straight lines (segments). As a result you will get the segments' endpoints, so finding all the intersections and its type will be quite simple maths.
The problem will be rather with finding the proper and universal parameters.

Theoretically you could use a morphological operation called hit-and-miss transform, but this will be difficult in your case, as the lines are hand-drawn and even after thinning the crosspoints may be difficult to be properly classified (as 'T', '+' or 'corner').

I would suggest you to make some morphological preprocessing to improve the binary image you have got.