Why is one of the modulated frequency spikes smaller? This is a simple amplitude modulation experiment. Both waves are sinusoids of the type

cos(2 x pi x f x t).

The base signal has a frequency 100Hz, and the carrier wave frequency is 1000Hz. The three figures below are the modulated signal in time and frequency domain for three different sampling periods.

For sampling frequency, Ts = 3000Hz

For sampling frequency, Ts = 4000Hz

For sampling frequency, Ts = 10000Hz

I expected both frequency components at 900Hz and 1100Hz to be of the same height. Except for the case when sampling frequency is 4000Hz, in the other cases this is not so. Why?

The reason is because of spectral leakage. Your frequencies do not fall exactly on the FFT bins- this causes the energy to appear in nearby adjacent FFT bins as sidelobes. The most spectral leakage occurs when the frequency of a signal falls exactly between two FFT bins. Another way to think of it, is that the signal is not periodic w.r.t. the length of the FFT you are using.

Your two frequencies are appearing at two different positions with respect to the adjacent FFT bins - so they are suffering different amount of leakage. Thus, the peaks have different magnitudes. In the case where the frequencies have the same magnitude, they are located the same distance away from the closest FFT bin.

You can try adjusting your sampling rate and FFT size so that the resulting frequencies will fall exactly on an FFT bin. With the 3000 Hz sampling rate use a 3000 point FFT - corresponding to FFT bins every 1 Hz. With the 4000 Hz sampling rate use a 4000 point FFT - again, corresponding to FFT bins every 1 Hz.

In your 4000 Hz case, you are probably using a 4196 point FFT. The FFT bins are close to being every 1 Hz, so the frequencies at 900 Hz and 1100 Hz are almost exactly on a FFT bin and thus have almost equal amplitude.


How is the energy conservation principle preserved in the case mentioned below? Suppose there are 2 sinusoidal signals $\cos(\omega t)$ and $\cos(\omega t)$, which implies that their individual power is $1/2$ unit and total power is $1$ unit. When added, the signal $2\cos(\omega t)$, however, has power $2$ unit.

Obviously,

$$\rm Energy_{signal} = Power_{signal}\times Time_{period}$$

So, how does the principle of energy conservation hold in this case?

Power is homogeneous to a second degree quantity. When the data is doubled, the power is quadrupled.

The energy here (classicaly) is defined in an Hilbert space as a squared $L_2$ norm, deriving from some inner product $\langle\cdot,\cdot\rangle$, assuming a real Hilbert space,  $\|x\|^2 = \langle x,x\rangle$. Hence:  $$\|x+y\|^2 = \langle x+y,x+y\rangle = \|x\|^2 + 2\langle x,y\rangle+  \|y\|^2 \,.$$

Now set $y=x$, you gracefully get $\|x+x\|^2 = 4\|x\|^2$.


Calculating the inverse DTFT of a signal There is a signal $y[n]$ with a differentiable DTFT $Y(e^{i\omega})$. How do I find the inverse DTFT of $i\frac{dY(e^{i\omega})}{dw}$ in terms of $y[n]$ (where of course $i =  \sqrt{-1}$)?

You can develop it relatively easily, anyways, here's a screenshot from wikipedia


Use samples of Fourier transform as DFT? Consider the LTI system given by: $H(z) = 1 - \frac{1}{2}z^{-1}+\frac{3}{4}z^{-2}$

Let $x[n] = (\frac{1}{2})^nu[n]$ be the input to the system. We want to find the output for $n = 0,1,...,N_a$, using the DFT of a relevant portion of the input and samples of $H(z)$ on the unit circle, at $z = e^{j\frac{2\pi}{N_b}k}$, for $k = 0,...,Nb-1$.

Now, because $H(z)$ is FIR, the samples of $H(z)$ on the unit circle are actually the values of the DFT of length $N_b$ of $h[n]$. I can work out the appropriate length of the DFT of $x[n]$ to use and it would work and
using the circular convolution property, find the actual values of the output.

What if instead of a FIR, I'm given a causal IIR such as $H(z) = \frac{1}{1-\frac{1}{2}z^{-1}}$?

Can I use the samples of this new $H(z)$ on the unit circle just as in the first case? My guess is yes, because the system is stable (because it has a pole at $\frac{1}{2}$ and it is causal) and then the Fourier transform exists, so I can obtain the samples of the Fourier transform and use them as DFT of $H(z)$, (I believe) it would be like truncating $h[n]$ and then using its DFT.
Now, what if the system was unstable and causal, for example:

$H(z) = \frac{1}{1-2z^{-1}}$?

The Fourier transform clearly doesn't exist, and I can't truncate $h[n]$

Check my answer on a similar question : H(z) as N-DFT

Only thing to add it that for both the systems, it is possible to compute $H[k]$ by figuring that the numerator and denomenator are as follows :

Numerator = $[1]$, Denominator = $[1, -\frac{1}{2}]$.
Numerator = $[1]$, Denominator = $[1, -2]$.

Take $N_b$-DFT of numerator and $N_b$-DFT of Denominator and point-wise division of these will give $H[k]$. Take IDFT to get $N_b$ samples of $h[n]$. Convolve with $x[n]$ to get output.

This is possible for stable and unstable both systems, because we are truncating $h[n]$, which will inherently make the system stable. The main requirement is that the system is causal.


Obtaining the magnitude of the frequency response by plugging $e^{jω}$ into the z-domain transform function? I am reading a text on discrete signal processing, which states that the frequency response of a signal can be obtained by plugging the value $e^{jω}$ into the z-domain transfer function $H(z)$. In other words:

$H_{IIR}(ω) = H(z)|_{z=e^{jω}}$

In the example given in the text, it is given that:

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

Next, the text presents the following equation for the magnitude of the frequency response when $ω = \frac{\pi}{2T}$:

$|H_{IIR}(\frac{\pi}{2T})| = \frac{|j|}{|j - \frac{1}{5}||j + \frac{1}{2}|} = \frac{1}{\sqrt{\frac{1}{25} + 1}\sqrt{\frac{1}{4} + 1}} = \frac{\sqrt{130}}{13}$

I don't understand how the author derived the above equation. That is, what kind of wizardry leads to the following equations?

$|H_{IIR}(\frac{\pi}{2T})| = |(\frac{1}{1 - \frac{1}{5}e^{-j\frac{\pi}{2T}}}) (\frac{1}{1 + \frac{1}{2}e^{-j\frac{\pi}{2T}}})| = \frac{|j|}{|j - \frac{1}{5}||j + \frac{1}{2}|} = \frac{1}{\sqrt{\frac{1}{25} + 1}\sqrt{\frac{1}{4} + 1}}$

Does it have something to do with magnitude and a special identity regarding the constant $e$?

It's  derived using Euler's equation http://en.wikipedia.org/wiki/Euler's_formula. You start with $$H(z) = (\frac{1}{1 - \frac{1}{5}z^{-1}}) (\frac{1}{1 + \frac{1}{2}z^{-1}})$$

Then plug in $e^{-jwT}$ for $z^{-1}$ and you get

$$H(z) = (\frac{1}{1 - \frac{1}{5}e^{-j\frac{\pi }{2}}}) (\frac{1}{1 + \frac{1}{2}e^{-j\frac{\pi }{2}}})$$

Now according to Euler we simply have $e^{-j\frac{\pi }{2}} = -j$ so the whole thing simplifies to

$$H(z) = (\frac{1}{1 - \frac{1}{5}\cdot (-j)}) (\frac{1}{1 + \frac{1}{2}\cdot (-j)})$$.

Utilizing $(-j) \cdot j = 1$ we can multiple each part of the fraction with j/j and get
$$H(z) = (\frac{j}{j - \frac{1}{5}}) (\frac{j}{j + \frac{1}{2}})$$

Then we can take the magnitude as

$$\left | H(z) \right | = \frac{\left | j^{2} \right |}{\left | j-\frac{1}{5} \right |\cdot \left | j+\frac{1}{2} \right |}$$

This isn't exactly the same formula that you have given but since magnitude j and j^2 are the same, the final result comes out to be the same.


How to Solve Image Denoising with Total Variation Prior Using ADMM? I was looking at some articles or Wikipedia on denoising images using the Total Variation norm. The setup is the Rudin Osher Fatemi (ROF) scheme, and the corresponding equation is:
$$
F(u)=\int_{\Omega}|D u|+\lambda \int_{\Omega}(K u-f)^{2} d x
$$
Some of the sources mentioned using the ADMM optimizer to solve this denoising problem. But I was hoping that someone might be able to direct me to some code to show an implementation of this approach. Code in MATLAB, Julia, or Python would be excellent, just something to get started with.
Thanks.

I have done a bit of this myself and you'd need to adapt.
There is a Douglas Rachford self implemented and a primal dual approach here implemented in Recovery of Fusion Frame Structured Signal via Compressed Sensing.
Note that Clarice Poon (Bath University) had some nice tutorials on it.
Another source is the Numerical Tours from Gabriel Peyre. See Denoising by Sobolev and Total Variation Regularization.


Complex Least Squares Approximation In the case of frequency domain FIR filter design, the error function given by :

$$E(\omega)=H(e^{j\omega})-D(e^{j\omega}) \tag{1}$$

is a linear function with respect to the unknown filter coefficients $h[n]$. Hence, the Least Squares frequency domain FIR filter design problem is a linear LS problem which can be solved by computing the solution to a system of linear equations.  (from this: Algorithms for the Constrained Design of Digital Filters with Arbitrary Magnitude and Phase Responses ).

I have a three questions :

How the error respect the unknown filter coefficient ?
Does equation (1) is the same : $E(\omega)=D(e^{j\omega})-H(e^{j\omega})~$?
Does this case is available with negative domain of frequency ?

Thank you in advance.

Let's assume we have $ m $ points in the Fourier Domain (Discrete Fourier Domain) for our reference frequency response. We'll designate their values by $ y \in \mathbb{C}^{m} $.

Also we need a filter of $ n $ coefficients to approximate the reference response as good as it can. We'll use $ h $ for that.

We need now the Fourier Matrix, $ F \in \mathbb{C}^{m \times n} $ in order to write the following linear complex least squares problem:

$$ \arg \min_{h} \frac{1}{2} {\left\| F h - y \right\|}_{2}^{2} $$

You could easily write this also in a Sum Form.
The nice thing about it is the solution is give by the classic Least Squares solution - $ \hat{h} = {\left( {F}^{H} F \right)}^{-1} {F}^{H} y $.

Of course, in practice, you need many point in frequency domain to add more and more restrictions. You could also write a continuous version of this problem using the DTFT.


How to detect start and finish of temperature control in temperature time series I have a huge dataset containing temperature data inside a building. I want to extract the time that the building starts and stops controlling the temperature (approximately around the vertical black lines in the pictures). In this case the temperature control starts and stops at the same time. However, I want to detect the the same for a building where the start and stopping time is variable, per day and per area inside the building. How would you go about finding a solutions for such a problem?
I can quite reliably detect the first line by detecting when the absolute velocity / acceleration reaches above a certain threshold, but this method seems less reliable for detecting the stopping point.

Smaller subset to better show the behavior of the data:

One option: It's pretty clear that in the uncontrolled range the temperature follows a  simple exponential heating/cooling curve. Something like
$$T_{n+1} = T_n + (T_{Ambient} - T_n) \cdot e^{-\Delta t / t_{Building}} $$
where $\Delta t$ is the time difference between samples and $ t_{Building}$ the thermal time constant of your building (which may be slightly different at different sensor locations).
You could simply test the difference of each new sample to the model prediction and set an threshold of the allowable error. The only two model parameters are ambient temperature and thermal time constant of the building.
The thermal time constant is a function of the thermal resistance of the building and it really doesn't change unless someone opens/closes windows or there is some actual structural work happening. Looking at your data that should be easy to estimate. Ambient you can try to get directly either with a sensor or looking at a reliable weather report. You can also just use the data from the first few sensor to estimate ambient and then test whether the other sensors agree or not.


Initial Rest Condition for LCCDE causal LTI systems I am self studying Alan Opennheim's course Signals and Systems. I am a math major and have no background in EE.
I understand that for a linear constant-coefficient difference equation (LCCDE) system to be linear its auxiliary conditions must be 0.
I also understand that for the system to be causal every output $y(t)$ corresponding to input $x(t)$ s.t. $x(t) = 0$ for $t<t_0$ must satisfy $y(t)=0$ for all $t<t_0$.
What i am having trouble understanding is why this imposes the initial condition $y(t_0)=0$ for the response to the unit step function. Couldn't we have set the initial condition to be $y(t_1)=0$ for $t_1>t_0$ and kept the system causal and LTI?
In the 6th lecture in the video series Alan Oppenheim finds the unit step function response to the causal LTI system described by $y'(t) +ay(t) = x(t)$. He imposes the initial condition $y(0)=0$, as i mentioned above i can't understand why this is imposed by the properties of the system. He then goes on to find the unit impulse response, and finds that it is $e^{-at}u(t)$ where $u(t)$ is the unit step. Clearly this function does not satisfy the initial condition which were imposed on the unit step response, I am having trouble understanding why this is valid.
Any help is appreciated, Thanks!

If you are given an input $x(t)$ with $x(t)=0$ for $t<t_0$, and you specify an initial condition $y(t_1)=0$ for $t_1>t_0$, then the resulting system is generally non-causal, because we already know the system's response at $t_1>t_0$, regardless of the input signal in the interval $[t_0,t_1]$.

For a system to be linear (in the sense used by Oppenheim) and causal, we need to specify the initial condition at the point $t_0^-$, i.e., we require that the left-sided limit $y(t_0^-)$ satisfies $y(t_0^-)=0$.

For the given example, if an impulse is applied at $t=0$, we require $y(0^-)=0$ for the system to be linear and causal. Note that there is no requirement on the right-sided limit $y(0^+)$.

Clearly, the response $y(t)=e^{-at}u(t)$ satisfies the initial condition $y(0^-)=0$.


"fred harris rule of thumb": I would like to cite  "fred harris rule of thumb" for estimating filter order.  Is there a book or paper where I can find this rule of thumb?  Thanks a bunch.

You can find the formula for (at least) one "Harris method" or "Harris approximation" in Multirate Signal Processing for Communication Systems, Fredric J. Harris, 2004, page 216, equation  (8.16), with details in Chapter 3. F. Harris calls it "harris approximation" (without capital for the name), namely:

$$ N = \frac{f_{SMPL}}{\Delta f}\frac{Atten(dB)}{22}$$

Also in Digital Communications with Emphasis on Data Modems: Theory, Analysis, Richard W. Middlestead, page 752, equation B.39.

Former questions might provide you with additional insights and references:

How many taps does an FIR filter need?
Filter Order Rule of Thumb


The Zero-Crossing rate threshold for a voiced/unvoiced decision I've implemented a function that calculates the zero-crossing rate for a given signal. I've used this same function to calculate the pitch. To differentiate voiced signals from unvoiced with reference to ZCR: A high ZCR means that the signal is unvoiced and a low ZCR means that it is voiced. My question is whether there is a threshold above which we can consider that a signal is unvoiced.

Plain ZCR criterion is not enough for robust and accurate voiced/unvoiced separation. That being said, your threshold should be adaptive, there is no fixed threshold that works well for all speech waveforms. It depends on the approach you follow but a statistical threshold should do the work most of the time. You can google your question to find relative papers and see if their approach suits your purpose.

Just a note, ZCR is not the best choice for pitch estimation - it just gives a rough estimation.


What is the effect of wiping the right half of Fourier Transform bins? I'm trying to change the pitch of a signal using a Fourier Transform (FFT) followed by an Inverse Fourier Transform (IFFT).
I've found many examples, some of which zero out the right half of the real and imaginary bins before changing the pitch. For example, if the signal was $8192$ bins, the real and imaginary parts from $4096$ to $8192$ are set to $0$. This seems to make the math for pitch changing easier, but reduces the volume by one half. This seems to be corrected by doubling the magnitude.
I'm wondering what effect wiping the right half of the bins has on the final signal, apart from reducing the volume. I am confused why these bins exist if they can be wiped without affecting the final signal too much.

The discrete Fourier transform (DFT) transforms a sequence of generally complex-valued numbers to another sequence of generally complex-valued numbers. If the input to the DFT is real-valued - as is the case for a sampled audio signal - then the right half of the DFT is redundant. If $X[k]$ is the length $N$ DFT of a real-valued length $N$ sequence $x[n]$, the following holds:
$$X[k]=X^*[N-k],\qquad k=0,1,\ldots, N-1\tag{1}$$
where $^*$ denotes complex conjugation.
For even $N$, the values $X[0],X[1],\ldots,X[N/2]$ completely represent the original sequence $x[n]$, so the DFT coefficients $X[N/2+1],\ldots,X[N-1]$ can be discarded. For odd $N$, the coefficients $X[0],\ldots,X[(N-1)/2]$ are sufficient.
Note that there are efficient routines for computing the DFT of real-valued input sequences, which only compute the non-redundant DFT coefficients.


I want to get impulse response of causal system I want to get impulse response of this causal system on first photo.

So I got frequency response on second photo from this system.

and now I am trying to transform this frequency response to impulse response on third photo.

but I can't invert it because there is a square of s on numerator.

How can I invert it to impulse response?

Several comments:

First, the usual symbol for the Z-transform is $z$ rather than $s$ (which is common for Laplace transform).

It is more common to express $H(z)$ in terms of $z^{-1}$. Note that the causal inverse transform for a single pole at $\alpha$: $H(z) = 1/(1-\alpha z^-1)$ is $h[n] = u[n] \alpha^n$.

So your transform $H(z) = \frac{3}{1 - 1/4 z^{-1} - 1/8 z^{-2}}$ can be expanded to the sum of two terms like the single pole above.


Cross-correlation of filtered random processes I have a wide-sense-stationary (WSS) process $\{x(t)\}$ and two linear filters with impulse functions $h_1$ and $h_2$.

Let $\delta(\omega)$ be the power spectrum of $\{x(t)\}$ and $$H_1:\omega\mapsto H_1(\omega)$$ and $$H_2:\omega\mapsto H_2(\omega)$$ the transfer functions of the filters.
The outputs of the filters are denoted $$y_1(t)=(x \star h_1)(t)$$ and $$y_2(t)=(x \star h_2)(t),$$ where $\star$ denotes the convolution.

How can we compute the correlation of $\{y_1(t)\}$ and $\{y_2(t)\}$ and when are these two random variable uncorrelated?

As an addition to Dilip's answer I'll show you how to derive that result:

$$\begin{align}R_{y_1,y_2}(\tau)&=E[y_1(t+\tau)y_2(t)]\\&=E\left[\int_{-\infty}^{\infty}x(\alpha)h_1(t+\tau-\alpha)d\alpha\int_{-\infty}^{\infty}x(\beta)h_2(t-\beta)d\beta\right]\\&=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}E[x(\alpha)x(\beta)]h_1(t+\tau-\alpha)h_2(t-\beta)d\alpha d\beta\\&=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}R_x(\alpha-\beta)h_1(t+\tau-\alpha)h_2(t-\beta)d\alpha d\beta\\&\stackrel{\gamma=\alpha-\beta}{=}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}R_x(\gamma)h_1(t+\tau-\alpha)h_2(t-\alpha+\gamma)d\alpha d\gamma\\&\stackrel{\zeta=\alpha-t}{=}\int_{-\infty}^{\infty}\underbrace{\int_{-\infty}^{\infty}R_x(\gamma)h_2(\gamma-\zeta)d\gamma}_{(R_x\star h_2^-)(\zeta)}\; h_1(\tau-\zeta) d\zeta\\&=(R_x\star h_1\star h_2^-)(\tau)\qquad\qquad\qquad (1)\end{align}$$

with $h_2^-(t)=h_2(-t)$.

The cross-spectral density $S_{y_1,y_2}(\omega)$ is the Fourier transform of the cross-correlation function:

$$S_{y_1,y_2}(j\omega)=S_x(j\omega)H_1(j\omega)H_2^*(j\omega)\tag{2}$$

where $S_x(j\omega)$ is the power spectral density of $x(t)$, and $H_1(j\omega)$ and $H_2(j\omega)$ are the frequency responses of the two filters.

Using $(2)$ it is straightforward to define a condition on $H_1(j\omega)$ and $H_2(j\omega)$ such that the cross-spectral density, and, consequently, the cross-correlation function become zero.


Why phase information is usually ignored after Fourier transform? Why phase information is usually ignored after Fourier transform?

I have read in a tutorial [ref: FFT Tutorial] that interpretation of phase is usually challenging, so scientists do not show it or use it for their analysis. I would like to know when one can confidently ignore phase information without losing key information in his/her judgements.

Simply put:

Because in interesting cases, it remains difficult to interpret, and people tends to use other tools
Yes, potentially a lot of non-stationary features, timings, onsets, discontinuities

For a basic use, like measuring the relative energy of different frequency component, the phase is of little use. As waves tend to last over several samples, sliding/overlapping versions can help retrieving timing/localization of features.

However, phase contains a lot on information, especially on images, and perhaps plays an important  role on perception (e.g. The importance of phase information for recognizing natural images). The importance of phase is recognized in many works, see for instance The Importance of Phase in Image Processing, and it is still an active topic. Multiscale phase might detect onsets, contours, while the distribution of phase coefficients could be used for texture retrieval. The phase computed is wrapped  between $-\pi $ and $\pi $, and unwrapping it in 2D or 3D seems to remain challenging. Hence, some resort to complex transformations that are more localized than Fourier bases, like filter banks or "wavelets", to capture phase information at the correct scales.


Given a noise-free channel with limited bandwidth, why is the channel capacity infinite? In the Shannon Hartley theorem, if N = 0 but B is limited, the theorem tells us that C is infinite.  Why is this intuitively true?

If bandwidth is limited, I can only post so many pulses per second over the channel.

Why is this intuitively true?

It's a purely hypothetical mathematical formula result. So it doesn't even need to be intuitive, at all.

It's very normal for humans to look for intuition in abstract things, but an observable noise-free channel doesn't exist, so meh.

But if you insist: in a noise free channel, a single symbol could take one of $N$ values, with $N$ becoming arbitrarily large, because, due to the absence of noise, you could tell arbitrarily fine voltage differences apart. With that, a single (equally probable) symbol $x$ would have Shannon info $I = -\log_2 p_x = -\log_2 \frac1N = -(\log_2 1 - \log_2 N) = \log N\rightarrow \infty$ information.

EDIT: I saw you just got the same answer not an hour ago:

  But in terms of pure math, there is no ($\aleph_0$) finite limit on the amount of information that a single real number (one sample) can contain.


What reference power should I use? My professor asked me to convert following powers into decibels:
15W, 25W, 33W, 450mW, 24MW, 44 GW,120pW.
I know that decibels are about ratio power/reference_power.
I don't have reference power specified here. Is there something like default value? What should I do?

According to https://en.m.wikipedia.org/wiki/Decibel_watt

dBW is relative to 1 W and dBm is relative to 1 mW.  I’d just make sure to provide the correct units, but maybe you could ask your teacher if it matters.


Vocoder carrier wave bleed I made a terrible vocoder earlier and got some help on it here:

My vocoder sounds bad

My new implementation is much better. I'm using the band pass filter described in the document suggested on the old question. It seems to work much better, but the carrier wave is bleeding heavily into the output signal to the point where I can faintly hear some good vocoder sound effects behind the sound of the carrier wave. This bleed gets louder when the input volume is higher. Is this a common and easily fixed problem or should I give more info for you guys to help me debug this?

I have, since my troubles with vocoders, learned a lot about signal processing in school and gone back to revisit some of my older projects and fix them. The problem with my vocoder implementation was that I was defining bandwidths in octaves, but linearly spacing my filters. There would be more overlap in higher frequencies and less overlap in lower frequencies. I'm not sure how I got carrier bleed, perhaps a bad envelope follower, but re-implementing with better filters was the real solution.


How is channel frequency response/channel state information related to free space path loss? As we know, the path loss can be described by the equation as
$$
\frac{P_{r}}{P_{t}}=D_{t} D_{r}\left(\frac{\lambda}{4 \pi d}\right)^{2}
$$
In the OFDM system, such as Wi-Fi, there are infinite paths exist, assume it is in free space, and only one path exist, for the specific path (LoS path), the channel frequency response(CFR)/channel state information(CSI) describes how is the channel affect the transmitted signal, it calculated by:
$$
H = \frac{Y}{X}
$$
where Y is the received signal at receiver and X is the transmitted signal at transmitter. the signal power is square of the amplitude of the signal, that is
$$
P_{t}=|Y|^2
$$
$$
P_{r}=|X|^2
$$
I have two questions:

Is the following equation holds?
$$
\left|H_{}\right|^{2} = \frac{P_{r}}{P_{t}}=D_{t} D_{r}\left(\frac{\lambda}{4 \pi d}\right)^{2}
$$
According to the first equation, is the following equation holds?
$$
\frac{X}{Y}=\sqrt{D_{t} D_{r}\left(\frac{\lambda}{4 \pi d}\right)^{2}}
$$

Thanks in advance.

In the OFDM system, such as Wi-Fi, there are infinite paths exist, assume it is in free space, and only one path exist,

Exactly. So if you use OFDM in a scenario where FSPL describes your transmission sufficiently, then you're doing something wrong, because you don't need OFDM, as you have a flat channel.
We can still describe large scale fading in multipath channels using models similar to your FSPL formula, but neither of your equations hold! Then, we're describing distance-dependent stochastic properties of the random variables that make up your CSI.
You need to cross over from your deterministic "this channel has a defined attenuation" model to a stochastic "and here's how we describe the properties of random variable, whose realizations are channels".


Temporal Orthogonality I'm studying Digital Communication and I couldn't understand the relationship between temporal orthogonality and sampling theorem.

Definition of the Temporal Orthogonality is as follows:

$$\mathrm{ACF}= \int s_i(t)\cdot s_l(t-kT) \ \mathrm{d}t= E_{il}\cdot \delta [k]$$

In the text, it is stated that:

  [...] This means that with respect to the represented symbol sequence
  lossless discrete-time signal processing is possible using only one
  sample per symbol interval, even if the sampling theorem is not
  fulfilled.

We are at the output of the transmitter.

Say you have two pulses, $p(t)$ and $q(t)$, such that the following three statments are true:

\begin{align}
\int_{-\infty}^\infty p(t) p(t) dt &= 1,\\
\int_{-\infty}^\infty q(t) q(t) dt &= 1,\\
\int_{-\infty}^\infty p(t) q(t) dt &= 0.
\end{align}

Then these pulses are said to be orthonormal. Now, let us say you want to transmit two information-bearing symbols $a_1$ and $a_2$. (The symbols carry information when the receiver does not know their values in advance). You can acomplish this by transmitting the signal $s(t) = a_1 p(t) + a_2 q(t)$.

How does the receiver recover the symbols? By performing these operations:

\begin{align}
a_1 &= \int_{-\infty}^\infty s(t) p(t) dt, \\
a_2 &= \int_{-\infty}^\infty s(t) q(t) dt. \\
\end{align}

I think it's clear that this scheme can be extended to any number of symbols; all you need to do is find a set of orthonormal pulses of the required size. The most common approach is to find a "prototype pulse" $p(t)$ with the property that, when time-shifted $lT_p$ seconds with $l$ integer, it is orthonormal to all other pulses time-shifted $mT_p$ seconds, for any integer $l \neq m$.  So, you can transmit any number of symbols with the signal

\begin{align}
s(t) &= a_1 p(t - T_p) + a_2 p(t - 2T_p) + a_3 p(t - 3T_p) + \cdots \\
&= \sum_k a_k p(t-kT_p).
\end{align}

The symbols $a_k$ are recovered in the same way as before, but this operation can be implemented with a filter "matched" to the prototype pulse $p(t)$ (see numerous questions related to matched filters on this site). When the matched filter's output is sampled at time $t=kT_p$, you obtain back the symbol $a_k$.

Note that the filter's output is not necessarily sampled at a rate that meets Nyquist's sampling theorem (I have said nothing about the bandwidth of $p(t)$ or $s(t)$), but that is OK because you're not interested in reconstructing $s(t)$ from its samples, all you want is the symbols $a_k$.


Sum of cosines as squashed as possible I have 10 cosine waves of the form $\cos(2\pi(x f_n+\varphi_n))$, where $f_n$, the frequency, is an integer from 1 to 10 and $\varphi_i$, the phase, is a number from 0 to 1:

$$s(x) = \sum\limits_{n=1}^{10} \cos\left(2\pi (x n+ \varphi_n)\right),
\quad \varphi_i \in [0,1]$$

If I set all the phases to 0, $m = \max\limits_x \lvert s(x)\rvert$ peaks as high as 10. All phases at 0.25 peaks at 7.596. Setting the phase to any linear function of freq just moves the waveform, keeping the peak the same. But making the phases random can sometimes get the peaks below 5!

Is there a more systematic way (other than picking random phases and hoping for the best) to minimize the peak of these cosines?

Edit:
I think the peak has to be bigger than $\sqrt{5}$ (amplitude of square wave with same RMS as signal above). Also, after 3 million attempts of brute force guessing, the lowest peak was 3.576. The phases still don't have any sort of pattern to them, though.

You can get pretty good results with an approximate linear chirp spread out over the period:

$$s(x) = \sum\limits_{n=1}^{N} \cos\left(2\pi \left(x n+ \frac{n^2}{2\times N}\right)\right), \quad N=10$$

$$m \approx 4.28$$

Group delay is the derivative of the phase with respect to frequency. For a linear (frequency as function of time) chirp, the phase should be a quadratic function of frequency. With discrete frequencies we don't get exactly that, but all frequencies seem to roughly find their slot in time if we sample the phase function. I wonder if a ping pong chirp would be possible here and if it would give a flatter envelope.


Inverting Sensor Transfer Functions? When you see software packages that read in data from sensing equipment, does the software do something to inverse the transfer function?

Correct me if I'm wrong in my understanding, but this is how I see it:

The sensor takes in input signal $x(t)$ and the sensor has some frequency response governed by the Transfer function of the sensor which alters the signal $x(t)$ to some deformed signal $y(t)$. So if anybody wanted to know about the actual signal would they inverse the Transfer function? And is this on top of other calibration values? Gain etc..?

I think your question has nothing to do with "software", so I'll ignore this part, and only discuss the calibration part of it. Then again, your question is very broad and cannot be answered with yes or no: It really depends.

In some applications you might want to try to invert as much of the transfer function of your sensor as possible. In other applications you might not care about it that much and leave it in or only partially compensate it.

What's though true is that in general you cannot completely eliminate the effect of your sensing equipment. By the way it filters the "true" signal $x(t)$, some of its components might be filtered out and be lost completely. As a simple example: practically all real systems have some sort of low-pass characteristic so that very high frequency components will not do anything to it and thus never show up in $y(t)$. If you excite a 1 kHz low-pass system with a 1 GHz input signal, you cannot hope to get it back by inverting the low-pass characteristic, you'd have to divide by zero essentially.

There may be other reasons you would not want to compensate. Say your system does have a very small gain at some particular frequency. Then if you try to invert this behaviour, you'll need to amplify this particular frequency a lot. This would mean that any unwanted component, such as noise, would be strongly amplified at this particular spot. In applications like this you might want to do a partial calibration, where you do correct for the transfer function in a certain frequency region of interest only, and ignore the rest of it.


Verifying Linear Time Invariance I have a system of the form:

$$T(x(n))=x(n)+3x(n-2)-5x(n-3)x(2n)$$

I claim that

$$T(x(n-k))=x(n-k)+3x(n-k-2)-5x(n-k-3)x(2n-2k),$$
$$y(x(n-k)= x(n-k)+3x(n-k-2)-5x(n-k-3)x(2n-2k)$$

and

$$T(ax_1(n)+bx_2(n))=ax_1(n)+bx_2(n)+3\left( ax_1(n-2)+bx_2(n-2)\right) - 5 \left(ax_1(n-3)+bx_2(n-3) \right)\left( ax_1(2n)+bx_2(2n)\right),$$
$$aT(x_1(n))+bT(x_2(n))) = ax_1(n)+3x_1(n-2)-5x_1(n-3)x_1(2n) + ax_2(n)+3x_2(n-2)-5x_2(n-3)x_2(2n)$$

Hence, the system is not linear, but is time-invariant.  I'm I applying the formulas correctly?

The term $x(2n)$ is not time invariant: it gives the even samples of $x$. Shift the input by 1 sample and you'll get the odd terms, which are different from the even terms shifted 1 place.

Your reasoning above is not correct. You should compare $T(x(n-k))$ to $T(x(n))|_{n-k}$ (that is, the original output, shifted by $k$ samples).

Edit:
These type of exercises are a common source of confusion because of the following abuse of notation: $x(n)$ is not a signal, but the nth sample of the signal, and having $n$ inside of the expressions is formally wrong. Even Oppenheim's book makes this mistake when working with LTI systems and impulse response.

The correct way for expressing transforms between signals is as follows:

Let $x$ and $y$ be signals; $T\{\cdot\}$ be the transformation (system) under scrutiny; $R_k\{\cdot\}$ be a delay of $k$ samples.

So, in general you would write the relation between input and output signals as $y = T\{ x \}$, where the whole signals are involved. Using $n$ inside these expressions collapses the signal to a single sample and this is wrong. And having the same variable $n$ on both sides of the equation only leads to confusion.

Then, a system is time-invariant if $R_k\{T\{ \cdot \}\} = T\{ R_k\{ \cdot \}\}$ for all integer $k$.

This is why I gave the example above with even-odd samples, in order to avoid the usual confusion with $n$ inside the expressions.

Now, going back to the system $T$ that evaluates to $y[n] = x[2n]$ for each sample $n$:

Let $x_k = R_k\{x\}$ and $y_k = R_k\{y\} = R_k\{T\{x\}\}$.

Evaluating all these signals at sample $n$: $x_k[n] = x[n-k]$; $y_k[n] = x[2n - 2k]$.

Now, $T\{x_k\}$ at sample $n$ is $x_k[2n] = x[2n - k]$ which is different from $y_k$.


Speech recognition: words vs. phonemes One could train a speech recognizer with words or phonemes. I guess that they both have their pros and cons:

with words there would be a lower error rate as phonemes can be put together building a non-existing word (this case would be impossible with word training)
using words would require a larger database and more space and processing time
phonemes are way more general: with all phonemes you can build all existing words, but with all the words in the database you can't recognize all the words in the language

I'm wondering if there is a reason to prefer one over another, or if in some cases one is better than the other. Do both methods have practical applications? Which ones? What are their pros and cons (apart from the ones I named, which may be wrong so feel free to correct me)?

Most modern speech recognition systems are based on Hidden Markov Models. These are probabilistic models that require a previous training.
So here is simple comparison:

By using phonemes you have a reduced number of possibilities (around 50 depending on the language), you get a lot of data to train each phoneme model, by combining several models you can have word models.
By using word models since the beginning, you would need to have thousands of models to train and probably not enough data.

If you are working on a very simple recognition system, with a very limited set of words, for example, recognition of numbers, then maybe you can use words models that don't require HMMs.


Complex normal Gaussian noise I would like to create complex normal Gaussian noise with dimensions $(M,N)$. The noise should be zero mean and the variance should be equal to one.

How can I do so?

You should use the the standard formula:

s = randn(m, n) + 1i*randn(m, n);

And as pointed out by MBaz, the output should be scaled accordingly by $\frac{1}{\sqrt{2}}$

s = s/sqrt(2);

More on that topic you can find for example in a book Digital Media Processing by Hazarathaiah Malepati:

It should generate exactly what you want.


DFT shift theorem proof Learning DSP on my own time. Can't figure out the proof for DFT shift theorem which states the following:

Given, $x[n]$ to be a periodic with period $N$, $\text{DFT}\{x[n]\} = X[k]$, then
$$
DFT\{x[n-a]\} = e^{-j\frac{2\pi}{N}a}X[k]
$$

I found a proof here, but I can't figure out how did they leap from $$\sum_{m = -\Delta}^{N-1-\Delta}{e^{-j\frac{2\pi}{N}k(m+\Delta)}x[m]}$$
to
$$
\sum_{m = 0}^{N-1}{e^{-j\frac{2\pi}{N}km} e^{-j\frac{2\pi}{N}k\Delta} x[m]}
$$

Specifically how did they change the limits of the summation without changing the variable?

For an $N$-periodic sequence $y[n]$, the sum over one period is always the same, no matter how you choose the first index:

$$\sum_{n=0}^{N-1}y[n]=\sum_{n=\Delta}^{N-1+\Delta}y[n]$$

for any integer $\Delta$. Try to verify that the sequence you sum over is indeed $N$-periodic.


Spectral Leakage, Phase and other questions about FFT I'm currently implementing a voice recognition algorithm using FFT and now some questions came up to me:

Spectral Leakage: I know what it is, why it appears and how to solve it using a Windowing Function, but is there another work around for this issue?
Phase: I know what it is, but where can I use it and how?
Getting the real amplitude (I know it has something to do with spectral leakage, but I will ask this separately): Let's say I have a frequency spectrum with a frequency resolution of 250 (sampling rate = 32000 Hz, frame size = 128) and the input I sent to FFT is a sinusoidal function whose frequency is 567.5 Hz (thus, the number of cycles are non-integer, causing the spectral leakage). My spectrum would be something like this:

So my question is: How can I retrieve the original amplitude of the frequency 567.5 Hz?

Windowing artifacts can be considered as due to a loss of information about a larger data set (e.g.what the window cut out or zeroed).  The "spread" of the Sync is related to the shortness of the data window (what's left after the information loss).  Thus other information about the original data from outside the FFT results can be used to help "solve" this problem (more data, perhaps used with longer FFTs, results in a narrower Sync main "hump").

Phase information within each FFT aperture plus the distance between adjacent or overlapped frames of FFT data is used in to refine phase vocoder frequency estimates or to aid in signal reconstruction.  Phase differences between frames taken from differently placed microphones can be used to help triangulate sound source location.

Peak amplitude can be estimated by interpolation equations, such as parabolic interpolation from 3 points, or better, from (windowed) Sinc kernel interpolation with successive approximation to find the zero derivative point.


Bandwidth of PSK modulated signal What does it mean that the bandwidth of PSK modulated signal without applying any filtering or shaping of the pulse, is equal to the sampling rate with "1st lobe" equal to two times of symbol rate? What is "1st lobe" here? And how it equals to two times of symbol rate?
I've also read, that with applying pulse shaping techniques, "the null-to-null bandwidth is approximately equal to two times of symbol rate", but I don't quite understand what the null-to-null bandwidth is here, and how it equals two times of symbol rate (although I read that null-to-null bandwidth for PSK is $ 2/T $ but I don't know what it comes from).

With no filtering or pulse shaping what you have is your signal with the default rectangular pulse shape
$$
\Pi(t) = \begin{cases}1 & \text{for} & -\frac T2\leq t\leq \frac T2\\0 & \text{otherwise}\end{cases}
$$
Where $T$ is the symbol duration. The Fourier transform of the such a pulse result in a sinc function of the form
$$
H(f) = T\operatorname{sinc}\big((f-f_c)T\big) \equiv T\frac{\sin\big(\pi (f - f_c)T\big)}{\pi \left(f-f_c\right)T}\tag{1}
$$
The zeros in Equation $(1)$ are at
\begin{align}
(f - f_c)T &= k\quad \big\vert \quad k\in \mathbb{Z}_{\neq 0}\tag{with $f\neq f_c$}\\
\implies 2(f - f_c) &= 2kR_s\tag{2}
\end{align}
The sinc function in Equation $(1)$ dies slowly resulting in infinite frequency content. Hence the choice of other practical pulse shaping filters such as the SRRC for bandwidth conservation (and ISI mitigation). In brief, having rectangular pulses one has remnant lobes outside the bandwidth of interest whilst with practical ones such as the SRRC you get high attenuation outside the main bandwidth.
For bandpass modulated signals (e.g. PSK) the required double-sided bandwidth is defined as
$$
W_{\rm DSB} = (1 + \alpha)R_s\quad,\qquad 0 \leq \alpha \leq 1
$$
Where $\alpha$ is the roll-off of the raised-cosine pulse, a measure of excess bandwidth, and $R_s$ the symbol rate.
For bandpass systems, the null-to-null bandwidth is defined as
$$f_2 − f_1\tag{3}$$
where $f_2$ and $f_1$ are the first nulls of the magnitude spectrum $|H(f)|$ above and below $\arg\max_\limits f\left(|H(f)|\right)$ respectively. So, from Equation $(3)$ and $(2)$ you have that
\begin{align}
\arg\max_\limits f\left(|H(f)|\right) &= f_c\\
\left(f_2 - f_1\right) = 2(f - f_c) &\equiv 2R_s\tag{with $k = 1$}\\
\end{align}


Smoothing data by using Kalman filter I would like to ask about smoothing data by using Kalman filter. Due to quantization, I have data that is not smooth. How can I smooth this data by using Kalman Filter. For your information, the data is not constant and have some gradient.

Using the same state transition information as this answer to another question, but using:

y(t) = round(H*x_truth(:,t) + rand(1,1,"normal")*sqrt(R));

as the signal model's output equation, we can apply the same Kalman filter.

This is not really accurate, because the round function is a nonlinearity sort of like quantization. However, quantization can also be modeled as an additive noise, so we'll proceed.

The results are shown in the plot below.

Here,

the black line is the true position,
the red + signs are the quantized, noisy position measurements, and
the green line is the Kalman filter's estimate of the position.

Is this the sort of "smoothing" you're interested in?

And the error between the true position and the Kalman filter's estimate is:

Including full scilab script for reference.

// Signal Model
DeltaT = 0.1;
F = [1 DeltaT; 0 1];
G = [DeltaT^2/2; DeltaT];
H = [1 0];

x0 = [0;0];
sigma_a = 0.1;

Q = sigma_a^2;
R = 0.1;

N = 1000;

a = rand(1,N,"normal")*sigma_a;

x_truth(:,1) = x0;
for t=1:N,
    x_truth(:,t+1) = F*x_truth(:,t) + G*a(t);
    y(t) = round(H*x_truth(:,t) + rand(1,1,"normal")*sqrt(R));  // <== Changed line!!
end

// Kalman Filter
p0 = 100*eye(2,2);

xx(:,1) = x0;
pp = p0;
pp_norm(1) = norm(pp);
for t=1:N,
    [x1,p1,x,p] = kalm(y(t),xx(:,t),pp,F,G,H,Q,R);
    xx(:,t+1) = x1;
    pp = p1;
    pp_norm(t+1) = norm(pp);
end

// Plots
figure(1);
clf;
plot(x_truth(1,:),'ko');
plot(y,'r+');
plot(xx(1,:),'g.');

figure(2);
clf;
plot(x_truth(1,:) - xx(1,:));


An Intuitive Way to Understand the Sharpen Convolution Filter in Image Processing Can someone explain sharpen in a non-technical way. So, I don't mean that you need to multiply every pixel with 5 and subtract it with the pixel left, under, above and right.

Consider this filter in 1D (for the sake of simplicity):

[-1 1 0]

This filter responds to pixel differences, i.e. edges (left pixel * (-1) + center pixel * 1 + right pixel * 0) which is (center pixel - left pixel).

Response of this filter is positive on rising edges and negative on falling edges.

The filter is not symmetric and hence it will "shift" the image to the left. To make this filter responding to "edge on the left" and "edge on the right" equally, we have to compensate it by similar filter:

[-1 1 0] + [0 1 -1] = [-1 2 -1]

Now this can produce "edge image". We will enhance the edges by adding them to the original image:

[-1 2 -1] + [0 1 0] = [-1 3 -1]

This kernel have sum of 1 and does not produce negative values. It is the desired sharpening kernel.

Summary: Sharpening =  edge detection + original image.


Determining Power of White Noise Component of an Incomplete Time Signal I want to determine the power of the noise component of a time dependend signal x(t) (similar to picture below)

With constant sampling rate I use Parceval's Therorem and calculate the standard deviation with

$$
\sigma^2 = \frac{\text{abs}(\text{FFT}(x))^2}{N}
$$

where N is number of data points and fft matlab dft implementation. For the standard deviation according to white noise I remove the discrete frequencies of the spectrum and use above formula with the mean value of the remaining spectrum.

In some signals there are quite a few NaN data points (~ 10%), so said method can not be applied anymore.

Is there a mathematically correct method to determine the power of the noise component of the signal? Interpolation does not seem to be a good solution, because this would lead to an underestimation of the noise.

I assume you have information about the signal itself (Not the noise).

What you can do is use LS to estimate the signal (Easily compensate for the missing samples in the Model Matrix).
Once you have it estimated you can do 2 things (Which are the same):

Recreate the samples using the model. Subtract the measurements. You'd be left with IID Samples of the noise -> Use the empirical STD / Variance estimator.
If it is Linear Least Squares use its internal estimation of the noise Variance.

Good Luck!


Faster way to find the quantization level of a sample Given a sequence of real numbers that represent a sound signal i.e [1.8, 2.2, 2.2, 1.9, -1.5, -0.7], we must quantize this sequence by dividing the range [-4, 4] to 32 equal parts.

If we had to do this with the traditional way we would need to write 32 different ranges

[-4   , -3.75]

[-3.75, 3.5]

.
.
.

and check in which one the real value is in.

Is there an another way to achieve this? Maybe with a normalization function or what?

As hinted by Florian, to quantise any $x[n] \in \mathbb{R}$ to a $y[n] \in \mathbb{Z}$ with $K$ levels you can apply the function:

$$y[n] = round\left( \frac{x[n]-min(x)}{range(x)} (K-1)  \right)$$

Where $n \in [0,1,2..|x|]$ and $|.|$ denotes the length of the sequence.

This will map your $x$ from its $min(x)..max(x)$ range to a $[0..K)$ discrete steps range.

For example (in Octave):

Fs = 100; #Sampling frequency (in Hz)
f = 4; # Frequency of a simple sinusoidal signal (in Hz)
T = 1; # Timeframe length (in seconds)
K = 5; # Number of levels to quantise to
t = 0:(1./Fs):(T-(1./Fs));
p = 2.0 * pi * t;
x = sin(f*p);
y = round(((x-min(x))/2.0) * (K-1));
plot(y);
xlabel("Discrete time (sample)");
ylabel("Amplitude");
grid on;

Which produces:

You might also like to see uencode

Hope this helps.


Glitches in certain frequencies I am trying to produce some simple audio waves using WASAPI (Windows Audio Session API). If I skip all the abstractions, my program basically receives a buffer from the API, then writes results of a sine function (or some other simple function) of a selected frequency to the buffer sample by sample, while counting the samples to keep track of time. When it goes to the next buffer it continues with an offset of the samples passed. So in theory, it should be a simple quantization of the sine wave.

The problem is, that some frequencies work and some don't. In particular, multiples of 50 work perfectly fine, while all the other frequencies sound glitchy.

Any ideas on what problem this could be related to? Is it some quantization error maybe? Should I look into some concepts like signal-to-noise ratio? Or are these problems a matter of the DAC and not my concern.

Edit:

I agree with the comments bellow, that my code being wrong is more probable, but after trying to inspect every detail of the code and also noticing this strange behavior I just thought I might ask some more experienced people, and hope that this behavior might be suspicious of some more general problem. If it sounds like a problem specific to my implementation, then I'll try to post a simplified example of my code. (Only if this community is a right place to post such questions of course.)

I resolved my issue. The problem was with timing. My samples were rendered correctly, but the problem was with waiting for the buffer. I switched to event driven approach instead of manually calculating the buffer length and waiting, and the problem was resolved.

I don't know whether this question might be useful for someone. If not or if the source of the problem is not related to DSP, I will delete it. But for now I'm keeping it here with my answer.


midpoint smoothing spline Suppose we have a time series which have peaks and troughs.(The red curve below)

I would like to get an algorithm which is able to identify the peaks and troughs locations, then find the midpoint between every peak and trough in both X,Y coordinates (the orange points below)

Then connect those midpoints by spline. (the blue line below)

How can I accomplish it using matlab or so? or there is a common solution for such problem?

Thanks in advance,

This is not in matlab but in R.  It seems to do what you need.  The plot below mimics your original in colors (red is the original; circles are peaks and troughs; orange are the interpolating points; and blue is the interpolating signal).

R Code Below

#31438

#install.packages("quantmod")
require(signal)
require(quantmod)
T <- 100
filt <- butter(5,0.1)
red_line <- filter(filt, rnorm(T, 0, 1))
peaks <- findPeaks(red_line)
valleys <- findValleys(red_line)

plot(as.numeric(red_line), col="red", type="l")
points(peaks + 1, red_line[peaks], col="red")
points(valleys + 1, red_line[valleys], col="red", lwd = 3)

peak_idx <- 1
valley_idx <- 1
counter <- 1

values <- array(0, c(2,T))

while ( (peak_idx <= length(peaks)) && (valley_idx <= length(valleys)))
{
  mid_idx <- floor((peaks[peak_idx] + valleys[valley_idx]) / 2 )
  mid_value <- (red_line[peaks[peak_idx]] + red_line[valleys[valley_idx]])/2
  points(mid_idx, mid_value, col="orange", lwd = 5)

  values[1,counter] <- mid_idx
  values[2,counter] <- mid_value
  counter <- counter + 1

  if (valleys[valley_idx] > peaks[peak_idx])
  {
    peak_idx <- peak_idx + 1
  }
  else
  {
    valley_idx <- valley_idx + 1
  }
}

values[1,counter] <- T
values[2,counter] <- red_line[T]

interpolated <- spline(values[1,], values[2,], n=T)
lines(interpolated, col="blue")


Difference between Nyquist rate and Nyquist frequency? So I've been searched online and can't seem to find a clear cut answer to this question.

From my understanding, the Nyquist rate is double of the maximum frequency of a signal which Nyquist frequency is half of the Nyquist rate.

Which would conclude that Nyquist rate is the lower bound of sampling and Nyquist frequency would be the upper bound where Nyquist rate is where you can preserve the original signal and anything in between the Nyquist rate to Nyquist frequency would cause some aliasing.

Is this correct?

Harry Nyquist invented/discovered/proved a lot of things; it can be hard to keep track of them all. The three most important for signal processing and communications are probably these:

If you sample a (real) signal $s(t)$ at $f_s>2B$ samples per second, then $s(t)$ can be reconstructed from its samples, where $B$ is the bandwidth of $s(t)$. The lower bound $2B$ is often called the Nyquist rate.
If you sample a signal $s(t)$ at $f_s$ samples per second, then $f_N=f_s/2$ is frequently called the Nyquist frequency (also folding frequency). All frequencies in $s(t)$ higher than $f_N$ will be aliased back to a frequency between 0 and $f_N$.
Assume you have a baseband channel of bandwidth $B$ (that is, signals with bandwidth $B$ or less are not distorted by the channel). This channel allows transmission of a maximum of $2B$ orthogonal pulses per second. $2B$ in this context is also called the Nyquist rate.


Understanding Magnitude Spectrum of Images I am facing problem in reading Fourier domain of a given image. I don't understand what to interpret from it. For instance, consider this image.

Ok, so there is a dot in the middle, with some striking vertical and horizontal lines and there are also some circulating white lines which are low in intensity. What do I interpret from them? I basically converted the image to a monochrome image and then took the fourier transform of the monochrome image and displayed it along with original RGB image. Please help me to understand this. I tried every other resource but nothing got into my head. What are the uses of this interpretation? Also, in one of the answers here on DSP stackexchange, the images of sine and cosine were taken and then converted to spectrum which resulted white dots in it. I don't want to know that. I know that these dot represents the frequency component present in the sine and cosine. I want to understand everything with respect to my image only.

Thanks in advance!!

Many people face the same issues. Your image is far from showing stationarity: quite edgy, with circular features, so the Fourier analysis might not be so enlightening, without showing the unwrapped phase, which might contain many interesting features for this class of images. However:

so there is a dot in the middle: standard images have integer values between $0$ and $255$, hence the average is not zero. The average is the center dot. Remove it, and the global contrast may be better
some striking vertical and horizontal lines: FFT generally uses periodic extension, and the differences between the bottom and top, and right and left, might obfuscate details.

In the following pictures, you get space data and the Fourier spectrum. The first line is your image, with an horizontal/vertical cross. When you window the image (second line), those disappear, and the most luminous  decreasing slope is, probably, due to the index finger, and the parallel lines from the two biggest gears. On the third row, the windowed image is subtracted from its  median filtered version.  One begins to see circles, related to locally oriented teeth.

So, to interpret complicated data, remove first order artifacts, and slow background features, to better interpret details.


Confusion between PSK and QPSK What is the difference between PSK (Phase Shift Keying) and QPSK (Quadrature Phase Shift Keying)? Doesn't PSK need quadrature component for its working?

QPSK modulation is a PSK modulation with exactly 4 constellation points.  The constellation points for QPSK are normally at $\left(\pm\dfrac{1}{\sqrt{2}}\pm\dfrac{1}{\sqrt{2}}j\right) $


spectral leakage for integral number of periods I'm attempting to understand spectral leakage and read the following quote on wikipedia:

  Windowing a sinusoid causes spectral leakage, even if the sinusoid has an integer number of cycles within a rectangular window.

This confuses me, since I thought using an integral number of periods caused the leakage to disappear. Moreover, my attempts to quantify/simulate this corroborate my beliefs. For example, take the following sampled, 1Hz simple sinusoid.

import numpy as np

f = 1
ts = f / 3
fs = 1 / ts
t = np.arange(0, 5.1 * f, step=ts)
s = np.sin(2*np.pi*f*t)

with fft:

fft = np.absolute(np.fft.rfft(s))
fft /= np.sum(fft)
bins = np.fft.rfftfreq(len(t))

Now if I add an additional sample,

import numpy as np

f = 1
ts = f / 3
fs = 1 / ts
t = np.arange(0, 5.1 * f, step=ts)
s = np.sin(2*np.pi*f*t)

As expected, I get spectral leakage.

Is the Wikipedia quote incorrect, or did I misinterpret it?

Note: I think a better tag for this is spectral-leakage, but I don't have enough reputation to add it.

If you have an integer number of cycles (exactly) within a DFT frame, and the signal is not windowed (aside from the rectangular window extending to the edges of the DFT fram that would result), there is no spectral leakage. What the article is referring to I believe is more commonly called zero padding which would then cause the DFT to approach the DTFT, interpolating frequency samples between the bins. I personally wouldn't call this spectral leakage although you could mathematically get to the same conclusion in that energy in a DFT bin simply represents correlation to the input signal. If energy is there, and the bin isn't your true frequency, then sure, it is "spectral leakage". In that paradigm however of zero-padding, I prefer to associate it with the concept of frequency interpolation (again the math and the results are the same, so no disagreement - just convention).

Spectral leakage occurs when you have a non-integer number of cycles such that your true frequency is mid way between DFT bins. Each DFT bin actually has a very wide frequency response with nulls at every other bin. So a frequency that is in between bins will show up within the response of the other bins, while a frequency that is exactly on bin center (integer number of cycles) will be in the null of the responses of the other bins. There is therefore no leakage under this condition (that is not an illusion as the article suggests- there is simply no leakage - but I do see how the author presented both options and stated "depending on your point of view"). And for the other conditions the signal "leaks" into the response of the other bins.

I explain this further here at this post.

Intuition for sidelobes in FFT


Load and Play Wave File (`.wav`) on MATLAB I created a wav file by "audacity" app on my mac and imported in MATLAB. This is the simple code I wrote to read the file. But as it reads, it plays so much slower than its actual speed. Does anybody know how to solve this problem? Thanks.

d = audioread('HelloWorld.wav');
sound(d)

When you load Audio file to MATLAB you basically load array of numbers.
When you play Audio there is a parameter which is needed - Sampling Rate.
This is what tells the player how fast to play the data.

Think if you load a file with 1000 samples. How long should is the audio?
Unless you supply some information, the player can not know (MATLAB just use a default which doesn't match your data).

If you look at the documentation of the sound() function you'd see there is an option to use it as sound(audioSignal, samplingRate).

In order to get the Sampling Rate of a Wave file use: [audioSignal, samplingRate] = audioread(filePath). Read more at audioread().


Properties of a Kalman filter of a non controllable system Let's say I have a standard system
\begin{align}
x(t+1)&=Ax\\
y(t) &=Cx(t)
\end{align}

As you can see $B=0$, so the system is not controllable.

For the steady state Kalman filter I'd say that if $A$ is Hurwitz then it converges.

But what if $A$ is say,  $\mathbf I_n$ (i.e. eye(n)), what can I say about the optimality of the SS.  filter varying $C$?

In this setting you are considering an observation problem. Controllability does not play a role. In geenral, the condition under which you can obtain a convergent estimate is that the couple $(A,C)$ is detectable, i.e. that you can design a gain matrix $L$ so that $A+LC$ is Hurwitz. Under this hypothesis the eigenvalues of $A+LC$ will be with negative real-part and they are assigned in a optimal way, in the kalman sense.
The kalman gains depend on $C$, therefore as you vary $C$ you should vary the gains accordingly to obtain optimality


Phase response of moving average filter -- how to interpret? There are many articles on the frequency response of the moving average filter but they all seem to focus on magnitude. However the phase response is intriguing and I find it hard to interpret. The phase appears to "wrap", but it wraps within the [-pi, pi) interval rather than at its edges. Example below:

A phase unwrapping algorithm would not resolve this, so it is really a "pseudo-wrap". Further, if I add taps to the moving average, it flattens this process out, so I suspect that mathematically, the moving average filter will never reach 0 or 2 pi, though I have never seen an explanation why. Example of an 11-tap:

I find this behaviour fascinating and would be interested in the interpretation of an expert. Does this suggest that features will be distorted at certain "weak spots" in the frequency response? Is it correct to call the phase of a moving average filter "piecewise-linear" rather than linear? I suspect not, given that symmetric FIR filters are analytically shown to have linear phase, but I have a hard time calling this linear.

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

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

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

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

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

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

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


The Gradient of Least Squares of 2D Image Convolution Given the objective function:

$$ \frac{1}{2} {\left\| h \ast x - y \right\|}_{2}^{2} $$

Where $ h $ is the 2D convolution kernel and $ x $ is the 2D convolution image and $ y $ is a given 2D image.

What would be:

$$ \frac{\mathrm{d} \frac{1}{2} {\left\| h \ast x - y \right\|}_{2}^{2} }{\mathrm{d} h}, \; \frac{\mathrm{d} \frac{1}{2} {\left\| h \ast x - y \right\|}_{2}^{2} }{\mathrm{d} x} $$

The easiest approach would be writing each case using Matrix Form of the convolution.
In this answer we assume the discrete convolution is applied only on valid support (Matching MATLAB's valid parameter for the convolution).
Namely, given $ x \in \mathbb{R}^{m \times n} $ and $ h \in \mathbb{R}^{k \times l} $ then $ h \ast x \in \mathbb{R}^{ \left( m - k + 1 \right) \times \left( n - l + 1 \right) } $. Needless to say $ m \geq k $ and $ n \geq l $ as otherwise the operation isn't well defined. Pay attention that this form of convolution isn't commutative.

Remark
For full convolution the problem is easier as the operation is commutative hence no difference between the gradient with respect to $ h $ or $ x $ as one could switch them.

Gradient with Respect to Convolution Kernel $ h $

The matrix form is given by:

$$ f \left( h \right) = \frac{1}{2} {\left\| X h - y \right\|}_{2}^{2} $$

Where $ X $ is the 2D Convolution Matrix Form of the image. Then:

$$ \frac{\mathrm{d} f \left( h \right) }{\mathrm{d} h} = {X}^{T} \left( X h - y \right) $$

The $ {X}^{T} $ forms a correlation (Versus Convolution) with full support of the operation (Equivalent of the full convolution shape in MATLAB syntax).

Hence we have:

$$ \frac{\mathrm{d} \frac{1}{2} {\left\| h \ast x - y \right\|}_{2}^{2} }{\mathrm{d} h} = x \star \left( h \ast x - y \right) $$

In MATLAB Code:

conv2(mX(end:-1:1, end:-1:1), (conv2(mX, mH, CONVOLUTION_MODE_VALID) - mY), CONVOLUTION_MODE_VALID);

Gradient with Respect to Convolution Image $ x $

The matrix form is given by:

$$ f \left( x \right) = \frac{1}{2} {\left\| H x - y \right\|}_{2}^{2} $$

Where $ H $ is the 2D Convolution Matrix Form of the kernel. Then:

$$ \frac{\mathrm{d} f \left( x \right) }{\mathrm{d} x} = {H}^{T} \left( H x - y \right) $$

In MATLAB Code:

conv2((conv2(mX, mH, CONVOLUTION_MODE_VALID) - mY), mH(end:-1:1, end:-1:1), CONVOLUTION_MODE_FULL);

MATLAB Code

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


Bandpass filter with very high sampling rate Consider a bandpass filter with Low Cut 17Hz, High Cut 22 Hz, Fs = 45000 Hz and Order = 6. When I pass a mixture of multiple sinusoidal waves through this filter (with a sine wave of frequency 20 Hz), I noticed that Scipy's Butter Bandpass filter, upon filtering, returned values that were unexpectedly high. I did not get a sinusoidal output, but rather, ended up getting an almost straight line that was scaled up by a factor of 1e305.  Since I haven't done a course on DSP, and hail from a CS Background, any help would be greatly appreciated.

Here is my code:

from scipy.signal import butter, filtfilt

def butter_bandpass(lowcut, highcut, fs, order=5):
    nyq = 0.5 * fs
    low = lowcut / nyq
    high = highcut / nyq
    b, a = butter(order, [low, high], btype='band')
    return b, a

def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):
    b, a = butter_bandpass(lowcut, highcut, fs, order=order)
    y = filtfilt(b, a, data)
    return y

Limited numerical precision makes your filter unstable.

In this case you can fix it by implementing the filter in second order sections using the zeros and the poles of the filter.


Length of Hamming Window to Resolve a Harmonic I am trying to answer the question:

"Given a single-sourced pitch signal with a fundamental frequency of 440Hz and a sampling rate of 44100Hz how long should the Hamming window of the STFT analysis be to resolve the harmonics of the signal"

My understanding of resolving a harmonic is to ensure the frequency is represented in the frequency domain e.g. the top graph plots the sum of 4 sinusoids in the time domain. The bottom graph shows 4 distinct frequency peaks in the frequency domain (4 harmonics are resolved).

For the question, I use the below code to calculate the STFT of a sinusoid with 440hz frequency.
Trying length 10 vs 100, based on the logic above I would say N=10 does not resolve the harmonic while N=100 does.

Is this correct? Is there a mathematical intuition for calculating the length that will resolve a harmonic?

Thanks.

import numpy as np
from matplotlib import pyplot as plt

SR = 44100 #Sampling rate (sometimes denoted as fs)
T = 1/SR #Time step in secs
t = np.arange(0,1,T) #A 1 second array arrange intp 44100 samples

##Sinusoid parameters
f = 400 #Frequency
A = 1 #Amplitude
x = A*np.sin(2*np.pi*f*t) #sinusoid equation

N = 100 #Hamming window of length N (this is what we want to optimize)
win = np.hamming(N)
ind = np.arange(0,N,1)
Window_sig = x[ind]*win #multiplying our first N samples by our window

Xwin = np.fft.fft(windowSig)
NposFreqs = round(N/2)+1
freqHz = ind[0:NposFreqs]*SR/N

plt.figure()
plt.plot(freqHz,np.abs(Xwin[0:NposFreqs]))
plt.xlabel('Frequency [Hz]')
plt.ylabel('Magnitude')
plt.title('Signal magnitude specturm, N = '+str(N))
plt.show()

Around 1000 samples to resolve the psychoacoustic pitch of concert A.

The fundamental period of 440 Hz is just over 100 samples at a sample rate of 44100.  But humans typically need somewhere around the range 6 to 10 pitch periods to hear a pitched sound instead of noise or clicks.  Since FFT windows attenuate the edges of the waveform, using the larger number of periods (10) is likely more reliable in matching human perception.

So I'd use a 1024 length FFT for initial testing.

To inspect higher harmonics if you already know the exact pitch frequency, you could use a shorter window, at least 2 pitch periods for a Hamming window, and zero pad after windowing to use a much longer DFT that is an exact integer multiple (say 16X) of the pitch period in length.


Digital Equations problem - The fundamental ones that everyone should know I think I am looking at conflicting equations from a few sources, or maybe I just dont understand it. But please can someone help me with these fundamental equations that everyone should know

$  \Omega $ is Discrete frequency
$  N $ is Discrete Period in number of samples
$  \omega $ is Continuous frequency
$  T $ is Continuous time period
$  F_{s} $ is Sampling frequency
$  T_{s} = n $ is Sampling Time Instant
$  k $ is Discrete harmonic index which goes from 0 to N/2 for a real signal

First equation forms I have seen
$$
\Omega = 2\pi f/F_{s} =\omega / F_{s} = \omega T_{s} = \omega n = 2\pi n / T
$$
Second equation I have seen
$$
\Omega N = 2\pi k
$$
$$
\Omega  = 2\pi k / N
$$
Third equation I have seen
$$
\Omega_{k} = k F_{s} / N
$$

First and second look the same with if $k=T_{s} $ and if $N=T$ and first and third seem contrary.

Can someone explain please?
Also if you can point in in a direction where all this is neatly derived and defined?

Thanks

D ear N atalie, as you also know that these are just letters of English (org Greek) alphabet which by themselves have no meaning or explanations other than what you have arbitrarily imposed on them and the following is by convention what's being imposed on them in the mathematics, physics and the standard DSP literature as accepted unless otherwise stated explicitly...

$\omega$ can represent both continuous-time and discrete-time frequencies in radians-per-second and radians-per-sample respectively.
When continuous-time and discrete-time signals are used simultaneously (such as in sampling operations) then $\Omega$ represents continuous-time radians per second frequency and $\omega$ represents the discrete-time radians per sample frequency.
$\Omega = 2 \pi f$ interprets $f$ as the continuous-time frequency in Hertz.
$\omega = 2 \pi f$ interprets $f$ the same (as frequency in Hertz) if $\omega$ was used for continuous-time radians per second frequency, otherwise if $\omega$ is used for discrete-time radians per sample frequency, then $f$ has no such interpretation.
If a discrete-time signal $x[n]$ is periodic, then its (fundamental) period is indicated by the capital letter $N$ indicating a minimum integer larger than $1$ for which $x[n] = x[n+N]$ holds.
if a continuous-time signal $x(t)$ is periodic, then its (fundamental) period is indicated by the capital letter $T_0$ indicating a minimum real number larger than $0$ for which $x(t) = x(t+T_0)$ holds. This $T_0$ can also be taken as $T$.
When a continuous-time signal is uniformly sampled by an (ideal) periodic impulse train, then the period of sampling is indicated by the real number $T_s$ seconds per sample.
The reciprocal of the sampling period $T_s$ is denoted as the sampling rate: $F_s = 1/T_s$ in samples per second. If $F_s$ considered as the sampling frequency, then $F_s$ is interpreted in Hz.
Letter $k$ is mostly used as the frequency integer index for the DFT or DFS $X[k]$ of a discrete-time sequence $x[n]$. It's also used to indicate the continuous-time Fourier series coefficient indexing. It's also used to indicate convolution sums index in discrete-time.


Help Understanding the Robert Cross Operator I have been reading this Wikipedia article:  https://en.wikipedia.org/wiki/Roberts_cross
And it says this:

Now I have also been reading this for the Robert Cross operator (from Feature Extraction and Image Processing 3rd Edition):

I can understand The books explanation (kinda) but not wikipedia.
So in Code I will have something like this:
EdgeValue(at x,y) = Max[(Px+1,y - Px,y+1), (Px,y - Px+1,y+1)]

So from the 2 values I choose the max. And in the new image I just draw that value in x,y. Makes sense.
What I do not get is Wikipedia. Why does it go on finding the gradient?
Also that looks like the Pythagorean Theorem not a gradient equation.

What I do not get is Wikipedia. Why does it go on finding the gradient?

You subtract two nearby pixel values. That is a gradient, or as close as you can get on a discrete image.

  Also that looks like the Pythagorean Theorem not a gradient equation.

As the book text say, "it is often better to consider the two templates as providing the components of an edge vector". If you want the length of a vector, you use the Pythagorean Theorem.

The better question is: why use the max or the sum of the absolute values, or: what's the difference?

To show that, it's best to try the filter on a rotationally symmetric template:

that way, we can see how it acts on different edge angles:

As you can see, adding the absolute values of the gradient vector (this is sometimes called the "L1-Norm") gives larger edge strengths for horizontal and vertical edges than diagonal edges. Taking the max of the two absolute values of the gradient vector (also known as the "L-infinity Norm") gives slightly larger edge values for diagonal edges. Using the ordinary Euclidean norm (i.e. the Pythagorean Theorem) gives exactly the same edge strength for any angle, which is of course the right result.

The only downside (I'm aware of) is that this requires two multiplications and a square root for every pixel, which used to be very expensive operations when the Roberts operator was invented. They're still expensive on mobile or embedded platforms today, but on a PC, you might not even notice the difference.


Kalman Filter to estimate 3D position of a node Code given on this link works for 1D:
Kalman filter for position and velocity: introducing speed estimates

In my problem I need to estimate 3D position.What is the criteria ?

How F, G ,H,Q and R change in 3D case.
This problem is restricted to estimate position only  (No velocity,No acceleration).What measurement vector should contain.Only positions?
Please respond!I am sorry for any inconvenience

As I said in the comments, just follow the 1D case from Wikipedia and augment it with the extra $y$ and $z$ dimensions (and velocities):
$$
\mathbf{x}_k = \left [
\begin{array}{c}
x\\ \dot{x}\\ y\\ \dot{y}\\ z\\ \dot{z}
\end{array}
\right]
$$

You will also need to augment $\mathbf{F}$ and $\mathbf{G}$:
$$
\mathbf{F} = \left[ \begin{array}{cccccc}
1 & \Delta t & 0 & 0 & 0 & 0\\
0 &        1 & 0 & 0 & 0 & 0\\
0 & 0 & 1 & \Delta t & 0 &  0\\
0 & 0 & 0 &        1 & 0 & 0\\
0 & 0 & 0 & 0 & 1 & \Delta t \\
0 & 0 & 0 & 0 & 0 &        1 \\
\end{array}
\right]
$$
$$
\mathbf{G} = \left[ \begin{array}{c}
\frac{\Delta t^2}{2}\\
\Delta t\\
\frac{\Delta t^2}{2}\\
\Delta t\\
\frac{\Delta t^2}{2}\\
\Delta t\\
\end{array}
\right]
$$

And then $\mathbf{H}$ is just
$$
\mathbf{H} = \left[ \begin{array}{cccccc}
1 & 0 & 0 & 0 & 0 & 0\\
0 & 0 & 1 & 0 & 0 & 0\\
0 & 0 & 0 & 0 & 1 & 0\\
\end{array}
\right]
$$
so that $\mathbf{z}_k$ is
$$
\mathbf{z}_k =  \left[ \begin{array}{c}
x\\
y\\
z
\end{array}
\right]
$$


Is this an error in Oppenheim and Schafer's Discrete-Time Signal Processing? In Discrete-Time Signal Processing by Alan V. Oppenheim and Ronald W. Schafer (3rd Ed.), in Figure 4.47 the input of D/A converter is $\hat{y}[n]$ but later in Figure 4.64 the input of D/A converter is $\hat{x}[n]$. Is this a mistake?
Normally, based on Figure 4.47 $\hat{y}[n]$ is the output of the discrete-time system with input $\hat{x}[n]$.

It is not a mistake. In section 4.8.4 the authors introduce the reconstruction of $x_r(t)$ from $x[n]$ in equation $(4.140)$; and in the paragraph below, the authors consider the input sequence and output signal to the ideal D/C converter as $x[n]$ and $x_r(t)$ respectively, as shown below
$$x[n]{\longrightarrow}\boxed{\textit{Ideal D/C converter}}{\longrightarrow}x_r(t)$$
Which will be more in line with the last input-output block of Figure 4.47(a) shown as
$$\cdots \quad y[n]{\longrightarrow}\boxed{\textit{D/C}}{\longrightarrow}y_r(t)$$
To clarify, the extension from ideal D/C converter to the D/A and their corresponding maths with inputs-outputs relations and the figures, please consider reading the full section 4.8.4 with emphasis on the change of notation just after equation $(4.140)$.


Finding System transfer function, $H(z)$ or, equivalently, its impulse response, $h[n]$ I just started on DSP and I have a question that I would like to ask.

I have a zero-mean uncorrelated wide sense stationary discrete-time random process {$x[n]$, $n$ is a set of integers}. $x[n]$ is input to a causal linear time-invariant system below, denoted by its transfer function $H(z)$.

$$x[n]{\longrightarrow}\boxed{H(z)}{\longrightarrow }y[n]$$

Since the output is a random process, suppose that its autocorrelation function, $R_y[k]$, is measured without approximation error.

From the information given, is there a method to identify the system transfer function, $H(z)$ or, equivalently, its impulse response, $h[n]$?

From only the process statistics, you cannot derive the exact impulse response/transfer function. However, by the Wiener-Khinchin-Theorem you know, that the power spectrum density (PSD) of an input signal to an LTI system is transformed with the magnitude-square of the transfer function:

Let $S_x(\omega)$ be the PSD of the input and $S_y(\omega)$ be the PSD at the output of the system. Note that

$$\begin{align}
 S_x(\omega) &= \mathrm{DTFT} \{ R_x[n] \} \\
             &= |X(e^{j\omega})|^2 \\
\end{align}$$

where $X(e^{j\omega})= \mathrm{DTFT} \{x[n]\}$ and similarly for $S_y(\omega)$.

Then, you have

$$
S_y(\omega)=|H(e^{j\omega})|^2 S_x(\omega)
$$

and you can at least calculate the magnitude of $H(e^{j\omega})$ by inverting the equation.


Difference between $\tt freqz()$ and $\tt plot(abs(fft()))$ in MATLAB I've plotted the response of a high-pass filter x=[1 -1];

freqz(x) gave me

where as the magnitude plot for the fft gave me

Why are the two plots different? Both are frequency responses, then why are my plots are diffferent. What does freqz() do?

You calculating FFT only from two samples. You need to pad your impulse response with zeros to get a valid result. So in MATLAB that would be:

N = 1024; % Number of points to evaluate at
% Create the vector of angular frequencies at one more point.
% Filter itself
b=[1,-1];
[h_f, w_f] = freqz(b, 1);
figure
grid on
hold on
plot(w_f, abs(h_f), 'or') % MATLAB
h = [b, zeros(1,N-2)];
HH = abs(fft(h));
HH = HH(1:length(w_f));
plot(w_f, HH); % Manual calculation
legend({'MATLAB freqz', 'Manual'})

As you can see it matches first and last value from fft you calculated. Please keep in mind that it is shown in linear - not dB scale.

For more info you can see my previous answer.


Does Zero Padding Distort the Spectrum of a Signal? It's said to "sample the DTFT", revealing what "DFT fails to see". And I fail to see how this sampling isn't distortion.
The "spectrum" aims to provide the sinusoidal frequencies comprising a signal. If we're getting a nonzero distribution from 1 to 63Hz for a pure 1Hz tone, that's by definition different.
Further, adding zeros changes the signal; who is to tell that these zeros didn't occur naturally as a part of it? If they did, the spectrum must look different - and it does. That it mathematically happens to be equivalent to sampling a related but different altogether transform is only a coincidence.
What's the deal?

DFT vs DTFT, Fourier Transform
The problem appears rooted in viewing DFT as a 'special case' of the continuous Fourier Transform, and of its input as some signal with legitimate frequency contents. That's a fallacy.
The DFT is a standalone mathematical transformation which does not condition itself upon anything but the input being finitely long, and having finite values. It's an information encoding algorithm. It operates on any and all finite information - be it signals, Shakesperean sonnets, or a cup of water - as long as represented numerically.
DFT can be shown to be equivalent to solving $N$ simultaneous equations defining an $N$ long input; it's a complete, lossless transform.

What the spectrum "should be" is not what we get

Should be: amounts of sinusoidal frequencies comprising the input data.$^1$
Is: amounts of sinusoidal frequencies that sum to the input data.$^1$

See here. The input can be a pulse, or a single nonzero value and the rest zeros; there's nothing of periodic nature about either of these, nor do they need to stem from a physical source which is sinusoidal. The only truth is that the coefficients we find give us the sinusoids that sum to the input; with another transform, it can be triangle waves.
1: for a real-valued input. For complex, something similar may hold, but it's still a sum of bases (discussion here).

DFT coefficients are not same as DTFT or FT coefficients
Both the forward and inverse transforms to acquire respective coefficients imply very different things for what the coefficients "mean". Under certain assumptions they can be treated equivalently, but not fundamentally.
This is clear from the inverse-transform perspective, or reconstruction; assume for now a rectangular-windowed DTFT:

DFT: sum of $N$ bases, each $N$-long, finite in value.
DFT coefficient: describes one such basis. $N$ coefficients in total.
DTFT: infinitesimal sum (integral) of infinitely many bases, each $N$-long, finite in value.
DTFT coefficient: describes one such basis. infinitely many coefficients in total.

From forward-transform:

DFT coefficient: obtained by adding (discretized complex phasor at integer frequency) * (input), pointwise products.
DTFT coefficient: obtained the same way. Except not really; see below.

Result? If you use one transform's coefficients in another transform, you get nonsense.

Generator vs Data
A generator is a complete descriptor - example: function, continuous. For describing a 1D signal, the signal is infinitely resolved; we can determine every possible intermediate value. From a generator we can get data, but not vice versa; even a truly continuous segment has infinite interpolations.
The DTFT (no windows) works with generators, not data. Its coefficients are described by a function (continuous), which are also determined by an operation between functions (one describing an infinitely long but discrete input, another the complex phasor at a continuous $\omega$).
The rectangular-windowed DTFT, on the other hand, works with both, but only "with permission" from generators - that is, the resulting spectrum happens to be the same whether we compute for $n=0$ to $N-1$ or $-\infty$ to $\infty$.
In contrast, the DFT works with data. No generators. But we can use generators if we wish - it's now the reverse case, generators "with permission" from data.

Zero-padding: DFT coeffs treated as DTFT coeffs -- all code.
Before the meat, last piece of background. Take 1Hz over two seconds; what's the DFT (which freq is nonzero)? 2. Now take 0.1Hz over 10 seconds; what's the DFT? 1. The DFT correlates against basis functions spanning the length of the input signal, with no regard to "actual time". $k=1$ corresponds to 1 cycle over the length of the signal.
We can use this to interpret "treat DFT coeffs as DTFT coeffs". Suppose
$$s[n] = \cos(2 \pi n / N),\ N = 128,\ n=[0, ..., N - 1].$$

Pre-pad, the spectrum is just what we expect: nonzero everywhere but at $k=1$ and $k=128$. Post-pad - a madhouse. What's the deal?
Recall how "0.1Hz over 10 seconds" $\rightarrow k=1$? By the same token, the original $f=1$ is now "seen" as $k=2$ - i.e., it correlates with the $k=2$ basis. Then what does $k=1$ "see"? Let's see:

(Only real component of basis shown for clarity; can add the imaginary component for completeness)
Put differently, the "normalized" frequency of the DFT is divided by two as a result of doubling the original signal's length with the padding, so $k_p=2 \leftrightarrow k = 1$, and $k_p=1 \leftrightarrow k = 0.5$. Or, integer frequencies with a padded signal's spectrum correspond to fractional frequencies in the unpadded's spectrum.
Is there more evidence? Yes; first, see this. If above is true, then we
can correlate the original signal with a $k=0.5$ basis directly, and yield the same coefficient as with $k_p=1$. And indeed:

Why one window for both padded and unpadded? ... because they look identical - check the code. Mean absolute difference: 6e-14 (within float).

In what sense do we "sample the DTFT"? What's the role of zero-padding?
Former should be clear now; the DTFT evaluates itself for every $\omega$. For any given $\omega$, this  exactly the same as taking DFT with an equivalent $k$ basis. In other words, DTFT correlates with neighboring fractional frequencies, just like DFT correlates with integer frequencies. Simpler, "DTFT is DFT with $k$ spanning continuously" (but not quite; see caveats in above sections).
So why does zero-padding "result in" the DFT sampling DTFT's spectrum? The answer is closer if we reformulate the question: Why does zero-padding result in the DFT correlating with fractional frequency bases?
Yes, they're indeed one and the same question, which was answered in the previous section. A DFT with $2N$ coefficients is equivalent to a DFT with $N$ coefficients at integer frequency bases, and another DFT with $N$ coefficients at fractional frequency bases$^1$.
But why zero-padding? Why not one-padding, or pi-padding? Because it's the only way for an $M>N$-DFT with integer frequencies to be equivalent to an $M$-DFT with fractional frequencies$^1$. Refer to the blue/orange figure; note how the sum of products of the orange with the blue for $k_p=2$ is the same as for $k=1$ with the unpadded blue (not shown; just imagine). If the zeros were anything else, the values might have been proportional, or somehow similar, but not equal. In technical terms, zero is orthogonal to all basis functions over any interval.
1: if $M=3N$, then the integer-DFT will correlate $k=[0, 1, ..., 3N-1]$ bases with the padded signal, and the "mixed" DFT will correlate $k=[0, 0.33, 0.66, 1, ..., N-1, N - .66, N - .33]$ with the original signal. Note that if $M$ is not an integer multiple of $N$, then the padded's integer-DFT may not have any integer equivalents with the unpadded's DFT.

Does zero-padding distort the spectrum?
The grand-prize. The answer to which is... yes and no:

No: zero padding is equivalent to correlating with neighboring, non-integer frequencies; there is no "distortion" here, the coefficients describe exactly what they intend to.
Yes: if "the spectrum" is supposed to describe the frequencies in the input signal (e.g. $f=1$), then correlating with bases neighboring $1$ yields non-zero coefficients, implying non-zero frequencies around $1$, whereas the signal has no such frequencies. Further, only integer-ratio frequencies are completely decorrelated (e.g. 1 and 2, 1.5 and 3), so we'll get many high and low non-zero frequencies, whereas the input has only one frequency.

So what's the deal? The deal is - DFT does NOT describe input's spectrum. This is the entire point of this answer, which was so warmly-received because of a technicality I was quick to retract with a NOTE: on top. Ignore the downvotes, they mean nothing without explanation.
Since it's not the spectrum we expect, it's meaningless to speak of a "distortion" as if it ever gave the actual spectrum. The input could be $f=1.1$, and then the unpadded DFT would be "distortion", and the only true spectrum would stem from incrementing the bases by $1.1$.
This also explains why we can't do much better than the DFT; we don't know what the "actual" frequency of a signal is without its generator (function). The DFT is simply our "best guess", one which must be interpreted carefully.

Possible to pad without distorting?
Yes, but in a specific sense: if the number of zeros added is an integer-multiple of the signal's length, then the DFT will contain the unpadded's spectrum. For example, pad by $2N$; then, the DFT, from original's point of view, is doing $k=[0,.33,.66,1,...,N-1,N-.66,N-.33]$. Notice how $[0, 1, ..., N-1]$ are all captured.
In other words, we can pad without losing the original spectrum.

Re: other answers

The fact that the input is discrete or finite has virtually nothing to do with this - which this upvoted answer contradicts. Though, a little more work is due to show this; short version is, imagine the "winding machine" winding forever, and why that yields an impulse, and how that relates to the finite spike in DFT. Whatever the answer, it's way besides the point.


What do we achieve by using PAM symbols symmetrical with respect to 0? I can't seem to find the answer to this one anywhere online. We want them symmetrical with respect to 0 and also at even distance between each other. Why is this important?

The probability of error for PAM depends on the distance between symbols. This means that the following two PAM constellations, both with distance $2$, will have the same probability of error:
\begin{align}
\mathcal{A}_1 &= \lbrace 2, 4 \rbrace \\
\mathcal{A}_2 &= \lbrace -1, 1 \rbrace
\end{align}

When using $\mathcal{A}_1$, the average energy spent per transmitted symbol is $(4+16)/2=10$ joules per symbol. When using $\mathcal{A}_2$, however, the average energy is $(1+1)/2=1$ joule per symbol: ten times less.

Even if transmitting a signal with constant DC (as pointed out by @AlexTP) were not a problem, you still want to use constellations that are symmetrical around zero because they minimize the energy expenditure for a given BER.

A good and simple excercise is to prove that this is true in general, and not just for the example I provided.


Why use PSD to analyze a signal instead of its Fourier Transform? Let's imagine I have a perfectly known signal $s(t)$ and I want to analyze its frequency components.
I noticed that mostly everyone would use the PSD of $s(t)$ to do that, instead of simply use the magnitude of the Fourier transform of $s(t)$. Is there a reason for that ?

The power spectral density (PSD) is a natural measure of the signal's power content with respect to frequency. A central part of non-parametric signal processing is to provide a "best" estimate  of the "true" PSD from knowing only one or some "realizations" with finite length. By taking into account the influence of stationary random processes, it should be noted that the Fourier transform is not defined for processes with infinite energy.
However, by looking at the second-order properties (autocorrelation), the true PSD can be well-defined, and the Fourier transform of  the realization  (on a finite horizon) can be used to provide an estimate of it.


Column operations in LDPC generator matrix I'm trying to produce a generator matrix from a starting low density parity check matrix. There are lots of references on the topic (including this Signal Processing Stack Overflow answer here). They make sense, with the exception of one problem I'm having.

My resource is the book Error Correction Coding: Mathematical Methods and Algorithms by Todd K. Moon. In it, there is a section that describes this that I will replicate below, and then highlight the part I'm stuck on.

  For a code specified by a parity check matrix $A$, it is expedient for encoding purposes
  to determine the corresponding generator matrix $G$. A systematic generator matrix may
  be found as follows. Using Gaussian elimination with column pivoting as necessary (with
  binary arithmetic) determine an $M$ x $M$ matrix $A^{-1}$ so that

  $$H = A_p^{-1} A = [I \quad A_2].$$

  (If such a matrix $A$, does not exist, then $A$ is rank deficient, $r = \mathrm{rank}(A) < M$. In this
  case, form $H$ by truncating the linearly dependent rows from $A_p^{-1} A$. The corresponding
  code has $R = K/N > (N - M)/N$, so it is a higher rate code than the dimensions of $A$
  would suggest.) Having found $H$, form

  $$G=\begin{bmatrix}A_2 \\ I\end{bmatrix}.$$

  Then $H G = 0$, so $A_p H G = A G = 0$, so $G$ is a generator matrix for $A$. While $A$ may be
  sparse (as discussed below), neither the systematic generator $G$ nor $H$ is necessarily sparse.

Here's the problem I'm running into:
I start with $A$ and do Gauss-Jordan elimination in $\mathbb{G_2}$ to get it into reduced row echelon form.

Sometimes the matrix I've been given (I'm not designing it, have to use what I'm given) can't be put into that form with solely row operations.

It requires actually swapping the order of some columns (is that what he means by "column-pivoting" in the reference above?). Everything in me screams that's wrong, but it turns out if I just do the same thing on the decoding side it works out.

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

Basically I have to re-do the generator matrix, which is by far the slowest part of my process (I'm going to make some improvements there though, using the QC-LDPC structure of the parity matrix).

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

Is there something I'm missing? For what it's worth I also have the example matlab code from the aforementioned book and this is what he does (column reordering and then re-applying that re-orderying on the decode side).

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

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

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

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

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

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

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

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

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


Is there any method to find K value in KNN? I am new to Matlab. For me the task is to perform data classification on some data set. The given data set consists of numeric values. One thing that bothers me is I have to perform all KNN, ANN, and SVM. In performing KNN, I faced a problem of how to find the value of K. If it is user defined? From which value we have to start and when to stop it?

Basically too large k may tend to underfit, while smaller k may tend to overfit. You may try to run the classification with different k within the preset range with cross-validation. Observe and find out the k that gives you the least cross validation sum-of-square error.

Another possible metric besides the sum-of-square error is called Silhouette applied to interpretation and validation of clusters of data. The technique provides a succinct graphical representation of how well each object lies within its cluster.


Decision directed TED (timing error detector). What does it mean? I am reading about symbol rate recovery and there is a term that I didn't grasp - 'decision directed' TED. Can somebody explain what does it mean or provide me with accurate reading? Thanks

"Decision directed" means that the inputs to the timing error determination use decisions (best guesses) of what the ideal symbol is based on the sample read.
This is in contrast to "symbol-aided" where a known symbol is transmitted and the replica symbol can be used in the receiver. You can visualize how this would allow for much further error determination due to a much wider range before the results alias.
I have diagrams specific to phase detection, but well demonstrates the concept and utility of decision-directed approaches applicable to timing recovery as well. As an example of this for the purpose of phase detection (used in carrier recovery) below shows a Decision Directed Phase Detector for QPSK and showing at the ideal sampling location in the symbol the actual sample measured at $V_2$ from which we assume (best guess) that the symbol transmitted must have been $V_1$. This is just an assumption based on nearest neighbor and only serves to drive the correction loop to that state, even though the whole symbol may end up being rotated once that correction converges (and then other techniques such as pilot symbols or higher level inspection or the data itself can be used to determine if such rotation occurred).

In contrast below shows a QAM Constellation where a symbol-aided approach is used for phase detection. The lock range as limited by the nearest-neighbor feature of a decision directed approach is greatly extended as depicted in the figure where this shows an actual symbol received as the dark circle and based on the known location (the outer symbol at +45 degrees) a much greater error can be determined. In a Symbol-Aided approach the known symbols are transmitted over a portion of the symbol such as in a pilot which is then used for carrier and timing recovery, and can operate in much lower SNR conditions. Decision-directed approaches out-perform  other approaches at higher SNR levels where all samples are based on the received signal itself (such as the Gardner TED which typically does not use a pilot or decisions) since a decision when good removes all noise for that input to the timing detector.