How does modulation work in simple language? I am reading modulation from this site (am a beginner). In the end, they've given this example:

I am having trouble understanding how exactly they arrived at the two modulated results for AM and FM. I have attempted to understand it, and this is my attempt (I am assuming the $x$-axis is time axis):

Amplitude Modulation: At each time $t=t_0$, they have taken an increase factor as $k=\frac{\text{amplitude of modulating wave at } t_0}{\text{constant amplitude of carrier wave}}$. Hence, to obtain the modulated result at each instant, they've done $\text{Final displacement}=\text{Original displacement of carrier wave}\times(1+k)$. Negative amplitude of modulating wave implies decrease of carrier wave's displacement, and positive amplitude of modulating wave implies increase.
Frequency modulation: I'd say a similar increase factor is at work here, but instead, that factor is now affecting the frequency of the carrier wave instead of its amplitude.

I am sorry if these are beginner knowledge but I don't have any good book to study from. I checked Wikipedia for these two as well. My guess for amplitude modulation's working seems reasonably close to their mathematics, especially the $y(t)=[1+m(t)]\cdot c(t)$. But, I could not understand why an integral is employed in frequency modulation. What is it supposed to do?

I hope I have a reasonably detailed and specific question here. I am looking for simple language explanations. Thank you!

In order to understand the integral occurring in frequency modulation you need to understand the relation between phase and frequency. A sinusoid with frequency $f$ can be written as

$$x(t)=\sin(2\pi f t)\tag{1}$$

where the argument $2\pi f t$ is called the phase $\phi(t)$. The frequency is the derivative of the phase divided by $2\pi$:

$$f=\frac{1}{2\pi}\frac{d\phi(t)}{dt}\tag{2}$$

If the frequency is not constant, the relation $(2)$ still holds. Assume now that the frequency is a function of time $f(t)$. From ($2$), the phase can be determined from the frequency by integration:

$$\phi(t)=2\pi\int_{0}^tf(\tau)d\tau\tag{3}$$

In frequency modulation, the frequency of the carrier is changed according to the message signal $m(t)$ around a given (average) carrier frequency $f_c$:

$$f(t)=f_c+\Delta f\cdot m(t)\tag{4}$$

where $\Delta f$ is the maximum frequency deviation, assuming that $|m(t)|\le 1$. From $(3)$ and $(4)$ we can obtain the phase

$$\phi(t)=2\pi f_ct+2\pi\Delta f\int_0^tm(\tau)d\tau\tag{5}$$

From $(5)$, the frequency modulated signal can be written as

$$s(t)=\sin(\phi(t))=\sin\left(2\pi f_ct+2\pi\Delta f\int_0^tm(\tau)d\tau\right)\tag{6}$$

Your understanding of amplitude modulation seems to be correct. You simply multiply the carrier with a non-negative signal that is obtained from the message signal plus a constant. The constant must be chosen such that the sum remains positive.


Feature for exact matching of images I am trying to match images from a large dataset which exactly match the given input image.(images with deformations/transformations are treated as different)

I have tried euclidean distance but it isn't very scalable. I have also tried matching gradient magnitude and orientation pixel by pixel, but again it seems to work in approximate manner and is slow.

I believe, pixel by pixel matching wouldn't be apt for matching from such large dataset. Moreover, even if certain approximation is allowed, can the threshold be decided dynamically rather than prefixed ?

Make use of perceptual hashes. It is very fast to compute and is very lightweight, both in terms of memory and cpu consumption. They are represented by simple long integers and can be indexed using many types of data structures such as VP Trees:

http://www.phash.org/

If that doesn't work, you can extract SURF features, quantize them into visual words using a vocabulary tree or BOW and index these with a database structure such as inverted file index. Check out OpenCV. It has many of the tools I mentioned here.


Prove that a specific auto-correlation function is even Background (from digital signal processing):
The $\delta$ function is defined as below:
$$\delta[n]=\begin{cases}
1 \quad \mbox{if} \hspace{.4em} n=0\\
0 \quad \mbox{o.w.}
\end{cases}$$
In an LTI (Linear Time-Invariant) system, we denote the response to the $\delta$ function by $h$.

Question:

Assume that $h$ is defined as below:
$$h(n)=\sum_{k=1}^p a_k h[n-k]+G \delta[n]$$
where $a_k$'s are some defined coefficients, and $G$ is another fixed
coefficient like them. $p$ also a fixed value.
Also, let the auto-correlation function $\hat{R}(m)$ be defined as
follows:
$$\hat{R}(m)=\sum_{n=0}^{\infty}h[n]h[n+m]$$
Prove that $\hat{R}(m)=\hat{R}(-m)$

My try:
I managed to define a new variable $v=n+m$. Since $n$ was from $0$ to $\infty$, $v$ will be from m to $\infty$. Re-writing $\hat{R}(m)$ using the variable, we have:
\begin{align}
\hat{R}(m)&=\sum_{v=m}^{\infty}h[v-m]h[v]\\
&=\sum_{v=0}^{\infty}h[v-m]h[v]-\sum_{v=0}^{m-1}h[v-m]h[v]\\
&=\sum_{n=0}^{\infty} h[n-m]h[n]-\sum_{n=0}^{m-1}h[n-m]h[n]\\
&=\hat{R}(-m)-\sum_{n=0}^{m-1}h[n-m]h[n]
\end{align}
Now I need to somehow prove that $\sum_{n=0}^{m-1}h[n-m]h[n]$ is equal to $0$. However, I do not know how. It looks like a deadend.

Note: Denoting the functions like $f[n]$ instead of $f(n)$ in this context is due to the fact that we are working with digital signals where we assume that the domain consists only of the integers.

Let me get into more additional detail from what Peter wrote.
First, note that $h[n]$ is the output of an IIR filter with input $G\delta[n]$. Now, note that such system is causal (uses only past information), and that the input is the impulse function. We know that if the system is causal, the impulse response of such system is causal too. Thus, $h[n]$ is a causal signal: $h[n]=0,\forall n<0$.
Now, let's get into the autocorrelation function thing. From your procedure you have the following:
\begin{align}
\hat{R}(m)&=\sum_{n=m}^{\infty}h[n-m]h[n]\\
\end{align}
So, you want to write this summation starting from $n=0$ right? If $m=0$ we are done. There are other two case: $m<0$ and $0<m$. In your procedure you did:
\begin{align}
\sum_{n=m}^{\infty}h[n-m]h[n] = \sum_{n=0}^{\infty}h[n-m]h[n] - \sum_{n=0}^{m-1}h[n-m]h[n]\\
\end{align}
However, this basically assumed that $0< m$. Otherwise, the second summation does't make sense since the lower limit would be greater than the upper limit. What if $m<0$? Then you would have:
$$
\sum_{n=m}^{\infty}h[n-m]h[n] = \sum_{n=0}^{\infty}h[n-m]h[n] + \sum_{n=m}^{-1}h[n-m]h[n]
$$
Lets inspect both cases under the causality property of $h[n]$. Hence, for the first case (with $0< m$) we have:
$$
\sum_{n=0}^{m-1}h[n-m]h[n] = 0, \ \ \text{since } h[n-m]=0 \text{ for } n-m<0
$$
this is true, since from the summation index range we have $0\leq n \leq m-1$ or equivalently $n<m\implies n-m<0$ (thus $h[n-m]=0$).
Now, for the second case (with $m<0$) we have:
$$
\sum_{n=m}^{-1}h[n-m]h[n]=0 , \ \ \text{since } h[n]=0 \text{ for } n<0
$$
this is true, since from the summation index range we have $m\leq n \leq -1$ or equivalently $n<0$.
Thus, no matter what value of $m$ you have, the result is
$$
\hat{R}(m)=\sum_{n=0}^{\infty}h[n-m]h[n] = \hat{R}(-m)
$$
Note that this works due to the causality property of $h[n]$: this same reasoning works for any causal signal, regardless of its particular definition (as the one you gave with $h[n]$).


Can Morse Code be encoded using "silence" for dots and dashes? Morse code dots and dashes are often displayed as a sound, or light that is on for a short or long period of time.

What if the only way to communicate was by knocking.... (or similar encoding) that can only encode using extended (or short) durations of silence.

Is it acceptable to use morse code in this way (where short or long durations of silence are the dots and dashes)?

i used to be a ham radio person (like more than 40 years ago) and still know Morse code.  it's like learning to ride a bicycle; you don't forget.  between the elements ("dits" and "dahs") of a single character, the space of silence is supposed to be the same as the ON time of a single "dit".  call that a single unit of time.  so a string of dits would be at exactly twice the frequency as a stream of dahs.  but the space of silence between characters is supposed to be three units, the same length as a dah.

so, if you encode only the onsets of each element, whether it's a did or dah, you can't tell the difference between a lot of dahs in a row in the same character (like zero or "0") or the same number of dits in separate letters (like five adjacent "E" in sequence).  from just the onsets of the beeps, a single "0" (which is "dah-da-dah-dah-dah") sounds the same as "EEEEE" (which is "dit-dit-dit-dit-dit" at the same rate).  nor could you hear the difference from "EENE" ("dit-dit-dah-dit-dit" at the same relative onset times).

if you make the space between letters long enough that there can be no confusion between separated characters and the closer elements within a character, you can figure it out from just knocking.  but, given the element and character spacing set forth by the format of Morse Code, you cannot tell "0" from "EEEEE" from "EENE" with just knocks.


Calculate the Derivative of Gradient Field of an Image I meet a confusing thing in image processing recently....

Assume the image $x \in \mathbb{R}^n$, with its derivative (difference) matrix: $D^+ = \begin{bmatrix} D_h  \\ Dv \end{bmatrix} \in \mathbb{R}^{2n\times n}$ ($+$ means forward difference), also equal to $\nabla$. Therefore, it is natural to define the divergence: $\triangle = \nabla \cdot \nabla$.

I have seen some papers use $div = \triangle = D_h^-D_h^+ + D_v^-D_v^+ \in \mathbb{R}^{n\times n}$, where $-$ denotes the backward difference.

Here is my question: assume I want to calculate the $\frac{\partial\|\nabla x -p\|_2^2}{\partial x}$ where $p\in \mathbb{R}^{2n\times 1}$ is a vector not related to $x$, what is the result? I have seen some authors use $\nabla\cdot (\nabla x -p)$.

However, if writing the $\nabla$ as matrix form $D$ as I have introduced before, $D^T$ is exactly adjoint of gradient, not backward difference. Hence $-\triangle x$ would appear! So what is the right formula? Could anyone tell me?

Consider the expansion of the term below
$$\|\nabla x -p\|_2^2 = (\nabla x -p)^T(\nabla x -p)$$
$$\|\nabla x -p\|_2^2 = (x^T\nabla^T -p^T)(\nabla x -p)$$
 $$\|\nabla x -p\|_2^2 = (x^T\nabla^T\nabla x - x^T\nabla^T p  - p^T\nabla x +p^Tp)$$

Now consider the following basic definition:

$$\frac{\partial(A x)}{\partial x} = A^T$$
$$\frac{\partial(x^TA)}{\partial x} = A$$

Now applying the above two definitions together with the expansion of the objective above to differentiate the objective we have

$$\frac{\partial\|\nabla x -p\|_2^2}{\partial x} = 2\nabla^T\nabla x - 2\nabla^Tp$$

$$\frac{\partial\|\nabla x -p\|_2^2}{\partial x} = 2\nabla^T(\nabla x - p)$$

since $\nabla^T = \nabla$, therefore we have

$$\frac{\partial\|\nabla x -p\|_2^2}{\partial x} = 2\nabla(\nabla x - p)$$ the constant 2 is just a constant, so the result we have is consistent with the ones that authors are using, its simply a consequence of vector differentiation


Why would we simulate baseband rather than passband My question is, Why do we usually use baseband in simulation rather than passband? however, if we want to test our simulated idea in real environment, it must be sometimes in passband !!

So, should simulation in baseband give the same results if we made it in passband? for term of BER, MSE and so on

Here is an example as asked in the first comments,

In baseband, Suppose you have ofdm system, with $s$ is the $n^{th}$ block with N symbol, when performing N-by-N FFT to get $x = FFT^Hs $, after adding a CP as guard interval with length P, we will have $M = x + P$ is the transmitted signal as baseband.

But in Passband, we will not transmit the $M$ as it's, we will have additional process when having $M$, this process is summarized  conversion --> matches filter --> carrier modulation by multiplying (.*exp(-j*2*pifct), fc is the carrier frequency t is the time period, and them transmit real part of this results, then in the receiver side, we follow the inverse of these steps till getting the baseband signal $M$ then, both passband and baseband will follow the same steps.

Thank you

One very obvious advantage of baseband simulation of communication systems is the low sampling rate it requires.

Considering a narrowband communication system having a few MHz of bandwidth to a few GHz of carrier frequency, it would be much easier (in terms of computation) to simulate it in the baseband sampling rate.

The reliability and the quality of a baseband simulation depends mostly on how accurate and realistic your mathematical modeling of the actual physical channel and the modulation system using it.


Non Coherent BPSK with no Channel State Information at the receiver Where or what is the proof for the failure of BPSK in fading channel when we have no knowledge of the fading coefficient $h$ at the receiver in the equation $y=hx+w$, where $x$ is BPSK symbol equally likely to be $A$ or $-A$, $h\sim \mathcal{CN}(0,1)$, and $w\sim \mathcal{CN}(0,1)$.

In any PSK modulation, the information is only in the phase.
An unknown $h$ randomly rotates your phase. You know nothing about the phase that $x$ had when you observe $y$. Hence, you have zero mutual information.
One could write that down as a very short formal proof based on conditional entropies, but I don't think that'll help you; the more intuitive understanding "without knowing $h$, the output phase is fully random" is sufficient as proof.


PyWavelets CWT: resampling vs recomputing wavelet Related. The implementation pre-integrates a wavelet once, and resamples it at each scale, finally differencing to implement below:
$$
C_{a, b} = \frac{1}{\sqrt{a}} \sum_k s(k)\left( \int_{-\infty}^{k+1} \overline{\psi \left(\frac{t-b}{a} \right) }dt - \int_{-\infty}^{k} \overline{\psi \left(\frac{t-b}{a} \right) }dt \right)
$$
But it also does * sqrt(scale), whereas above we're clearly dividing; what's the deal? And how's it compare to actually recomputing the wavelet at each scale?

The normalization is indeed by 1 / sqrt(scale), and it's an L2-norm; the trick's in the scale wavelet.
I'll use wavelet='morl' throughout. Prior to integrating, we can inspect the wavelet here; it's returned by wavelet.wavefun, which is binary-compiled, but after some guesswork, I found to it to match exactly with
scipy.signal.morlet2(1024, 1024 / 16) * sqrt(1024 / 16) * np.pi**(.25)

which is, from source, using Wiki's notation, $\psi(t) = \psi_{\sigma}(t/a)$, where $a$ is scale, and
$$
\psi_{\sigma}(t) = e^{j\sigma t} e^{-t^2/2} \tag{1} \label{1}
$$
(scale and $\pi^{-1/4}$ cancel). This is what's integrated via cumsum(psi) * step, and then resampled for all scales.

Resampled vs. recomputed
What exactly is this resampling doing in in terms of Eq 1? Is it just a higher resolution of the wavelet at the same scale, or is it equivalent to recomputing Eq 1 at each scale? Conveniently enough, latter, but only approximately, and approximation degrades substantially for small scale (-- code1):

Notice from code1, however, the recomputed wavelet:
Ns = len(int_psi_scale)
w = real(morlet2(Ns, Ns / 16) * sqrt(Ns / 16) * np.pi**(.25))  # repeats first blob
w /= scale

The recomputing includes 1 / scale, which, together with * sqrt(scale), equals 1 / sqrt(scale). Mystery solved.

Your code is wrong, where's * step?
Replaced by 1 / scale. How?
In MAE code, notice that for scale=64, we have int_psi_scale == int_psi, which == cumsum(psi) * step. For w_int, we do cumsum(w) / scale. And 1 / scale is... == step. Thus, the pre-integrated wavelet, psi, is just w at scale=64 (in morlet2 code above, 1024 / 16 == 64, checks out), and step happens to be ... conveniently? == 1 / scale when integrating.
Then, why is 1 / scale there? Unclear. Two possibilities in mind: (1) preserving the wavelet's norm at integration; (2) scaling the wavelet, independent of integration.

If the wavelet was L1 or L2 normalized before integration, either will be preserved. This is from chain rule; just replace $f$ with $\psi$, and $k$ with $1/a$:

$$
\int f(k x) dx = \frac{1}{|k|} \int f(x) dx
$$

This seems likelier, as the diff later is closely equivalent to undoing the integration, defeating the purpose of (1). Why rescale the wavelet? Normalization - see here.

You cheated earlier; there is no w /= scale
True, the code actually shows w_int = cumsum(w) / scale, but the two are exactly the same. It's thus the earlier question of where 1 / scale "belongs", or "comes from". This is answered here, and in other part below.

Why step == 1 / scale at integration? (-- for reference, from here (in code1, $n$ is x)):

Just a coincidence, or is step, along $n_i$, carefully crafted to yield the convenient resampling properties, which in turn demand step = 1 / scale? Don't know, might update answer later.


To find the filter Co-efficients of the following equation Given $y[n]=h_1[n]x[n]+h_2[n]x[n-1]+h_3[n]x[n-2]$ is an LTI system with unity gain at $\omega =0$ and zero gain at $\omega =\pi$. $h[n]\neq0$ .Also given that the system has a linear phase. Compute

$h_1[n]$

$h_2[n]$

$h_3[n]$

How do we approach a question like this? Any help is appreciated. I want to confirm my working for $h_1[n]$ and I found other two co-eficients as 1/4

You're overcomplicating things here. There's no need for sines and cosines and squares. Note that $\omega=0$ corresponds to $z=1$, and $\omega=\pi$ corresponds to $z=-1$. From the definition of the $\mathcal{Z}$-transform you should be able to figure out that the DC term of the transfer function $H(z)$ is given by

$$H(1)=h_1+h_2+h_3\tag{1}$$

and the value at Nyquist ($\omega=\pi$) is

$$H(-1)=h_1-h_2+h_3\tag{2}$$

Equating $(1)$ and $(2)$ with the desired values at these two frequencies gives you two equations for three unknowns. The third equation comes from the linear phase requirement. I'm sure you can take it from here.


Should the transmitted signal be continuous in FMCW? I am a novice in this field and working on respiration detection through near ultra-sound between $18 \textrm{kHz}$ and $20 \textrm{kHz}$.

In the referece article and other radar use articles, the generated signal was transmitted continuoulsy such as the picture below:

My question is what if the interval would be allowed between the signals such as the picture below:

I would like to know why the signla be tranmitted continuously and whether it is possible to be transmitted with intervals.

Thank you in advance.

Let us say the sweep repetition interval is $T_s$, it is the part with frequency modulation. If interval $T_i = T_s/N$ for large N, then drawbacks are: loss of average power, decrease of maximum measurable unambiguous Doppler frequency and it might cause amplifiers to introduce unwanted transients in the signal. Therefore, there is no advantage of adding an interval.

However, if interval $T_i = N*T_s$ for large integer value of N, then it is no longer called an FMCW, rather it will be called LFM, linear FM modulated pulse signal. The processing of which is more complicated than that of FMCW.


Question about fading Could you please help me to answer this question ?

We have a Shannon formula for Capacity as
\begin{align}
   C=B\log_2(1+\frac{h*P}{B*N_0})
\end{align}
where h is a random variable which denotes the channel gain with channel mean power
\begin{align}
   \Omega=E[h]
\end{align}

When I do the simulation such as for Rayleigh fading. I have to set a value for $\Omega$. But I don't know if $\Omega > 1$ or it must be smaller than 1. I have read many papers from IEEE. Some of papers set value of $\Omega > 1$, some of them is smaller than 1.

Could you please explain for me?

1- When it is greater than 1 and when it is smaller than 1.

2- In practice, how the range of $\Omega$ is ?

Consider two random variables $x \sim N(0, \sigma^2)$  and $ y \sim N(0, \sigma^2)$.
Then $h = \sqrt{x^2 + y^2}$ is Raylegh-distributed with parameter $\sigma$, whose p.d.f is

$$
f(h) = \frac{h}{\sigma^2} \exp(-h^2/2\sigma^2)
$$

According to Wikipedia,

$$
\Omega = E[h] = \sigma \sqrt{\frac{\pi}{2}}  \implies  \Omega \propto \sigma
$$

Hence, the range of $\Omega$ depends on the variance of your input signals $x$ and $y$. If you think of $x$ and $y$ as "signal + AWGN", then $\sigma$ is essentially the signal-to-noise ratio. So, whether $\Omega < 1$ or $\Omega > 1$ depends on whether you are modelling low SNR conditions or high SNR conditions, respectively.


Proper way to detect and decode transmission in gnuradio I am receiving AM modulated data over USRP. I successfully receive the data using an IF frequency (i.e., offset from the carrier), apply a highpass to remove RXLO phase noise and then downmodulate the signal to DC to get I/Q samples (I apply additional filtering and processing).
I am interested in the received data stream and the phase (i.e., the argument of the complex I/Q samples).

Now this is all straight forward if the data is coming in continuously. It works and I see the bits, buried in noise, rushing over my screen.

However, I am receiving the data asynchronously. So I first need to detect what is actually a transmission. Then decode the data and then process them.

What is the proper way to do this?

A simple power measurement of the signal does not seem robust to me since power- and noise levels are not static.

It sounds like you are talking about a "bursty" transmission. This is common in real-world communications systems (for example, WiFi and Bluetooth operate this way).

What you say is exactly right - relying on power levels alone is not a reliable way of detecting a transmission. A real-world bursty communications system solves this problem by inserting some amount of completely "known" signal into every transmission. This carries no payload information whatsoever - it is specified by the communications standard (so both the transmitter and the receiver know it) purely to aid reception.

In many cases, this known redundant signal is inserted somewhere at the start of the transmission, in which case you might see it referred to as a "preamble". (However, in unusual cases, like TETRA which is what the police radios use in much of Europe, a "midamble" is used instead). In other contexts, this kind of known signal might be referred to as a "pilot" or "training sequence".

Although it might seem wasteful to transmit this redundant signal, it turns out to be very useful - not only for detecting a transmission, but also for performing other tasks that improve the reception (such as timing synchronization, frequency correction, gain control and channel estimation).

Therefore, if you have control over the data that is being modulated, then you should try putting some known data at the front of every transmission. Then, at the receiver, (in the simplest case) you can correlate using this known sequence in order to reliably detect the transmission. There is quite a lot of study behind how to choose a sequence that will produce really reliable detections in different scenarios. You could read about M-sequences, Gold sequences, Hadamard codes, Barker codes and many more.

If you don't have control over the transmissions, then you need to find some information about them. If you really know nothing about their structure, and can't estimate anything, then unfortunately power level detection may be the best you can do.


Calculating THD where peaks are not exact harmonics I have captured an FFT of a 1 kHz signal that I am generating using a WM8904 codec running at 48k samples/second connected via I2S to a ARM microcontroller (Atmel SAME70).  It appears to have some jitter in it, which I believe is due to a small delay in restarting the DMA after each block (which I know how to fix, but I was more interested at this point in the particulars of taking this measurement).  I am using a Agilent MSO7054B to make the measurements.

The fundatmental is clearly visible at 1.0 kHz.  But the next peak is at 2.1 kHz, not 2.0; the one after that at 3.2 kHz, etc.  Is that due to the jitter?

Can I still use those peaks to calculate a THD value?  When doing so, I got a value of 4.12%.  Does that seem reasonable?

Jitter should pretty much be zero-mean, otherwise it'd be a frequency offset.

I've zoomed into your screenshot and it seems the horizontal distance between two points in the FFT plot is about 5 pixels, and your picture happens to be 1000 pixel wide – wild guess, this is an 256-point FFT.

In that case, the frequency bin spacing is $\frac{f_{sample}}{256}=\frac{111\,\text{kHz}}{256}\approx433.6\,\text{Hz}$; that explains why you see an "unsharp" peak for your fundamental, and why your harmonics don't happen to land in exact multiples of where you presume they should be. (The wide peak is just leakage in action – oscillations that don't fit in one DFT window of observation leak energy in multiple DFT bins.)

Either, use a signal frequency that is a multiple of the resolution calculated above, or change the sampling rate or DFT length. Simple as that: make your oscillation fit an integer time in the 256 sample points that your DFT is long.

Generally, you'll get finer frequency bin spacing and a more accurate PSD estimate with larger DFT sizes – so, if your scope allows you to somehow export a couple thousand points to your PC, just do the frequency analysis there; it's really not that hard using python these days :)


Recovering phase from spectrogram In this answer, an iterative algorithm is given for recovering a signal from its spectrogram, assuming overlapping windows:

$$
x_{n+1} = \operatorname{istft}(S \cdot \exp(i \cdot \operatorname{angle}(\operatorname{stft}(x_n)))),
$$

where $S$ is the spectrogram and $\operatorname{(i)stft}$ is the (inverse) short-time Fourier transform. I tried this out using NumPy, SciPy, and librosa:

import librosa
import numpy as np
import scipy

def mse(phase_pred, phase_true):
    '''
    Calculate the mean square error between the true phase and the
    predicted (reconstructed) phase.
    '''
    return np.mean(np.angle(phase_pred/phase_true)**2)

# Load an audio file and calculate STFT.
x, sample_rate = librosa.load('audio.wav', sr=44100)
D = librosa.stft(x)
mag, actual_phase = librosa.magphase(D)

# Try to reconstruct the phase using the iterative algorithm above.
phase = np.exp(1.j * np.random.uniform(0., 2*np.pi, size=actual_phase.shape))
x_ = librosa.istft(mag * phase)
print('iter {} mse {}'.format(-1, mse(phase, actual_phase)))
for i in range(100+1):
    _, phase = librosa.magphase(librosa.stft(x_))
    x_ = librosa.istft(mag * phase)
    print('iter {} mse {}'.format(i, mse(phase, actual_phase)))
    if i % 10 == 0:
        scipy.io.wavfile.write('recons{:05d}.wav'.format(i), 44100, x_)

As far as I can tell, it isn't converging:

However, the audio after the 100th iteration definitely sounds closer to the original than the one with randomized phases.

Why isn't this algorithm converging to the correct phase? Am I misunderstanding its purpose?

It should be clear from the definition of the STFT that the result is invariant under global phase changes. With the restriction to real input signals that still provides for a sign ambiguity, or a phase difference of $\pi$. That means, even if your reconstruction is perfect, the result may differ by that phase from the original.

Another reason for the failure to converge may be the existence of nearly disconnected components in your signal. Say you have two sound components that are temporally separated so that there is no STFT frame that overlaps with both. Each component may be reconstructed with a different sign, giving you four different reconstructions with identical STFT. The same works for components that are separated well enough in frequency direction. Each component doubles the ambiguity of the reconstructed signal and will increase the probability of a deviation from the original signal.

The final reason I can see lies in the nature of the iterative algorithm. Fixed point iterations, of which kind this algorithm is, can fail to converge to a proper fixed point, keep oscillating or even diverge. Unless the maps involved are well behaved in a certain sense, predicting convergence is hard. It seems that your iteration doesn't settle, and it's unclear if adding iterations will help. So this may apply here, together with the other two problems.


How can I compute a log-spaced power spectrum? I would like to compute a power spectrum in which the frequencies are logarithmically spaced.

In Welch's method there is a trade-off between the frequency resolution of the resulting power spectrum and the number of averages (i.e. error in the result).  I would like this trade-off to be dynamic, i.e. do fewer averages for low-frequency points in order to have a finer resolution at low frequency.

Is there a standard way to do this?

I suppose one way would be to initially do pwelch with a very high resolution (low number of averages), and then rebin the resulting spectrum using logarithmic binning.

I found a paper that addresses this question directly:

Michael Tröbs and Gerhard Heinzel, "Improved spectrum estimation
 from digitized time series on a logarithmic frequency axis," in
 Measurement, vol 39 (2006), pp 120-129. (Or free pre-print.)

The first few figures in the paper nicely illustrate the problem that this algorithm solves, and the references contain a useful bibliography of other approaches (constant-Q transform, tempered Fourier transform, a survey article, etc).

Their approach is not to re-bin the output of an existing FFT-based power spectrum estimation, but to only compute the discrete Fourier transform at the (logarithmically-spaced) frequencies of interest.  For each frequency to be estimated, they basically implement Welch's algorithm, but with a transform length (and hence also, number of averages) specifically chosen for each frequency. The computation of each frequency bin uses the entire time series, but segmented differently. The results have the desirable property that the resolution (bin width) is a smooth function of frequency, and the results can be calibrated as either a power spectral density or a power spectrum.

Matlab implementation here: https://github.com/tobin/lpsd

Disclosure: The authors of this paper are at the same institution as me.


unit step signal as sum of impulses from MIT OCW video From MIT's OCW on Signals and Systems by A. Oppenheim, at 4:30 mark https://www.youtube.com/watch?v=jGk3w1b7UXQ&index=3&list=PLLNp7XoiSLQYygYw8Mzt763zZCQdzCZcd

He said that by forming the sum from minus infinity to some value n of a unit impulse function is equal to the unit step function.
$$u[n]=\sum_{m=-\infty}^n \delta[m]$$

At n<0 the sum is accumulating nothing. We see that indeed it is equal to the unit step since the unit step is zero at n<0.

However, I am confused why at n>0, the sum can be equal to the unit step. From the diagram in the video, we see that there is only one impulse value that can be accumulated in the sum. But then the unit step is not composed of just one impulse, rather it is composed of an impulse train.

The expression is not adding an impulse train since the impulse function is defined only at m=0 which is 1. At m=1,2,3..the impulse function is already zero, so there is just one impulse that is contained in the expression.

What is my wrong assumption/understanding?

yes there is only one term in the sum when $n \ge 0$
$$ \delta[m]=\begin{cases}1 \quad m=0 \\
0 \quad \text{otherwise} \end{cases}$$

$$ \begin{align}\sum_{m=-\infty}^n \delta[m] &= \delta[-\infty]+ \dots +\delta[0] + \dots + \delta[n] \\
& = 0 + \dots + 1 + \dots + 0  \\
&= 1  \; \text{for} \; n\ge 0
\end{align}$$
The sum holds the value of 1 for $n\ge0$. His figure shows the terms inside the sum, not the sum of terms.

The definition of $u[n]$ is
$$
u[n]=\begin{cases} 1  \quad \text{for}\; n \ge 0 \\
0  \quad \text{for}\; n<0 \end{cases} $$
They are equivalent.  Defining $u[n]$ function as a sum of delta functions may seem strange but I'm sure he has a "trick" to show that will be useful down the road.

There is a similar relationship in continuous time relating the unit step function $u(t)$ as the integral of the Dirac $\delta(t)$.


How to make the error of the Gardner timing error detector not depend on signal amplitude? I have implemented a Gardner timing error detector with an interpolator and timing recovery in software. I noticed that the error term of the Gardner TED depends on the signal amplitude, which causes problems with the large variations of signal strength my receiver is expected to see.
With bursts (no timing sequence, only CW + sync word over 64 symbols), there are variations of 10 dB from the start of the transmission to after a few (15-20) symbols of frame data. For continuous signals (long CW + long timing sequence), this is less of an issue as the signal level tends to be relatively constant or varying rather slowly.
Is there a way I can make the Gardner error factor not depend on signal amplitude or should I instead focus on making sure the TED always sees a constant power level?

The maximum absolute error term is roughly 2 times the input power for the Gardner TED. Having a fixed (or predictable) input power, through the use of an AGC, means you can scale the error term to adjust the speed at which corrections will be performed.
With M-QAM modulations (for M > 4), the input power will vary depending on the symbol, and thus the maximum error term will also vary. This means that symbols with less energy will move the timing corrections more slowly than symbols with more energy. Given a proper initial synchronization sequence, this should not be an issue when maintaining the lock. Proper care should be taken to make sure the lowest energy symbol can still make the timing offset move if necessary.


Forward error correction with limited possible errors I am using a Reed-Solomon FEC implementation in GF(1024) to correct T=15 errors in a given codeword.

In the particular implementation though, I know for a fact that for each transmitted symbols, the corresponding received (potentially error'ed) symbol will always be one of a given set of symbols. The set depends on the value of the original symbol, and can contain anywhere from 32 to 243 members.

So for example if my original symbol is A0='0000000000', then the received symbol will always belong to a known set S(A0) which obviously contains A0 and some other symbols, never to exceed 243 symbols total per set. If A!=B then S(A)!=S(B) (but the intersection of S(A) and S(B) in general is non-empty).

Can you think of some optimization that I can use in order to reduce the parity of the code, given that I know that there is a limited number of possible received symbols for each transmitted symbol?

Thanks,
-Dimitri

Reed-Solomon codes are perfectly matched to a symmetric channel for which the probability that a given symbol $X$ is changed to a symbol $Y \neq X$ is the same for all choices of $Y \neq X$. That is, no particular symbol error is more likely than any other symbol error.

Thus, knowing that the channel is not symmetric and that certain symbol errors are impossible doesn't help in the decoding process at all. This knowledge might be usable in the post-decoder processing. If the decoder has corrected $Y$ into $X$ and $Y \notin S(X)$, then that is a dead giveaway that something is awry, right? So, the post-decoding processor can reject such purported "correct decoding" and ask for a retransmission or blank out that codeword to produce a blank area on the screen instead of gobbledygook.

Unfortunately, the most likely problem with Reed-Solomon decoders is not that they find an incorrect codeword, but rather that they fail to decode because the received word is not at Hamming distance $t$ or less from any codeword.

I don't think you are going to find any way of increasing the rate of the code based on your knowledge of the channel's asymmetry.


High / Low Pass Filters from Zero I everyone,

I'm a computer scientist but not from the signal processing community. I want to implement a real time digital low pass filter (to begin with) for a sound processing application. I've done some research and found some solutions:

Wikipedia gave me this filter: y[i] := y[i-1] + α * (x[i] - y[i-1]),
which is exactly the kind of formula I want.
I've found an online filter generator, that can generate C code from the description of a Butterworth filter. Great.

So my question is not really "give me a filter algorithm". Its more:
How do you make those formula?

What I want is not a ready-to-use algorithm that I can copypast. I want to understand what I'm doing and how it works. However! All the documentation that I've found until now, including the Wikipedia pages, are directed to people that already have a solid background in signal processing (a lot of unexplained notations, function names, etc.). This is not my case.

I need someone to explain to me how a low pass filter works, how the Butterworth filter (or something else I don't care) works, how I can extract a nice algorithm from it with the correct coefficients. And all that like I'm a 5 years old.

I haven't found anything like that yet. If you'r willing do it here you'r welcome. Or you can give me links to online literature, blog posts, anything that popularize all thoses notions. You'll make me happy.

Thanks <3

EDIT: I want to clarify that I'm not looking to know how a general low pass filter works, using a convolution and all that. That I already know. I mentioned the Butterworth filter and associated algorithm because it is a linear filter (algorithmically speaking). Something that I can reasonably use in a real-life application for real time processing.

I will try to be intuitive more than rigorous.
Low Pass Filter means a filter which limits the ability of the output data to change fast (Relative to a parameter called the Low Pass Cut Off Frequency).
Namely it make the signal transient to take longer (Smoother).

How do we do that?
We, make adjacent samples have similar values.

So for applying a Moving Mean Filter (Output is the mean value of few adjacent samples) means we applied a Low Pass Filter.

Another way is Auto Regressive model, which you pointed.
In that mode the output sample is the previous samples with a small addition which again, with proper ratio between the weight of the previous output to the new input means the speed of changes is limited.

This is the intuition.
If your needs are basic this might work (Just apply some weighing with non negative values, Let's say Gaussian Filter).
If you need more advance control you'll need to learn the frequency Representation of the sampled data.


Getting non smooth spectrum after performing fft on numerically obtained acceleration data I am trying to calculate the spectrum of Bremmstrahlung, which involves calculating the Fourier transformed acceleration. I am solving a non-linear ODE to numerically calculate the acceleration in the time domain. After taking the Fourier transformation using Numpy's fft, the resultant spectrum looks highly non-smooth and "non-physical" .
I cannot paste the entire code so I am posting what I think is relevant snippet. Can someone point out what I am doing wrong?

Note: My acceleration is a function of two variables (beta, and b, the impact parameter), and I want to plot it for the different b, and in order to factor out the beta term I am just summing over all the values of acceleration for the different beta, for a given impact parameter b.

Also my spectrum is Fourier transformed acceleration square times a constant factor (Larmor's formula)

'''
Fourier Transformation of the acceleration
'''
N = 2**8
sampling_frequency = 10
a_w_normalized = [[] for a in range(len(impact_parameter))]
a_w =[[] for a in range(len(impact_parameter))]
intensity_normalized = [[] for a in range(len(impact_parameter))]
acc_summed_over_velocity = []
intensity_summed_over_velocity =[]

for index,b in enumerate(impact_parameter):
    acc_sum =[0 for x in range(N)]
    intensity_sum = [0 for x in range(N)]
    for j in range(len(velocity_z_component)):
        #window_kaiser = signal.kaiser(N, 15)
        #window_hann = signal.hann(N,sym=True)
        window =1
        fft_input = acc_normalized[index][j]*window
        ft_acc_normalized = np.abs(np.fft.fft(fft_input,norm=None))

        acc_sum =np.add(ft_acc_normalized,acc_sum )
        intensity_list = [power_spectrum_factor * (a ** 2) for a in ft_acc_normalized]
        intensity_sum = np.add(intensity_list,intensity_sum)
        intensity_normalized[index].append(intensity_list)
        a_w_normalized[index].append(ft_acc_normalized)
    acc_summed_over_velocity.append(acc_sum)
    intensity_summed_over_velocity.append(intensity_sum *acceleration_factor**2)

#intensity_summed_over_velocity=  ma.masked_less_equal(intensity_summed_over_velocity,1e-5)

'''
Plotting acceleration for selected values of impact paramters in time and frequency domain
'''

if plot is True:
    plt.figure(figsize=(12, 8))
    for i in range(len(impact_parameter)):
        plt.plot(t,acc_timedomain_summed_over_velocity[i], label='b={:.3f}'.format(impact_parameter[i]), )
    plt.legend()
    plt.ylabel(r'$ a(\tilde t) $', fontsize=14)
    plt.xlabel(r'$ \tilde t $', fontsize=14)
    ax.spines['left'].set_position('center')
    ax.spines['bottom'].set_position('center')
    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)
    ax.spines['left'].set_smart_bounds(True)
    ax.spines['bottom'].set_smart_bounds(True)
    ax.xaxis.set_ticks_position('bottom')
    ax.yaxis.set_ticks_position('left')
    plt.xticks(fontsize=15)
    plt.yticks(fontsize=15)
    if screening is False:
        plt.title('a(t) vs time without screening',fontsize=15)
        plt.savefig('./Plots/acc_without_screening.png')
    else:
        plt.title('a(t) vs time with screening',fontsize=15)
        plt.savefig('./Plots/acc_with_screening.png')

    '''
    Fourier transformed intensity for different impact paramater
    '''

    plt.figure(figsize=(12, 8))
    freq_normalized = np.fft.fftfreq(N)*(2*np.pi*sampling_frequency)
    for index,b in enumerate(plasma.impact_parameter):
        plt.plot(np.abs(freq_normalized), intensity_summed_over_velocity[index], label=r'$ \tilde b={:.2f}$'.format(b), )
    plt.xlabel(r'$ \tilde \omega $', fontsize=14)
    plt.ylabel(r'$ I_{\omega} $', fontsize=16)
    plt.legend(loc="lower left")
    plt.xscale('log')
    plt.yscale('log')
    plt.xticks(fontsize=12)
    plt.yticks(fontsize=12)

    if screening is False:
        plt.title('Single particle spectrum without screening')
        plt.savefig('./Plots/spectrum_no_screening.png')
    else:
        plt.title('Single particle spectrum with screening')
        plt.savefig('./Plots/spectrum_debye_screening.png')

    plt.show()

Most of the time, a non-parametric PSD estimation based on FFT alone (called a Periodogram) will provide a random looking spectrum of the numerical data.

To reduce those random variations, or to get a smoother looking estimate, you can do one of the two things:

Use Periodogram averaging such as Welch's method.
Use a model based parametric PSD estimate instead.

Note that the first approach will trade-off spectral resolution for reduced variance, whereas the second method will be high resolution but requires a consistent knowledge of the associated physical model of the process.


Some problem with the expression of QAM basis and signal The expression of QAM basis is
\begin{align}
\psi_1&= \sqrt{\frac{2}{E_g}}g(t)\cos(2\pi ft)\\
\psi_2&= -\sqrt{\frac{2}{E_g}}g(t)\sin(2\pi ft)
\end{align}

I want to ask why is $\psi_2= -\sqrt{\frac{2}{E_g}}g(t)\sin(2\pi ft)$,
instead of $\sqrt{\frac{2}{E_g}}g(t)\sin(2\pi ft)$? Why is there a "-" (i.e. minus) before $\sqrt{\frac{2}{E_g}}$?
And QAM signal is $\ s_m(t)=A_{mi}\sqrt{\frac{E_g}{2}}\psi_1 +A_{mq}\sqrt{\frac{E_g}{2}}\psi_2$, does the $\sqrt{\frac{E_g}{2}}$ here have a relation with its basis $\sqrt{\frac{2}{E_g}}$ ?

Base choices are ambiguous. The $-$ is there because the author of the book you're reading decided it was better that way.

You'll probably find the cause for that in the quadrature mixer graphic that they use, where if the phase shift from $I$ to $Q$ is $+\frac\pi2$, since $\cos(x+\frac\pi2) = -\sin(x)$.


why don't we use Broadband technology to transmit digital signals ,and use Baseband technology to transmit analog signals? I learn the Broadband technology and baseband technology in this website: https://www.computernetworkingnotes.com/networking-tutorials/differences-between-baseband-and-broadband-explained.html
Broadband technology uses analog signals in data transmission
Baseband technology uses digital signals in data transmission
I want to ask why don't we use Broadband technology to transmit digital signals ,and use Baseband technology  to transmit analog signals ??
I have two thinkings about this:
1.we cannot use Broadband technology to transmit digital signals?if not,why?
2.In fact we can,but we will cost  lots of price
these are two different things,but i don'y know which is right

"Broadband technology" is not inherently coupled with "analog signal transmission", as "baseband technology" is not coupled with "digital signal transmission". Both wire cables and optical fiber cables transmit signals in the form of electromagnetic waves, in RF frequency range for a wire medium and in IR/light frequency range for an optical fiber medium.
For digital transmission, the source data are transformed to changes in amplitude/frequency/phase parameters of a generated output signal and then the pulse shaping and filtering operations are performed on this signal. Resulting voltage/current or amplitude/phase waveforms have no sharp rise/fall jumps like those you see on the logic analyzer display. In this way, the resulting waveforms look similar to an analog signal.
In a baseband transmission the signal is sent "as is"; in a broadband transmission a modulator device creates a passband signal from a baseband signal and a carrier. The passband signal frequency is shifted by modulation process and occupies the frequency band around the carrier frequency with the width of the baseband signal spectrum. If the entire channel is used for a baseband transmission, we can use a medium's full bandwidth for the baseband signal. On the other hand, passband techniques allow frequency-division multiplexing, and we can share the channel for multiple concurrent transmissions. The number of subchannels is defined by a medium's full bandwidth divided by a baseband signal spectrum width.
Cable television systems use "broadband technology" for distribution of multiple analog video channels. Digital networks also can use frequency-division multiplexing (the "broadband technology", as your reference calls it). In 1980s, 10BROAD36 Ethernet standard was developed in an attempt to use analog TV cables for Ethernet networking. The 10BROAD36 communication line was even able to share a cable with TV channels. The attempt failed for a number of reasons. What's more, the frequency multiplexing is not a single (and so unavoidable) solution for sharing communication lines in packet-switching networks of Internet. The sharing of a communication line to multiple data streams is performed with packet switching. Still, frequency-division multiplexing is used in wired (Discrete Multitone Transmission) and wireless (WiFi, OFDM) digital transmissions.
As for analog transmissions in rare survived analog telecom systems, I do not know of analog signals with spectrum that wide as to occupy the entire bandwith of transmission lines. So these systems can safely use frequency multiplexing to share cables for multiple transmissions ("broadband technology", as your reference calls it).


how much of DSP (digital signal processing) is related VLSI? Although,this question may not be very much related to the topic of discussion, but i just wanted to know how much of DSP (digital signal processing) is used in VLSI. Like in the industry or for the sake of knowledge, which all topics share common between the two Domains?
PS:I know this may be off topic, but still thought that I could get answers from knowledegeable people on this platform.

First, an example: designing FIR filters pops to mind, or for that matter any signal processing module. The FIR filter is a nice example, since it can easily be implemented in VLSI even at the lowest abstraction level - physical mask level (in Cadence for example).

If we recall an example of the FIR filter design:

We can see that we only need to design:

Delay
Multiplier
Adder

elements, which are quite easy to implement in VLSI, either at the high level (VHDL, Verilog), or lower level (Cadence, Magic IC, etc.).

So now we have a mixed problem - the DSP part is designing the FIR filter, which can be designed in various ways, and dealing with the transfer function/frequency response. The VLSI part is actually implementing the filter - creating the layout which can be printed on a real integrated circuit.

In general, the two fields are connected, but in industry usually one deals in only a single layer of abstraction. Either one deals with circuits/elements as black boxes, without caring so much about the physical implementation, or conversely focuses on the optimal design of a component, which would later be produced and used in projects where DSP is employed.


SVD is not applicable to time series analysis? I am studying a deconvolution problem for which observations include noise. I am considering using one of several common regularization methods, including Tikhonov's solution, truncated Singular Value Decomposition (SVD), or Wiener filter.

A manuscript I am reading, Riel and Berkhout, "Resolution in seismic trace inversion by parameter estimation" in Geophysics, 50(9), pp. 1440-1455, 1985, compares analysis in the Fourier domain with general SVD analysis:

  Examples where SVD-based techniques may be used while Fourier
  techniques are no longer applicable are: time-invariant filtering,
  space-variant wave field extrapolation [...], and linear parameter
  estimation problems. Furthermore, in many nonlinear problems
  linearization is meaningful, again allowing application of the SVD.
  The nonlinear problem treated is such an example. In our opinion many
  problem areas in seismic exploration stand to benefit from an analysis
  involving the SVD. However, there are some limitations compared to
  Fourier analysis. For instance, SVD is not applicable to time-series
  analysis. Furthermore, the eigenvectors and eigenvalues do not, in
  general, have the clear physical meaning the sinusoids and their
  amplitudes in Fourier analysis do. For these reasons SVD-based
  methods should, when possible, not replace Fourier analysis methods,
  but rather be used in conjunction with them.

I was surprised to read that "SVD is not applicable to time-series analysis." It is not clear to me that this is the case, and indeed the authors use SVD to construct an inverse filter for their convoluted time series.

I wonder if there is common knowledge in the DSP community regarding SVD's unsuitability for time series analysis in general or specific aspects, as this is the first time I have come across this statement. (It is also surprising to find no SVD or eigenanalysis tag.) Thanks in advance.

SVD of multi-variate time-series is of course very common.  The point I think the author is making is that SVD cannot meaningfully be directly applied to a univariate time-series, instead it is commonly used indirectly to explore the main eigenvalues of a transition matrix, estimated from observed trajectories of a certain length. see http://en.wikipedia.org/wiki/Singular_spectrum_analysis .  SSA assumes the time-series is stationary and can be modeled as Markov.


Generate loopable sine wave cycle for given frequencies I am somewhat new to DSP so bare with me. I am working on a simple wavetable oscillator. What I am trying to do is generate one cycle of a wave form for each note (frequency) and store it in an array for later playback. I know I could generate one frequency, store it and resample to the other frequencies in real time but I'm looking for performance and I'm still learning so I'm trying to keep it simple for now.

In my first attempts I was able to generate one cycle of a sine wave at each frequency but when I went to play back the stored wave forms the frequencies were off by 0.5 - 0.10 hz or more and had a noticeable clicking noise when the cycle repeated. I know this is because the start and end of the wave form are not lining up on 0. I've been struggling to solve this problem.

Below is my most recent C# code I use to generate and store the sine wave. This will produce more cycles but seems to fix Hz variance I was getting with just generating one cycle. The only issue that remains is the clicking when it loops back. I may be going about this completely wrong so any suggestions would be greatly appreciated.

double frequency = 440.0; //A4
double sampleRate = 44100;
double numSamples = sampleRate / frequency;
double phase = frequency / sampleRate;

int j = 0;
for (double i = 0.0; i < numSamples; i += phase)
{
    double sample = Math.Sin(i * 2 * Math.PI);

    sineWaveTable[j] = sample;

    j++;
}

I'm not sure if it's even possible to generate a single cycle of a sine wave for each note and have it start and stop on 0. If not, am I stuck with having to use resampleing to get the other desired notes?

Your clicks are coming from two sources. The wavetable and the hardware.

If you are looking to create a simple little instrument that produces tones by playing back a wavetable, then there is no way to fix the size of the "playback window" to reproduce tones without clicks.

The frequencies that compose the chromatic scale follow a geometric progression that achieves doubling of the frequency in 12 steps. The frequency for each tone is produced by multiplying the previous tone by something like 1.05. That factor is worked out by: Let's take $A_F = 440Hz$ and its octave at $A_O = 880Hz$. Each one of the tones is $A_{n+1} = A_n * r$, so $A_O = A_F \cdot r^{12}$. But we know that $A_O = A_F \cdot 2$. Therefore, $r^{12} = 2$ and $r \approx 1.059463094359$

So, you might fit $A_F,A_O$ neatly in a wavetable, but $A_F \sharp$ now needs $\approx 1.0594$ more samples to fit in neatly. Not only there is no integer there but even if you could reproduce the tones in neat windows you would get problems with the tempo, the duration of each note in a musical piece, because $A_F$ is shorter than $A_F \sharp$.

The "classic" way around this is to setup a high resolution wavetable of a single sinusoid which is then read back at "different speeds", with interpolation. So, $A_F$ is read back at $r^0 = 1$ speed (and practically no interpolation) but $A_F \sharp$ is read back at $r^{1} \approx 1.0594$ speed and a bit of interpolation when the "needle" has to land between known samples.

BUT!, this has nothing to do with your output stream. In other words, you have to set-up (or adopt) a stream system where samples are continuously being pushed to the sound card and you have a "global" sense of time.

Once upon a time, you did this by setting up three pointers (that's in C, C++) to three buffers. Let's call them $H,D,L$. You set up the sound card to record in $H$ and playback $L$ and while that was going on, your main code was processing buffer $D$. When the sound card called its callback to signify that the recording ended, you would "circularly switch" the buffers around so that they now pointed to $L \rightarrow D, D \rightarrow H, H \rightarrow silence$ and re-set up the sound card to do the same job while you were now processing the $D$ which is now your freshly recorded $H$. Those buffers had a length, say for instance 1000 samples and the task then was becoming to find a sweet spot that balanced length of buffer with responsiveness.

But as you can see, in this system, samples are being pushed out to the sound card continuously in blocks of $bufferLength$ samples and you don't really care if what is being played back fits the buffer. All that you do is iterate your models (e.g. the wavetable), produced a buffer of samples, freeze the models, playback that buffer and move again from the top.

Now, obviously, the above includes real time input from audio sources too. If all you are doing is a tiny little synth, you could have two buffers and switch around between them, working on one while the other is being played back.

But even if you did this with some off the shelf cheap sound card in the market, you still got clicks although, you have worked so hard to perfectly align the zeroes between blocks. This is because of something called latency.

In other words, just because your program called some Operating System Application Programming Interface function that said "Record NOW!" (or "Play NOW!), it doesn't mean that this command will propagate through the driver and eventually the audio hardware and recording will start NOW. Recording (and playback) will commence some $NOW! + latency$ units of time later and in off the shelf sound cards, this used to be an awful time later. So awful that you could hear the silence between calls.

So, the solution there was to move over to something like ASIO which eventually packaged this whole workflow into a complete solution all the way from the software calls to the hardware.

As an example, an 8-track external sound card can do 48kHz Stereo playback / rec over buffers of 256 samples (with absolutely no clicks). This means, that sound is processed in blocks that last approximately 5ms. The human ear begins resolving sounds as distinct when they are apart for more than 40ms, so this 5ms is really safe. It means that we can plug an electric instrument in the sound card and process it, in real time like it went through an effects box and the performer will not be noticing any difference between what they play and what they hear.

If you tried to make $H,D,L$ the size of 256 and try to do the same with an off the shelf (relatively popular) sound card, first of all you would g/et/ /an/awf/ul/am//ount/of/cli/k/s// (one every 5ms) and secondly, the driver and operating system would be so much overwhelmed by the amount of function calls that soon you would be greeted with Seg Faults, blue screens, system freeze, reprepreprepreprepreprepreprepreprereprepreprepreprepreprepeated waves out of synsynsynsynsynsynsynsynsynsynreprepresync and other amusing side effects.

So, bottom line, keep working on the DSP side of your synth but make sure that your sound hardware is up to the task as well if you really want to get rid of the clicks.

Hope this helps.


Adjust Matlab Gui parameter while playing an audio I would like to design a very simple MATLAB Gui program to adjust the magnitude of an audio while it is being played via a slider (slider1). However, the program any adjustment will not response until the audio is finished playing. How could I fix this?

The audio is played after a push button is pressed

function PlayFile_Callback(hObject, eventdata, handles)
%Reading the audio file
FrameSize = 2048;
AR = dsp.AudioFileReader('Filename','test.wav','SamplesPerFrame',FrameSize);
Fs = AR.SampleRate;
%initializing audio player
AP = dsp.AudioPlayer('SampleRate',Fs);
tic
Tstop = 30;    %playing time

 while toc < Tstop
    gain = get(handles.slider1,'value');  %Get the gain factor
    audioIn = step(AR);
    audioOut = audioIn*gain;
    step(AP,audioOut);
end

As far as I know, the simplest, very old, matlab function to play wave PCM audio is "sound(x,Fs,Nbits)" and it is a high level function that will buffer the entire piece before it is played. And it will not allow adjustment on "x" before it is fully played. If your play function is also similar in mechanism, then you shall consider either using a low level audio playback function or divide your audio into small pieces and send each piece in sequence. However this last method may create problems due to multimedia timing issues. And as a result played audio may have pops & cracks in it.


How to apply Prewitt Edge Detection algorithm on a given image? I am trying to figure out the answer to the following question:

As I understand it, I am supposed to apply the algorithm only on this part of the image

[2,5,6]
[4,2,7]
[1,3,5]

My question is how to apply the two gradients for x and y respectively on this part of the image? Thank you in advance!

I was able to answer my question after doing more research. The key is to slide the 3x3 kernel and perform Gx and Gy (gradient for the x and y directions) and calculate the output image values. For example, to calculate the output image top left value -

[x - -]
[- - -]
[- - -]

we have to calculate both gradients (x and y)->

Gx

[-1 0 1]
[-1 0 1]
[-1 0 1]

we take the first chunk of the original image (e.g. sliding the kernel over the top values of the input image)

[1 2 4]
[5 2 5]
[5 4 2]

Gx = 4+5+2 - 1-5-5 = 0

Gy = 5+4+2 - 1-2-4 = 4
We take the absolute values of both and sum them so |4| + |0| = 4

And this is the first top value of the output image.

[4 - -]
[- - -]
[- - -]

The same way we calculate the rest. I hope that will help anyone else facing the same problem.


Why is GPS 1575 MHz? I’ve read that that frequency is a ‘window’ in weather, impacted less by raindrops for example. I otherwise would’ve thought that simply increasing the wavelength would decrease impact. I don’t understand why 1575 MHz is a ‘window’. Also, what factors other than weather impacted the decision to use 1575 MHz?

The technical decision was to use L-Band which extends from 1-2 GHz. The specific frequencies to use within L-Band depended on political decisions including existing allocations and users of those frequencies and possible interference. Unprotected GPS signals are very sensitive to interference rather than GPS interfering with other users.
L-band is attractive for GPS for these reasons:

L-band has good penetration of foliage, rain, snow and clouds and signal propagation is less influenced by weather.

Frequencies below 1 GHz and significantly above 2 GHz are much more affected by the ionosphere (resulting in variable delays).

The GPS waveform itself requires a relatively high bandwidth (C/A code has a main lobe that is approximately 2 MHz wide), where higher carrier frequencies are more likely to support wider bandwidth allocations.

The antenna pattern in L-band can easily realize a wide pattern for whole earth coverage and reception; higher frequencies tend to be much more directional in practical realization.

In addition to the civil L1 signal at 1575.42 MHz, two additional frequency carriers are emerging as part of GPS modernization for civil use in L-band, one at 1227.60 MHz (L1C) and another at 1176.45 MHz (L5 for aviation safety services). There will also be an improved civil signal at 1575.42 MHz (L1C) that will work in conjunction with the existing civil L1 signal and military signals there. A dominant feature of the new carrier for normal civil use will be the ability to do ionospheric correction due to the differences in the receiver of the signals at each carrier.


Time Variant or Time-Invariant system? $y(t)=\mathcal{T} \{ x(t) \} = \int_{-\infty}^{x(t)} x(\tau) d\tau$

Is the above system time invariant ?

The system is not time invariant as can be seen by the fundamental test procudere:

Given the system with an input/output relation:
$$y(t)=\mathcal{T} \{ x(t) \} = \int_{-\infty}^{x(t)} x(\tau) d\tau$$

First check the new output for a shifted input $x(t-d)$ as
$$y_d(t) = \mathcal{T} \{ x(t-d) \} = \int_{-\infty}^{x(t-d)} x(\tau-d) d\tau $$
$$y_d(t) = \mathcal{T} \{ x(t-d) \} = \int_{-\infty}^{x(t-d)-d} x(\tau) d\tau $$

And then also check if the output $y(t)$ shifted by $d$ yields th same formula?
We know that output $y(t)$ is given in terms of input $x(t)$ with the formula
$$y(t) = \int_{-\infty}^{x(t)} x(\tau) d\tau$$

Hence the shifted output $y(t-d)$ is obtained by simply replacing every $t$ variable with $t-d$ in that formula which yields:
$$y(t-d) = \int_{-\infty}^{x(t-d)} x(\tau) d\tau $$

And now looking at the formulas for both $y(t-d)$ and $\mathcal{T} \{ x(t-d) \}$ given above, we see that since $x(t-d)-d \neq x(t-d)$ for all $d$ then $y(t-d) \neq y_d(t)$ and we conclude that the system is not time invariant.


Why is level of power spectrum dependent on FFT resolution? I created a sinusoidal wave in some noise, and plotted the power spectrum of the signal using two periodogram estimates (welch procedure). One estimate is 'high resolution' - i.e. it uses a longer window than the corresponding 'low resolution' estimate. The result and code used to generate it are shown below. For reference: the y-axis will have units of |signal|^2 and the x-axis are just FFT samples.

My question is: From the plot, it would appear that there is overall more power content in the low resolution (red) plot, although this is obviously not the case, the input signal is the same in both cases. This poses some questions about how we compare signals, it seems it is critical to use the same parameter values (like window lengths) when comparing the energy content of different signals.

For reference, I do understand that the length of the window in the welch estimate will determine the resolution of my spectrum, and that the number of FFT points used will determine the smoothness of my output, which has nothing to do with resolution.

Fs = 1000;
t = 0:1/Fs:1-1/Fs;
x = cos(2*pi*100*t)+randn(size(t));

nWindow1 = 100;
noverlap1 = 0*nWindow1;
nfft1 = 1000;
[pxx1, f1] = pwelch(x,hann(nWindow1), noverlap1, nfft1, Fs, 'centered', 'power');

nWindow2 = 1000;
noverlap2 = 0*nWindow2;
nfft2 = 1000;
[pxx2, f2] = pwelch(x,hann(nWindow2), noverlap2, nfft2, Fs, 'centered', 'power');

figure
hold on
plot(f1, 10*log10(pxx1),'r')
plot(f2, 10*log10(pxx2),'b')
legend('Low Resolution Estimate', 'High Resolution Estimate')

The observed phenomenon is simply due to the fact that when you set window length to 100 it is as if you were doing fft with 100 samples instead of 1000. So Nw=100 similar to NFFT = 100.

And why does FFT size affect results (level of noise VS level of tone)? Good question

Two types of spectral densities:

1) PWELCH-like PSD: "PSDW"

A tone grows in dB as FFT size increases. As in a true PSD the tone is a Dirac's delta for NFFT->inf. Larger NFFT means smaller dF, and larger A, in such a way that A * dF = constant
A unit power noise stays at 0dB.

This is really a [Watt/Hz] spectral density and total power given by:

totalPower = trapz(F,PSDW)  =  trapz(fs * F,PSDW/fs)  = sum(PSDW)/NFFT

2) OSCILLOSCOPE-like PSD: "PSDO", PSDO = PSDW/NFFT

A tone remains the same dB as FFT size increases
A unit power noise becomes less dBs as NFFT increases.

This is a [Watt/Hz/NFFT] measure, with power:

totalPower = NFFT * trapz(F,PSDO) =  NFFT*trapz(fs * F,PSDO/fs)  = sum(PSDO)


Types of rounding in coefficients quantization Suppose we have create an IIR filter with matlab function "ellip",  and then we want to quantize the coefficients using:

\begin{align*}
bq=Quantize('round',b,2^8); \cr
aq=Quantize('round',a,2^8);
\end{align*}

I have read that there are 4 major types of rounding:

truncate
round
convergent rounding
round-to-zero

What is the differences between of them and how i know which method is best to choose?

The various rounding methods have a computation vs. quantization error tradeoff.

Truncate

Truncation is the simplest method.  Everything after the decimal point is simply lopped off.  For instance, both 2.1 and 2.9 become 2.  This is very simple, but is the worst method in terms of quantization error.  It is particularly bad because when you are dealing with non-negative numbers it introduces a strong negative bias.  In some algorithms that can be bad.

Round/Convergent Rounding

Simply saying "Round" doesn't tell you enough to know what is meant.  What kind of rounding?  People usually mean convergent rounding when they say "round", so I will assume that that is what is meant.

Convergent rounding rounds down when the decimal place is $\le .499\overline{9}$ and rounds up when the decimal place is $\ge .500\overline{0}1$.  The question remains- what to do when you have exactly $.5$?  Some rounding algorithms always round down (this is called "round-to-zero"), but that introduces a very small amount of bias.  Convergent rounding tries to eliminate the bias by rounding to the nearest even number, the assumption being that around half the time that will be up and half the time it will be down.

This is the best algorithm in terms of quantization noise, but is the most computationally intense.  In many situations though, you don't care how computationally intense it is since the calculation isn't done at run-time.

Round-to-Zero

As mentioned previously, round-to-zero is the same as convergent rounding except it always rounds towards zero when the decimal portion is $.5$.  This is a compromise between truncation (easy computationally) and convergent rounding (low quantization error).


What's the computational complexity of MMSE equalizer in OFDM system As I know that in case of one tap OFDM equalizer, (ZF equalizer), we should have $N$ division complexity since we divide every subcarriers by its correspondent channel tap in frequency domain.

But I couldn't calculate how much exactly the computational complexity in case if we used MMSE equalizer for the same OFDM system.

Could you please help me how to calculate that?

Thank you

For a SISO channel, assuming you have estimated channel coefficients for each sub-carrier $k$, the MMSE equalizer is
$$
\hat{h}_k=\frac{h_k^*}{|h_k|^2 + \frac{\sigma_x^2}{N_0}}
$$
So you can see already there is a multiplication in the numerator ($h_k^*$), and then there is a division by a term in the denominator. So for all $N$ subcarrier, this itself will be $2O(N) \rightarrow O(N)$ for large $N$. There is additional $|h_k|^2$ computation which is basically multiplication with conjugate. So if you have precomputed SNR ($\frac{\sigma_x^2}{N_0}$) and $h_k$ for all sub-carriers, complexity should still be linear on $N$. But this is a very simplified assumption because in most practical OFDM packets, there a pilot sequence periodically present so you may need to update $h_k$ using available pilot sequence.


confused about time shifting property of Laplace Transform In book signals and systems 2 edition a question is given which is as follows:

$$      x(t)=e^{-3(t+1)}u(t+1)      $$

and we are asked to find the unilateral Laplace Transform of the signal. The method that is given in the solution manual is as follows:

Using Table 9.2 and time shifting property we get:

$$     X_2(s)  = \frac{e^s}{s+3}      $$

Now I am given a question which is as follows:

$$      e^{-2t}u(t-1)       $$
and asked to find the Laplace Transform. Now can I apply the method as used above for unilateral Laplace Transform and get:

$$             \frac{e^{-s}}{s+2}       \rightarrow A   $$

Or does that method only holds true for unilateral Laplace Transforms? Because the answer marked A is wrong when I use this method. Also tell me when can I apply the property?

If you have written the function correctly then its Laplace transform could be found very similary to your first example:

Given
$$x(t) = e^{-2 t} u(t-1)$$
its Laplace transform could be found as follows. First denote the signal

$$x_0(t) = e^{-2} e^{-2t} u(t) $$

then its obvious that
$$x(t) =x_0(t-1) $$

Using the tables and properties to conclude:

$$X(s) = e^{-s} X_0(s) $$
$$X(s) = e^{-s} \frac{ e^{-2} }{s + 2} = \frac{ e^{-(s+2)} }{s + 2}   $$


Gaussian noise generator algorithm trouble I stumbled upon the following algorithm for Gaussian noise generation and I just can't figure out how this is supposed to work at all.
The algorithm is as follows:

First of all: $j = \arg\min_i (q_1 - p[i])$ does not depend on $q_1$ and therefore yields the same result every time - I assume $\arg\min_i (|q_1 - p[i]|)$ is what's really meant here.
But then there's more: I can safely assume that $q_1$ is distributed linearly within $[0,1]$, otherwise I could just scale that, round afterwards and wouldn't have so much hassle.
In this case I got another problem: As $q_1$ is linearly distributed, choosing $j$ near to zero (which is what $\mu=0$ would implicate) is not particularly likely, so I suspect the resulting distribution is not a normal distribution.

This algorithm reminds me of inverse transform sampling ( http://en.wikipedia.org/wiki/Inverse_transform_sampling); however it does not appear to be a variant of it. I'm a bit confused here because that's somewhat easy to compute and if p was a cumulative distribution here, I'd understand this one right away.

So basically I'm stuck with understanding this one and looking forward to any hints as to where I'm wrong.

Thank you folks in advance!

I agree with you, this looks like a flawed description of an inverse transform generator.

The $p[i]$ shouldn't be like described, but should be the prefix sum of these, corresponding to the cumulative distribution.
The $\arg\min$ should indeed be taken on the absolute value (or allow positive values only).
As the pdf is computed on the positive $i$ only, the range of the uniform drawing should be $[0,\frac12]$ (or the $p[i]$ scaled by $2$.)

(When the $p[i]$ are not cumulated, they are in $[0,\frac1{\sigma\sqrt2\pi}]$, and the range $[0,1]$ doesn't make sense.)

For better efficiency, you can precompute a lookup-table of the $i$ indexes for sufficiently dense $q_1$.


Any way to reverse engineer professional audio algorithms or those in e.g. rack units? Any way to reverse engineer professional audio algorithms or those in e.g. rack units?

In order to study how they're implemented?

Not suggesting to steal ideas, but as education.

Not really.

Anything that's a linear and time invariant is of course easy to characterize but most pro audio stuff doesn't fall into this category. It's typically dynamic, time variant, non-linear, signal dependent, etc. These types of algorithms are extremely difficult to quantify if you don't know the internal structure and the parameter space.

Even if you get a full characterization, you may still not learn something about the way it's done. A time-invariant reverb, for example, can be measured through it's impulse response. However, you can't really tell how it's generated (unless it's a convolutional reverb in the first place).


Can a discrete sinusoidal signal represent only limited number of frequencies? If $\omega$ is the frequency of a continuous time sinusoidal signal, $-\infty < \omega < \infty $.

But when it comes to discrete time sinusoidal signals, the range becomes limited : $-\pi \le \omega_{discrete} \le \pi $.

Does that mean discrete sinusoidal signals cannot represent all frequencies possible ? What is the true picture here ?

It is important to realize that in the discrete-time case, $\omega$ is normalized by the sampling frequency $f_s$:

$$\omega=2\pi f/f_s$$

The inequality in your question is basically the sampling theorem, which says that a discrete-time signal can only represent frequencies lower than half the sampling frequency. All higher frequencies will be folded back into the frequency band $[0,f_s/2]$ (this is called aliasing). If you need higher frequencies then you need to increase the sampling rate.


Automatical skyline alignment on image I need for an algorithm which do steps, displayed on image below.
Are solutions for doing this exist? (Preferably on C# or JavaScript).

You could try detecting line segments in the image for example using hough transform, or maybe lswms, edlines or lsd.

After detecting the lines, you can compute their length and angle and find the longest almost horizontal line. If you are lucky, this line will be the horizon and you can rotate the whole image in the opposite direction to get the aligned image.

This might work but it will depend on how good the horizon is visible and how well it is detected.

Maybe you could find a number of potential horizon lines and have the user select the best one.

On second thought, if you involve users into this, just let them draw the horizon and then rotate the image. That's not really what you were asking for, but it will probably work better than what I suggested above. :)


Why Shannon Theorem has nothing to do with frequency? According to Shannon theorem
$$C= B \log _2(1 + S/N)$$
My question is, why there is no frequency in this formula?

Let's say we have 2 channels, both of them have same 20 MHz bandwidth,
one is from 100 MHz to 120 MHz,
another is from 1 GHz to 1.02 GHz.

Does this formula mean that these two channels have the same transmission speed?

Cause I think the one from 1 GHz to 1.02 GHz is faster, since it has a higher frequency.

Does this formula mean that these two channels have the same transmission speed?

Yes. That's exactly what you should take away from this: Channel capacity has nothing to do with center frequency; only bandwidth and SNR are relevant.

  Cause I think the one from 1 GHz to 1.02 GHz is faster, since it has a higher frequency.

But that's plain wrong. It's not based on an understanding of what it means to use the channel:

To transport a specific amount of data, you will need to transmit something; the faster you change what you're transmitting, the more data you can send per second.

But the more often you change what you're transmitting, the higher the bandwidth. Whether that happens in a channel around 100 MHz or around 72 GHz makes no difference whatsoever.


what is standard acoustic volume In Audio Experts book, it's given as

  Sound pressure meters must be calibrated at factory to a standard acoustic volume.

So what doth this standard acoustic volume means?

Normally it is 94 dB of Sound Pressure Level, where 0 dB is 20 µPa of acoustic pressure. This is the SPL you will find in most of the acoustic calibrators for measurement microphones.


What pattern maximises a JPEG file's size? As I understand the JPEG algorithm, it smooths (ignores) gradual variations in colour and brightness in favour of high frequencies.  This loss of information is what creates the characteristic high level of compression relative to lossless algorithms.

So the grey square would result in a small JPEG file, whilst the grids should generate larger files:-

Clearly such grids /patterns could be of any pitch right up to single  pixels.  I'm stipulating that all the images are grey scale types. Is there some way to determine the type of pattern to maximise the size of a JPEG file, or must this be done experimentally?

Pursuing on @OlliNiemitalo, JPEG allows some loss related to frequency quantization. So any file can be compressed to any size, hence you cannot maximize a JPEG size without setting some constraints, like quality. Let us assume a standard metric, like mean-squared error.

Now, let us check some of the JPEG features for image sparsification, and let us try to go the opposite direction:

JPEG assumes similarities in RGB color planes, correlated with a luminance/chrominance transform for instance. So you choose independent R, G, and B, you are good
8$\times 8$ patches are compressed almost independently, except for the average, predicted for the blocks before (raster order). If you have neighboring blocks with very different means, good
the more the high-frequency, the less a 8$\times$8 block as natural frequency, the better. They avoid long trails of zeros encoded by an EOB symbol.

One cannot easily fulfill all these properties together,  especially with RGB planes made of 8-bit integers.
Forgetting about the average-thing, generating  independent RGB images, with each pixel randomly $0$ or $255$ (ie high range). Quantization is likely to produced here a lot of false colors. This could give you very bad images for the JPEG in MSE. However, original and compressed ones could look as ugly.

You can find other fractal-based or fancy-patterned proposal in
What is the least JPG-compressible pattern? (camera shooting piece of cloth, scale/angle/lighting may vary)


What proportion of a padded FFT should be actual values For a given signal, I've been told that you can pad the vector with $0$s at the end to get a larger DFT, and as a result get more precision in frequency bins.

What are the limits to this approach?
Is there a loss of accuracy when you do this?
Or if not, why wouldn't you just use a small number of points and pad with zeros all the time?

Be careful with this in thinking that you would increase your frequency resolution- you won't! Zero padding is very effective in iterpolating more samples between the samples you have, but it does not add any more information about the frequency content of your signal between your original bins.  You do not lose anything either.

The frequency resolution in all cases is $1/T$ where T is the time length of your data.

Please see my answer at this post specifically for more details on zero-padding and how it works: What happens when N increases in N-point DFT


Confused about the period of a discretised time-series I have managed to confuse myself about the periods of discretised time series. It is all rather embarrassing.

This is one period extracted from my continuous signal (radians on y, time on x):

To use it in my calculations I use a discretised vector of N points. The above signal has a period of $T=0.95s$. The discretised version A has $N=100$ and the angle, over time, is indexed from 1 to 100.

Q1: is the period of my angle vector $[\theta_1,\ldots,\theta_{100}]$; $T=100$?

Now suppose I need to standardise my data as such, incl the index:

% standardise training inputs
m = mean(A);
s = std(A);
x(:,1) = 1/s(2) * (A(:,2) - m(2));  % index
x(:,2) = 1/s(1) * (A(:,1) - m(1));  % velocity

Q2: Is the new period then simply max(x(:,1))?

On the continuous setting, the period is the minimum $T\ne 0$ such that $\theta(t+T)=\theta(t)$. Once you decide to have $N$ regularly sampled points inside a period, you introduce a sampling period $T_N$ that is directly proportional to the period: $T_N=T/N$. Hence, you can describe your period as a quantity with a time unit (second) or as a unitless number of samples, equivalently.

Now, if for some reason you normalize the time axis with an affine function, the mean will only shift the whole axis (no change for the period), and the scaling only will dilate or shrink the interval between samples.

So here, the effect of  $s(1)$ will turn $T_N$ into $T_N/s(1)$, and also  $T$ into $T/s(1)$.

In your case, you can directly compute $s(1)=T_N\sqrt{\frac{N(N+1)}{12}}$. In this formula, the scale factor depends on $N$ only, and not on the start and end indices.
For $N=100$, the new period should be now a unitless quantity of
$$T/s(1) = \sqrt{\frac{12N}{N+1}}\,.$$


Easy way to get rid of noise in a hand drawing I would like to take a user's hand drawn input, and filter out "noise", aka small sketches or blobs, where the user might have screwed up and accidentally touched his hand to the drawing device, or something similar.

I'm using Imagick, a PHP wrapper that wraps an image processing library that does everything short of actual, high-level blob detection.  So I'm free to use convolution, erosion, dilation, and a slew of other image processing tools.  I can't write blob detection in the high level PHP code because it's too slow.  And before I go down a path of writing the C code using ImageMagick's API, and wrapping it in PHP, I was wondering if there's some sort of "trick" that I can do to accomplish my objective.

For instance, my initial thought would be to do an erosion on the image (in normal situations, when I think of removing noise from an image, I think to erode then dilate).  But since ALL the lines in the image are just a few pixels wide, everything would disappear.  Then I was considering some sort of blur so the longer lines would end up with higher intensities, then I could just threshold out the lower ones.  But again, in the context of lines, I'm not sure how well that works.

I'll attach an image of increasing "blobs".  I would consider the first 4 left ones as blobs that I want to remove, but then keep the 5th, the line, and the box.

You can use opening on your image http://php.net/manual/en/imagick.morphology.php with a kernel big enough to erode the small objects


Determine filter given transfer function I am given $$H(z) = 1 + \frac{\alpha}{1-\alpha z^{-1}}$$ where alpha is between $0$ and $1$. This is apparently is a low-pass filter with cutoff frequency $f_c$.

How can I see that? And how can I compute/plot the frequency response in Matlab? Neither freqz(1+alfa,[1 -alfa]) nor abs(H(z)) seems to yield a low-pass filter on the above transfer function in Matlab.

Furthermore, $1-H(z)$ should be a high-pass filter. But again, using freqz(-alfa,[1 -alfa]) in matlab it is clearly not a high-pass filter.

No $$H(z) = 1 + \frac{\alpha}{1-\alpha z^{-1}}$$ is not a lowpass filter. It's a lowpass filter in parallel with an all pass filter (kind of low-boost shelving filter). If you want a lowpass filter which rejects the high frequencies then you may take the following:

$$H_{lp}(z) = \frac{\alpha}{1-\alpha z^{-1}}$$ which is a lowpass filter for the range of parameter $\alpha$. And therefore $1 - H_{lp}(z)$ will be the complementary highpass filter in simplest terms.

So you can see it by the following code:

freqz(alfa, [1 -alfa])

Note that the gain of the filter will be dependent on $\alpha$! Note also that as Matt.L commented you might have $1 + \alpha$ in the numerator instead...


Audio Recognition and Learning I'm a beginner with absolutely no experience in signal processing, but it has been a long-term goal of mine to do something that mixes music (a hobby of mine) and computer science (my "expertise").  I've been thinking of doing a project that involves learning notes from music, but I don't know where to begin.

For one, there doesn't seem to be a name for this problem. However, it doesn't seem to be a solved problem.  I've searched terms including "Audio recognition," "note recognition," "chord recognition," but I've only come up with this paper: "Automatic Chord Recognition from Audio Using an HMM with Supervised Learning" by Lee and stanley.  It seems to have some good references to processing and even a pointer to a dataset I could potentially use, but one paper doesn't seem to be enough for a good understanding.  Alternatively, I could skip the processing part entirely but I can't find data sets using terms like "audio signal data set"

Signal processing seems to be a pretty broad field. To save time, I want to pick up a book more oriented towards my project.  Are there any books like this?  I've found this:
http://www.dspguide.com/ch1/3.htm but I want to know if this is relevant before I jump into this 34 chapter book.  Are they any other search terms I can use such as "audio feature extraction" to get me closer to the ML part of the problem?

Finally, to cut down on work, I want to use libraries for the signal processing part.  One library I found for Python is Pyo but are they any other alternatives for Python that people have used?

Basically, I'm a little overwhelmed by how much I don't know and I would greatly appreciate a pointer in the right direction.  Even google search words would help.

The Lyons book is gold for someone in your situation.  Have patience with the topic of audio DSP -- the more time you spend trying to learn it, the less you will know that you know.  :)


How to demodulate an OFDM signal I'm looking at decoding an OFDM signal which consists of 6 carriers (or tones) which are BPSK modulated and a pilot tone which aids tuning. This is the first time I have worked with OFDM so I need to know if I'm approaching this in the right way.

The way I'm thinking of decoding it is to use the pilot tone to calibrate (as the receiver may be slightly mistuned) then to have six band pass filters to separate out each carrier which is then demodulated in the usual way. Can anyone see any problems with this ? or can you suggest a better way of doing this.

One of the nice features of OFDM is that it allows a very simple structure for the modulator and demodulator: given a set of symbols (in general complex-valued, taken from a signal constellation such as BPSK, QPSK, or QAM) to map onto each carrier, the modulator can be implemented using an inverse discrete Fourier transform, typically implemented using an FFT. Each set of symbols (one per carrier) is transformed to yield an OFDM symbol, which is then sent to the channel. The DFT length will typically be chosen to be larger than the number of desired carriers to allow some "guard band" near the system Nyquist rate.

In addition to the above DFT-based structure, most OFDM systems also incorporate a cyclic prefix, which allows for simple implementation of an equalizer in the frequency domain. Equalization can provide improved link performance in multipath environments (e.g. many wireless communication scenarios). It can also be used to aid in synchronization, as described below.

The simple structure carries over to the receiver; an OFDM waveform can be demodulated using the inverse transform to that used at the transmitter, yielding the original symbol values. The inverse to the inverse DFT used at the transmitter is a "regular" (forward) DFT. Thus, you will often see OFDM receivers diagrammed with an "FFT" block at the front end. The output of the transform contains the symbol values mapped onto each of the carriers, including all of the unused ones that make up the guard band. The demodulator plucks out the (complex-valued) amplitudes of each of the carriers of interest and passes those on to any further decoding logic (equalization as described above, channel decoding, mapping to bits, etc.).

As usual, however, the answer isn't quite so simple; the above explanation overlooks some important issues that must be addressed for a practical system:

Timing synchronization: When you actually think about how you would build an OFDM receiver, one of the first problems you will run into is how to align the receiver's FFT frame with the stream of incoming samples. Synchronization with the OFDM signal's symbol timing is required to properly align the receiver's FFT operation with the appropriate time period in the observed sample stream.

This can be implemented using a correlation-based approach. As stated previously, most OFDM waveforms include a cyclic prefix, which is a scheme of forceably adding some circular periodicity to the transmitted waveform. This can be exploited at the receiver to obtain symbol timing; the timing detector simply calculates the sliding autocorrelation of the observed symbol stream using a lag commensurate with the known period between the transmitted signal and its cyclic copy. The magnitude of the result will reach a peak at the instant that corresponds to the start of each OFDM symbol.
Frequency synchronization: Fine frequency synchronization is also key to robust OFDM reception, as frequency error causes intercarrier interference. Correction of frequency error can also be estimated using the timing synchronizer's correlator output. As stated previously, the autocorrelation of the observed stream at a lag equal to the cyclic prefix delay has a large magnitude at the start of each OFDM symbol. The phase of the correlator output provides a measure of the amount of phase drift over the course of each symbol time. This measure of "phase drift per unit time" can be recast as a measure of "frequency drift" instead. If the receiver can safely assume that the frequency error is constant over the course of a symbol time (which is reasonable for many cases), then the bulk frequency offset can be removed before calculating the DFT.

There can be even more problems to tackle for each of your carriers, depending upon the modulation used on each. For the simple case of BPSK, you may also need to worry about phase synchronization if you desire a coherent receiver. However, synchronization of timing and frequency are the key implementation details that often seemed to be glossed over in discussion of OFDM receiver structures.


Fast Fourier Transformation using parallel processing coding for image processing I'm now working in parallel precessing -and i want a full code for fast Fourier transform(FFT)
written  parallel processing Code for any image processing application

Well, you should have a look at FFTW as their implementation allows you to use more than one thread to perform the FFT. One thing that should note with this library is that multi-threading is not necessarly beneficial for all problems sizes as stated here. For a software application that I'm working on, for example, using multiple threads slows down things (the biggest FFT being on 4096 points).


The plot of instantaneous power of the Dirac function I am very confused. I have tried researching this question for the last two weeks and I cannot get a conclusive answer.
I was wondering how would I go about plotting the instantaneous power in the time domain of the Dirac input signal.
The problem becomes as I try to calculate the power in time domain, $P(t) = \left\lvert\delta(t)\right\rvert^2$. I know how to calculate power by using Rayleigh's theorem to convert the the signal to frequency domain where delta(f) = 1. to get infinity over an infinite period which is undefined value for power, but I just want to plot power over time for the Dirac function.
delta isn't square-able since as x goes to 0 the delta function goes to infinity and you cannot square infinity, also something about the delta function being a distribution and not able to be square-able.
I was wondering if someone could point me in the right direction to be able to plot this function.
Thank you

The problem becomes as I try to calculate the power in time domain, $P(t) = \left\lvert\delta(t)\right\rvert^2$

As you figured out yourself, that's impossible: the power of the Dirac impulse isn't defined, it diverges (goes to infinity).
Hence, you especially can't plot the power "correctly"; it's zero anywhere but at $t=0$, where it diverges.
But quite honestly, this should come as no surprise: you can't correctly plot $\delta(t)$, either (the value at $t=0$ isn't defined just as much as its square isn't).
So, whatever question can be answered by that plot: I'm afraid you'll have to find a "better" question.


What kind of noise is city noise? White, pink, brown...? In general, to what category does it subscribe the spectrum of city noises, like traffic noise, ambulance/police siren, alarms, etc, put together?

If the question is too broad, then, what is the color for particular noises, like traffic?

Noise color isn't a rigorous model for traffic noise

https://hal.archives-ouvertes.fr/file/index/docid/506615/filename/Can_2010_Applied_acoustic_71_8_p.pdf


Enhancement by adaptive wavelet shrinkage denoising I was wondering does anyone know any programming library or a way that can help me implement this algorithm.

The cropped image stems from slides named: Introduction to Medical Image Processing. It illustrates to me edges emphasis and/or contrast enhancement (unsharp masking) by non-linear coefficient weighting in appropriate wavelet domains.
I would retrace the mammographic image and  enhanced versions  in papers by Andrew Laine and co-authors, like Enhancement by Multiscale Nonlinear Operators (Andrew Laine, Walter Huda, in Handbook of Medical Imaging, 2000), Contrast Enhancement by Multiscale and Nonlinear Operators (J. Fan and A. Laine, 1996) [Figure 7.14: (a) Oringinal mammogram image M73. (b) Nonlinear enhancement with adaptive wavelet shrinkage denoising]:

The bottom plot is unknown to me. I don't have the knowledge of such codes  shared by the authors. There is a recent code called Wavelet Unsharp-Mask and Clarity. You might look for "homomorphic wavelet filtering" that could provide similar results.


Controllable realisation of $\frac{s^4+1}{4s^4+2s^3+2s+1}$ is both controllable and observable? I am trying to find the controllable realization of the following transfer function:

$$H(s)=\frac{s^4+1}{4s^4+2s^3+2s+1}$$

I approach this by first using polynomial division to assure that $H(s)$ is strictly proper. This yields:

$$H(s)=1/4+\frac{-1/8\cdot s^3-1/8\cdot s+3/16}{s^4+1/2\cdot s^3+1/2\cdot s+1/4}$$

From this I can extract the controllable canonical form where the matrices A, B, C and D are:

$$A=\begin{bmatrix}0 & 1 & 0 & 0\\ 0 & 0 & 1& 0 \\ 0 & 0 & 0 & 1\\-1/4 & -1/2 & 0 & -1/2\end{bmatrix}$$

$$B=\begin{bmatrix}0\\0\\0\\1\end{bmatrix}$$
$$C=\begin{bmatrix}3/16 &-1/8 &0 &-1/8\end{bmatrix}$$
$$D=1/4$$

From this I am confused because when I construct the observability and controllability matrices I found that the system is both controllable and observable (both matrices have full rank). This seems like a contradiction. Can it be possible? What am I doing wrong?

Controllability is a property of $A,B$ and observability is a property of $A,C$  matrices.

Transfer functions don't exhibit this property since if there is such a problem then you will not get a nonzero transfer function as there will be no possibility of control so nothing to transfer or nothing to observe. Hence, this is also a state space property.

The main theorem says that a state space realization is minimal if and only if it is both observable and controllable.

These representations allows one to quickly see, if exists, the unobservable and uncontrollable parts of the system. Such representations are obtained via the procedure called Kalman decomposition.

Here you have only constructed the controllable canonical form and already found out that the system is also observable. Thus the represantation is minimal.


Meaning and unit of frequency in Laplace (Fourier) transform Imagine transfer function obtained by Laplace transform, for example:

$G(s) = \dfrac{1}{s+1}$

Now, I would like to do some frequency analysis, so I replace the $s$ with $\omega i$ (let's consider this operation valid for this example).

What is the unit of the $\omega$? So far what I have seen, the $\omega$ is noted as frequency or angular velocity. I asked my colleagues and I got various answers:

rad/sec
Hz
no unit

What is correct and why? Does it depend on real variable passed to transform (if somebody uses different variable than time)?

If you are dealing with the Laplace transform $G(s)$ of a time domain signal $g(t)$ and its evaluation on the imaginary axis to get the Fourier transform $G(j\omega)$ (assuming it exists) then the unit of your frequency $\omega$ is radians per second assuming the unit of the time was seconds.

It's relation to cyclic frequency is : $$\omega = 2 \pi f$$ where $f$ is the frequency in Hz (cycles per second).

On the other hand if the initial function was like $g(x)$ where $x$ was a spatial variable with unit of meters, then the transform domain frequency unit will be in radians per meter where its relation to space-frequency is still the same with $\omega = 2 \pi f$ where $f$ will have the unit of cycles per meter