Cost function for adaptive algorithms I am having little difficulty to understand that why most of the adaptive algorithms use error power or addition of error power as cost function/minimization criterion. I have read that minimization of error power criterion gives us best results from estimation mean point of view. I am not being able to grasp that either. So can someone please tell me in simple language what are the actual reasons to use above criterion? Why cant we simply use error as the minimization criterion?

In a very simple form: Gaussian-like noises are frequent, and derivatives of squared functions are easy to minimize.

In very simple words: error, directly, cannot be used as a criterion. You can minimize a signed function to $-\infty$, which is not really useful here. As least, you need a cost function with a lower bound.

Power (if by that you mean the squared error) is an increasing function of the error, and positive. So if you minimize power, you can expect to minimize, somehow, a measure of the error (more precisely its absolute value).

The first problem is that the absolute value is not differentiable around zero, and can be difficult to optimize (one may look for a zero in the derivative). And when can you very easily find a zero in a derivative? When it is linear, because since Laplace, Gauss, Seidel such systems can be solved efficiently (matrix inversion). What are the functions that yield linear derivatives? Squares, parabolas...

Basically, for a while, squared functions were the only ones one could solve easily. The standard mean is just the minimum of the sum of differences squared. Interestingly, in linear cases, unrelated disturbances can add with a Gaussian distributions. And Gaussian distributions are well-related to squared errors.

But you are interested, recent advances in optimization allow to optimize more easily other types of functions, including absolute value loss criteria   or penalties.


What technique could I use to smooth the right amount these signals? I'm developing an application which one of its capabilities is peak detection.
I expect to handle three types of the similar data (Acceleration vs Frequency) but measured in different laboratories and with different sampling and noise. The data is not equally spaced in frequency.

Here is the plot of each of them:

Detail of the previous one

The first contains about 2000 points and there are a lot of local maxima detected. I can't use an average-smooth because the last one contains only 100 points and if I do so I will be missing important peaks (each peak occupy only 3-7 points).

What can I do to improve the detection of peaks and how it could be done?

Thank you very much.

If you come up with a set of criteria that define what makes a local maximum become a "peak of interest", then you could test whether each local maximum meets the criteria you set. For example:

A local maximum at $X(f)$ is considered a peak of interest if for some chosen small value $\epsilon$ and some chosen height $\mu$ , we can satisfy the following

$$\exists \, \delta_1 , \delta_2 \lt \epsilon \quad s.t. \quad X(f)-X(f-\delta_1) > \mu \quad \text{and} \quad X(f)-X(f+\delta_2) > \mu$$

If one such criterion is not sufficient to describe your concept of a peak of interest, then you can logically combine multiple criteria.

You could approach the problem from a filtering and smoothing perspective as you mention, but notice that each filter or smoothing operation has an implicit criteria that it is enforcing.  No matter how you implement your peak detection, realize that there is no "right answer" as to what defines a peak of interest. So before you can achieve the desired result you must first decide what that result should look like.


RF signal identification algorithms There exist many, many analog and digital data transmission protocols. It is often difficult to tell which mode one listens to and tedious to do for many signals.

In Ham Radio, a Reed-Solomon Identifier may precede the signal, but this is not the case for every signal in existence.

Are there algorithms that automate the detection and mode differentiation from raw signal data that don't include metadata such as RSIDs?

Plenty! In fact, since the mid-1970s, people in security have been working on automatic transmission classification.
It's a hard problem, and you logically can't detect all transmission schemes, but it's something that's being sold as service and as product, especially for government (and three-letter government agencies).

These guys, for example, focus on HF/VHF/UHF comms, and I've seen the list of classifiable transmission schemes they had for a public demo. It's long.

Now, what these recognize are just common standards (in Ham Radio, that would be called modes),  and they are specific to one greater usage scenario (security-relevant communication).

There's a lot of underlying algorithms to automatically classify the number of observed channels, bandwidths, symbol rates, constellations, spreading, scramblers, error-correction codes, frame lengths, synchronization preambles.... but even listing these would take a book.

I can point you to what a student did during Google summer of code: It's gr-inspector, a GNU Radio-based signal detection and classification program.

It's pretty generic in that it detects bands with power in them (which is no good for spread-spectrum systems, but for anything else it's a viable approach), and it gets pretty specific once you have an OFDM signal in there. Part of his job was writing a classifier just to get OFDM parameters: carriers, spacing, length of cyclic prefix. That took someone who was dangerously close to being a M.Sc. in communications engineering nearly 3 months. You can imagine how much work it is writing a set of classifiers for all useful transmission schemes!

You might be interested in the two papers he cites in his blogpost I linked to above; and especially their references.


How to scale Phase Transition Diagram for Compressed Sensing? I want to compute a Phase Transition Diagram as shown here ($A \in \mathbb{R}^{n \times N}$ and $k$ is the sparsity: $\vert \vert x \vert \vert_0 = k $ )

My question is: For $n=1$ I can only compute $k=1$, for $n=2$ I can only compute $k=1,2$ and so on. Therefore the computable number of $\rho$ values is different for each $\delta$ value, but in the plots it is even for $\delta=0$ a continuous, changing plot into the y-direction.
My plots are currently looking like this (yellow means reconstruction with probability 1):

I recommend first fixing $N$. By the way, remember here that compressed sensing works asymptotically well; for $N \rightarrow \infty$ your phase transition is going to move up/left (better) in the $\delta, \rho$ diagram, see for example Donoho & Tannner 2010.

Your size of $N$ essentially determines the resolution you can achieve in your phase transition diagram. For example we can look at Maleki & Donoho 2010, Section IV: they describe how they select $N = 800$ and then calculate $n$ and $k$ from a grid of values $(\delta, \rho)$ between 0.05 and 1 in 30 steps. For each value $\delta$ or $\rho$, the corresponding $n$ or $k$ are calculated as $$n = \lceil \delta N\rceil,$$ respectively $$k = \lceil \rho n \rceil.$$

As you can see there is generally some imprecision involved during to the rounding. If you select a too fine grid of values $(\delta, \rho)$ for a given $N$, several of the resulting $n$ or $k$ values are going to round to the same number. You will have to either select a larger $N$ or choose a coarser grid for your $(\delta, \rho)$ to avoid this.

Finally, you should (if you want the same type of plot as Donoho et al.) plot your success ratio in each of your points in a $(\delta, \rho)$ coordinate system. Your figure above looks as if it was perhaps plotted with $k$ on the first axis?

Now, if you wish to estimate the actual phase transition boundary, you can also see in Maleki & Donoho 2010, Section IV, how they fit a logistic curve to the success ratios and can then read a specific percentage reconstruction boundary $\rho$ for each value $\delta$ from the fitted logistic function.


What is difference between terms $X(j \omega) ,X(\ e^{j \omega })$ and $ X(\omega)$? While studying frequency transforms ,I get confused with the terms like $X(j \omega) ,X(\ e^{j \omega })$ and $ X(\omega)$
 ,where $ \omega = 2 \pi f $.
 So what is the difference between them ?

they should have different subscripts.

bilateral Laplace Transform:

$$ \mathscr{L} \left\{ x(t) \right\} \triangleq X_\mathcal{L}(s) \triangleq \int\limits_{-\infty}^{+\infty} x(t) e^{-st} \ dt $$

continuous Fourier Transform:

$$ \mathscr{F} \left\{ x(t) \right\} \triangleq X_\mathcal{F}(\omega) \triangleq \int\limits_{-\infty}^{+\infty} x(t) e^{-j\omega t} \ dt $$

it's pretty easy to see that

$$ X_\mathcal{F}(\omega) = X_\mathcal{L}(j \omega) $$

likewise in discrete-time, the bilateral $\mathcal{Z}$ Transform:

$$ \mathcal{Z} \left\{ x[n] \right\} \triangleq X_\mathcal{Z}(z) \triangleq \sum\limits_{n=-\infty}^{+\infty} x[n] z^{-n} $$

and the Discrete-Time Fourier Transform (DTFT):

$$ \mathcal{DTFT} \left\{ x[n] \right\} \triangleq X_{\small{DTFT}}(\omega) \triangleq \sum\limits_{n=-\infty}^{+\infty} x[n] e^{-j\omega n} $$

it's also easy to see that

$$ X_{\small{DTFT}}(\omega) = X_\mathcal{Z}(e^{j\omega}) $$

now often the context is understood and everyone knows that $X(e^{j\omega})$ is the Z transform evaluated at $z=e^{j \omega}$, which is the DTFT and we don't bother with the clunky subscript notation.  we might leave the space for subscripts for another function, like if we have a vector of signals as with a state-variable system representation in control theory.


Doubling the frequency of a sound fragment in real time? I'm given given a stream of sound, as a serie of little buffers .

Is it possible to make this stream sound one octave up ? Or an arbitrary number of semitones up ? My goal is to add it to the original sound and to pass this through a reverb to create a "shimmer" effect .

My research so far was in the direction of doubling the frequency by creating a new buffer, and copying every other sample of the old one. This works, but I have a lot of noise and artefacts. I guess this is because the buffer length is not an exact multiple of the frequency but I can't change the buffer length and I have no prior knowledge of the frequency of the sound.

I don't think this is the place to post C/C++ code so I won't do it but I guess you get my idea.

I'm not looking for a high quality sound at all, i.e, a perfect, one octave up, reproduction of the sound,but at least something decent.

Any ideas please ?

After loads of research, I found a really good article about Short Time Fourier Transform (STFT), which can apply FT to a really small buffer, change it's frequency and apply IFT. There's even an implementation in C. You can find it here : http://blogs.zynaptiq.com/bernsee/pitch-shifting-using-the-ft/

As a "side effect", you also know the pitch of your signal in almost real time

HTH


MATLAB center frequency shift I am wondering how to shift the center frequency of the below complex data from 0, to say 20 MHz (10 MHz to 30 MHz- 20 MHz bandwidth).  I tried to change the indexing of the bins but I could not get it to shift.  I understand fftshift shifts frequencies about 0 Hz, any way to shift the whole spectrum to something other than 0?

nFFTSize = 64;
subcarrierIndex = [-26:-1 1:26];
nBit = 1000;
ip = rand(1,nBit) > 0.5;
nBitPerSymbol = 52;
nSymbol = ceil(nBit/nBitPerSymbol);

ipMod = 2*ip - 1;
ipMod = [ipMod zeros(1,nBitPerSymbol*nSymbol-nBit)];
ipMod = reshape(ipMod,nSymbol,nBitPerSymbol);

st = [];

for ii = 1:nSymbol

inputiFFT = zeros(1,nFFTSize);

% assigning bits a1 to a52 to subcarriers [-26 to -1, 1 to 26]
inputiFFT(subcarrierIndex+nFFTSize/2+1) = ipMod(ii,:);

%  shift subcarriers at indices [-26 to -1] to fft input indices [38 to 63]
inputiFFT = fftshift(inputiFFT);

outputiFFT = ifft(inputiFFT,nFFTSize);

% adding cyclic prefix of 16 samples
outputiFFT_with_CP = [outputiFFT(49:64) outputiFFT];

st = [st outputiFFT_with_CP];

end

st2 = st*exp(2*pi*20*i*.0000032);

scatterplot(st2), grid;

You can shift your spectrum in the time domain by multiplying it with a complex exponential. It is the same way you implement an AM modulator.

Let $A(t)$ be an unspecific baseband signal. You multiply in time domain with a complex exponential and the result will be the bandpass signal:
$$
y(t) = A(t)e^{j2\pi ft} = A(t)*\cos(2\pi ft) + j*A(t) \sin(2 \pi ft)
$$
Since cos and sin both contain only a single frequency component, a fourier transform (if you look at the real part) will result in two impulses at $\pm f$ and the spectrum of $A(t)$ shifted to $\pm f$.

The same applies to your problem: Multiply your time domain signal with a complex exponential of the desired frequency offset and you will observe your spectrum to be shifted down.


Optical flow equation and real images We know optical flow equation cannot be calculated without additional constraints. But what I cannot figure out is why it does not hold exactly for real images? This equation has gradient and so, it can be used to estimate the motion vectors by focusing on edges. Is it correct? How we add other constraints like smoothness into account?
Thanks

Writing the optical flow vector as $(u,v)$, then the optical flow equation only states that light is preserved.
Basically, it states that variations in the observed light patterns (measured by the image derivative with respect to time) are explained by the motion of objects (smaller light patterns inside the image).
Object motion is measured through spatial variations of the said light patterns:

$$\frac{\partial I}{\partial t} + u\cdot\frac{\partial I}{\partial x} + v\cdot\frac{\partial I}{\partial y} = 0$$

This problem is ill-posed for any image, be it synthetic or real.
This is easy to see since 1 pixel provides 1 measurement $\frac{\partial I}{\partial t}$ only, while at the same time there are two unknowns $\left(u(x,y), v(x,y)\right)$.

As you said, since this equation involves image spatial gradients $\nabla I = (\frac{\partial I}{\partial x}, \frac{\partial I}{\partial y})$ it can be solved for edges (by tracking them).
Other ways to overcome this issue:

assuming that the optical flow is constant over a square of several pixels and then solving the resulting system. Typically, using 3x3 squares yields 9 equations and a least squares solution to the optical flow can be computed;
more involved approaches rely on regularized problems. Regularized approaches look for a flow that will fit the observed data while having good properties such as smoothness, exhibiting thin discontinuities...


Why is the state space representation, $dx/dt = A x(t) + B u(t)$ linear? One thing puzzles me about the state space representation and that is why is it claimed that such systems are linear and time-invariant but particularly why are they claimed to be linear?  The right-hand side of $dx/dt = A x(t) + B u(t)$ does not seem to obey superposition? Is it because the solutions to these systems are linear?

For example:

$y_1 = A x_1 + B u $

$y_2 = Ax_2 + B u$

$y_3 = A (x_1 + x_2) + B u $

Therefore $y_3 \neq y_1 + y_2 $

I'm I doing something wrong?

What you have listed is a state equation.  Usually a state-space representation has two equations:

(i) $\dot{x}(t)=Ax(t)+Bu(t)$
(ii) $y(t)=Cx(t)+Du(t)$

One takes the input $u(t)$ and generates the state $x(t)$ and another that takes the input $u(t)$ and state $x(t)$ and generates the output $y(t)$.  You are missing the second equation.

That aside, the equation that you have is basically a 1st order linear differential equation.  It's important to note, here, that the input is not $x(t)$, it is $u(t)$.  $x(t)$ is the state variable.

I will examine just the first half (input to state) and leave the second half (input and state to output) to you.  I will remark, though, that if $x(t)$ is linear with respect to $u(t)$ then $y$ in equation (ii) above necessarily has to be linear by virtual of adding two linear terms together, so the second half is really trivial.

The first half is:

(1) $\dot{x}(t)=Ax(t)+Bu(t)$

Can you show that it is linear?

I will show additivity and leave homogeneity to you.  Say you have inputs $u_1$, $u_2$, and $u_3=u_1+u_2$.

One way to do this is to solve (1).  As it is a 1st order ODE, it is not too difficult, but let's be lazy.

Let $x_1$ be the state generated by input $u_1$ and $x_2$ be the state generated from input $u_2$, then we have

(2) $\dot{x_1}(t)=Ax_1(t)+Bu_1(t)$

and

(3) $\dot{x_2}(t)=Ax_2(t)+Bu_2(t)$

adding (2) and (3) gives

(4) $\dot{x_1}(t)+\dot{x_2}(t)=Ax_1(t)+Ax_2(t)+Bu_1(t)+Bu_2(t)$

From linearity of derivative, we get that:

(5) $\dot{x_3}(t)=A(x_3(t))+B(u_3(t))$

where $x_3 = x_1+x_2$.

So input $u_3 = u_1+u_2$ generates state variable $x_3 = x_1+x_2$, showing that additivity holds for the state equation.


Does an equivalent transformation of a signal to a spectrogram image exist in which the phase information is part of the resulting image? I'm working on a research project where we would like to apply convolutional neural networks to an image representation of a signal. However, it seems that if I would use a spectrogram, I would end up loosing the phase information. As a result, my question is:

Does an equivalent transformation of a signal to a spectrogram image exist in which the phase information is part of the resulting image?

Thanks a lot and have a great day, Maxime

From https://en.wikipedia.org/wiki/Spectrogram:

  Limitations and resynthesis

  From the formula above, it appears that a spectrogram contains no
  information about the exact, or even approximate, phase of the signal
  that it represents. For this reason, it is not possible to reverse the
  process and generate a copy of the original signal from a spectrogram.

If you create a volumetric (unflattened 3D) image, you can use 2 layers in the 3rd dimension to represent magnitude and phase, or real and imaginary components of a complex spectrogram output.

In a 2D color spectrogram, you could try using an orthogonal color mapping, for instance, in RGB space, red + green for magnitude and blue (or delta blue) for phase, or similar for real and imaginary components, or use YUV or HSL color spaces, etc.


periodicity of constant discrete time signals are constant discrete time signals periodic?
example \begin{equation}
e^{i10\pi n}
\end{equation}
my proffesor says that this signal is aperiodic, in the discrete sense. but it seems wrong, because
unlike in the continuous case, i can calculate the smallest time period , which is 1.

normally, "$n$" is the symbol we use here for discrete-time.  if your professor said that:

$$\begin{align}
x[n] &= e^{i10 \pi n} \\
 &= e^{i 2 \pi (5n)} \\
\end{align}$$

is not periodic with a period of $1$ (assuming $n \in \mathbb{Z}$) or a period of $\frac15$ (assuming $n \in \mathbb{R}$), then your professor is mistaken.


Can someone help explain what this FFT workflow is doing to my signal, and why it works? I am working with a piece of equipment that does acoustic analysis. It does this by acquiring a time domain signal and performing an FFT on it to extract the amplitude and phase at a single desired frequency, call it f0. Now, I understand the FFT. That part is easy enough. However, what I don't understand is the consequences of all the extra steps I found in this instrument's LabView code. Here's the procedure:

Acquire time domain signal and apply calibration coefficient to every point.
Subtract mean from signal to remove DC component
Mix with a pure tone at f0 by multiplying the complex time domain signal with cos(2πf0) - i sin(2πf0)
Window using exact Blackman
Perform FFT
Take DC component of output, divide by N/2 (where N is the number of samples in the time domain)
Get amplitude and phase at f0 from DC component

Specifically, I'm confused on how step 3 allows me to get the value at f0 by looking at the DC component of the FFT output. I've recreated this in python using numpy.fft.fft and some simulated signals and I get similar results.

At the end of the day, I guess this sounds like a theory/math question.

The answer from Bird has it right.

The described setup is an implementation of a lock-in amplifier.

Lock in amplifiers act as very narrow band filters for detecting the presence of a known signal that is buried in noise.

If you multiply two signals (both simple sine waves) you get a new signal that contains a frequency that is the sum of the original two signals and a frequency that is the difference of the original two signals.  (1000Hz * 200Hz gives an output of 800Hz mixed with 1200Hz.)  This process is used in radio receivers to bring the received radio signal down to a frequency where it is more convenient to do further processing.

Your LabView program makes use of this effect to detect the amplitude and phase of the signal in comparison to a reference signal.

The difference of the input signal and the reference signal is a signal of 0Hz - that's DC.  You get it out of bin 0 in your FFT results.

If the frequency of the input and the reference are not the same, then the difference lands in another bin and is ignored by your detector.

That also makes sense of removing the DC offset from the sample sets - any real DC in the input would cover generated DC from the heterodyning step.

As to why your device is doing that, it is probably so that it can recover a transmitted signal to tell how much of it is coming back to the receiver, and do it pretty much no matter how weak the received signal is or how much noise is present.

Since you say it also takes the phase of the output, it may be using it as part of measuring a distance.

You don't say what the result is used for, or in what conditions you use it so there's not really much more I can say about the "why."


Limits of convolution Say we have a function of time ($S(t)$) of the length $T$, and then a customized impulse response (say $I(t)$) of the length $T+N$. The question is, when $S(t)$ is convolved with $I(t$), what are the possible frequency domain (or alternately time domain) representation of all of the possible resulting functions?

For the case when the convolution is performed circularly, the solution appears immediately clear. If $S(t)$ has zero magnitude components in the frequency domain, they will stay zero no matter how the impulse response is configured. Thus, it's impossible to for example, convert a sinusoidal waveform into any non-sinusoidal waveform using circular convolution.

But the result is a lot less clear for the case when linear convolution is performed. If we limit our horizon to the period of $T$, and completely disregard the tail that the impulse response will generate (for cases where $T>1$, and length of the impulse response $ > 1$ ), what will the answer be then? Are there any limits, or can you transform any function into any other function (provided the tail is thrown away), using linear convolution?

(EDIT: for non-casual impulse responses the "pre-echo" is also thrown away, time shifting is allowed also for such cases).

If the function is zero-valued throughout then so will be the output. But if you exclude this trivial homogeneous case, then your question can be reformulated as:

Given any function $f$ of compact support and any list $l$ of frequency domain zeros, does a function $g$ always exist that equals $f$ within the support of $f$ and has a Fourier transform that has the zeros given in $l$, with other zeros allowed?

Starting at $g_0 = f$, it would be enough to demonstrate a procedure $p$ such that $g_{i+1} = p(g_i, l_i)$ that adds a single given frequency domain zero $l_i$ to function $g_i$ without modifying it within or to the left of the support of $f$ and without altering the existing frequency domain zeros. We choose a causal filter as that procedure. The procedure being a filter guarantees that existing frequency domain zeros are preserved. To not change the function within or to the left of the support of $f$, the filter must have an impulse response that has a Dirac delta function followed by a zero-valued segment of the same length as the support of $f$. A comb filter consisting of two Dirac delta functions, the second one inverted and located after the zero-valued segment, matches this description. The comb filter can be made to have a zero at any wanted frequency by placing the inverted Dirac delta to an integer multiple of the period of that frequency. Choosing a large enough multiplier moves the inverted Dirac delta to outside the zero-valued segment.

If the functions are discretely sampled sequences rather than functions of a real argument, a Dirac delta function that lands between the sample points gives a sampled sinc function that has no compact support. To circumvent this problem, instead of a sinc impulse, a fractional, unity gain delay at any given frequency can be created by setting to certain values the two samples surrounding a multiple of the period of that frequency. The delay and gain at other frequencies will be off but this matters not.

So (unless I missed something), the answer to your question is: No, there are no limitations if the (allowed to be long enough) tail is thrown away.

Note that successive application of the above procedure does not generate the impulse response of your original question but an output (one of many possible) that enables a realizable impulse response.


Detecting which digital number an image represents I photograph the dashboard of my car every day to register the mileage for work.
I use OpenCv to detect the region of the digital displayed mileage, crop to the numbers and match them. The method I use now to determine what numbers are displayed is working for (lets say) 95% of the cases. Finding the numbers is working fine (resize, rotate, correct perspective). Matching the numbers I find still a bit tricky.

The method I use now is determining the center location of each part/stripe that makes up the digit and match the difference of the center to predefined values. To find the center I use blob detection of open cv.

is then converted to a binary image

and blob detection is used to find the center of each blob (marked with red)

The center location of each blob is compared (squared distance) to the center location of a predefined set of locations. And the least square distance is related to what stripe of the digit it is. Then the combination of the stripes found is compared to a dictionary containing the stripe locations for each digit.

The problem I have now is that, although the images are normalized, sometimes more blobs are detected because the the threshold for making the binary image leave out pixels unexpectedly in the center of a blob. And blob detection sees that are separate blobs.

Is there a better approach to determine what digit it is using the particular display of the digits shown above? 'Better' would be: more robust, simpelere, less logic/code required, more universal applicable. I don't want to pull in a whole OCR library for this, if possible.

A simple count of the pixels where the binary color is matching between the image to be recognized and the images of the ten reference digits might be good enough.

For best results, it is advisable to create the reference images using exactly the same method as that for the target image.

If it turns out that binarization doesn't work repeatably enough, you can switch to finer comparison metrics such as Sum of Absolute Differences, Sum of Squared Differences or Normalized Grayscale Correlation.


Understanding Q formats I have some misunderstandings with Q format. Please help me in understanding this.

If I say Q8 [unsigned] then there will be $8$ fractional bits. Is it ? How can I represent a $1$ and $1.5$ in this case?
If I say Q8 [ signed ] then there will be $1$ sign bit and $7$ fractional bits. Here how can I represent a $1$ and $1.5$?
When I say Q1.7 means I have $1$ integer bit and $7$ fractional bits. And the minimum number I can store in Q1.7 is $-1$ and maximum number is $1-1/2^7$. Is this understanding is correct ?

I prefer to use a more explicit notation that consists of the number of integer and fractional bits, like in your last question. I'll take them one at a time:

By Q8, I assume you mean unsigned Q0.8 (no integer bits, 8 fractional bits). In this format, there is no representation for the numbers $1.0$ or $1.5$. The largest value that an 8-bit unsigned integer can hold is $255$, which in this representation would correspond to a fixed-point value of $\frac{255}{256} = 0.9961$.
Similarly to the above, in signed Q0.7 format, there is no representation for $1.0$ or $1.5$. The largest 8-bit signed number is $127$, which would represent $\frac{127}{128} = 0.9922$.
You need to be explicit about whether your number is signed or unsigned. If you have a signed Q1.7 format, that would imply that your values are 9 bits long (1 sign bit, 1 integer bit, 7 fractional bits), which is unlikely. However, assuming that you had such a 9-bit signed number with 1 integer bit and 7 fractional bits, you can represent fixed-point numbers in the range $-2$ to $(2 - 2^{-7})$.

Recall that the conversion between the integer value and the corresponding fixed-point representation is merely a division by two raised to the number of fractional bits. A 9-bit signed integer can hold values in the range $-256$ to $255$, which leads to the fixed-point range I described above.


Are two jointly stationary white noise processes independent? I am currently dealing with a problem concerning beamforming, where two "jointly stationary zero-mean white noise processes" form the input of an adaptive system. One of those processes resembles the actual signal $s[n]$, the other is an interference signal $v[n]$. Is it clear from the stated definition above that they are independent from each other, i.e. $E(XY)=E(X)E(Y)$? Why?

Stationarity (WSS in particular) does not imply independence, even if they are white noises.

That's to say, there exist (at least one) jointly stationary zero mean white random processes which are also dependent. A simple example is this.

Consider a zero mean, wide sense stationary (WSS), white, real random process $x[n]$ with auto-correlation sequence (ACS), $r_{xx}[k] = \sigma_x^2 \delta[k]$, which is passed through a linear time invariant (LTI) filter with real impulse response $h[n] = M \delta[n]$, to produce the output $y[n] = h[n] \star x[n]$ :

$$ y[n] = \sum_k h[k] x[n-k]  = M x[n]$$

Now it can be shown that $y[n]$ will also be zero-mean, WSS white noise with its ACS given by
$$ r_{yy}[k] = h[k] \star h[-k] \star r_{xx}[k] = M^2 \sigma_x^2 \delta[k] $$,

and further, they are jointly WSS and white, as their cross-correlation depends on the lag $k=m-n$;
$$ \mathcal{E}\{x[n] y[m] \}  = r_{yx}[k] = h[k] \star r_{xx}[k] = M \sigma_x^2 \delta[k] $$

However it can be shown that they are not independent, since the necessary condition for their independence does not hold:
$$ \mathcal{E}\{x[n] y[m] \} \neq \mathcal{E}\{x[n] \} \mathcal{E}\{y[m] \} $$

Where the left hand side is the cross-correlation sequence $r_{yx}[k] = M \sigma_x^2 \delta[k] $, as defined above, which is not zero for all $k$, whereas the right hand side is identical to zero for all $n,m$ as both processes are zero-mean WSS.

Hence these two jointly WSS and white random processes $x[n]$ and $y[n]$ cannot be independent and are in fact dependent, in line with intuition as the filtering creates a dependency between them.


Audio Watermarking in MATLAB - How to choose frames to embed in? Hello I am working on audio watermarking in MATLAB and am trying to choose specific frames of the audio signal where the watermark will be embedded. Through research, I have read that the best way to do this is to choose parts/samples of the audio signal that are the loudest, or above a certain threshold. By using wavread in MATLAB, the result is an array of the audio data in type double precision, which represent the sound pressure at a given sample. My question is, since the numbers represent sound pressure, will the loudest portions of the audio signal be the largest numbers? If not, is there another way to achieve what I am trying to do? Thank you!

Loudness is defined as human perception of volume. In general, there is a very good correspondence between the two. Therefore, you can often use one as a proxy for the other.

In this case, you can easily measure volume: the portions of the audio signal with the highest signal intensity (or volume) correspond to the parts with the greatest absolute value (or, if you prefer, square value). With a few assumptions, you can say that the loud sections are also the sections with high volume.


Transfer function and deconvolution Forewords

This question is about methodology references and numerical application. I am posting on Signal Processing because I think this question belong to this place. I am new to the stack, feel free to edit if you feel the Original Post should be improved.

Experiment Setup

Basic experiment setup: inject a dye into a reactor to analyze its dynamic. There are two detector cells (sensitive and selective to dye substance) at in entrance and another at the exit of the reactor. No reaction occurs, the dye is a passive tracer in this experiment.

Two signals are recorded:

$x(t)$, a finite impulse (eg.: we realized a gaussian-like impulse, at the entrance);
and $y(t)$, the response of a system to this impulse (at the exit);

Goal is to compute $g(t)$: the system Retention Time Distribution (RTD).

Signal Processing

For the MCVE sake, I created a synthetic dataset to perform computation on well controlled scenario without noise.

Where $y(t)$ is a convolution such as:

$$y(t) = g(t)*x(t) \xrightarrow{\mathcal{F}} Y(\nu) = G(\nu)\cdot X(\nu)$$

Now, it might resume to get the inverse transform of the transfer function:

$$ G(\nu) = \frac{Y(\nu)}{X(\nu)} $$

But numerically, even for this smooth example, it does not work well. Mainly because $X(\nu)$ has zeros and small numbers that make the operation unstable.

Alternative

There is a processing (adapted from this example) that seems to solve my problem:

$$ G'(\nu) = \frac{X^*(\nu)\cdot Y(\nu)}{X^*(\nu)\cdot X(\nu) + \lambda^2} $$

Where $\lambda$ is related to Signal-to-Noise Ratio:

$$ \lambda^2 = \frac{1}{\mathrm{SNR}} $$

I feel comfortable with this solution, because the operation tends to the former when $\lambda$ goes to zero:

$$\lim\limits_{\lambda \rightarrow 0} G'(\nu) = G(\nu)$$

Then I can estimate RTD by taking the inverse transform of $G'(\nu)$:

$$\hat{g}(t) = \mathcal{F}^{-1}[G'(\nu)]$$

Numerical Example

For the MCVE dataset, it renders as follow:

Derived RTD (after inverse transformation) quite fits the original function:

And error on this method is acceptable for the MCVE:

Even if there are artifacts at the very beginning and end of the recovered signal.

Observations

Setting $\lambda$ to higher values change the error behavior, but the RTD looks fine:

Setting  $\lambda$ to smaller value may affect the shape of the RTD:

Questions

Reading on the internet I had difficulties to find references where the transfer function is the unknown. Most of the articles speaks about deconvolution and estimation of $\hat{x}(t)$ knowing $y(t)$ and $g(t)$.

I also found the Wikipedia page of Wiener Deconvolution which shows a formulae that is very similar to $G'(\nu)$ in my OP. But I could not relate them together. I miss something to get the full understanding on why this methodology is correct, if it is, when performed as I did.

My questions are:

How can I formally justify the expression of $G'(\nu)$? Can it be shown there is a relation with Wiener Filter?
Is there a good reference to apprehend this class of problem?
How is called the oscillating effect at the ends of the signal?

Your help is greatly appreciated.

Code

Below the Python code of the MCVE:

import numpy as np
from scipy import stats, signal, fftpack

# Setup:
fs = 10   # Sampling Frequency [Hz]
dt = 600  # Experiment runtime [s]
l = 1e-4  # 1/SNR in term of amplitude [-]
t = np.linspace(0, dt, dt*fs+1) # Time [s]

# Signals [AU]:
xt = stats.norm(loc=20, scale=5).pdf(t)
gt = stats.lognorm(0.45, loc=-20, scale=100).pdf(t)
yt = signal.convolve(xt, gt, mode='full')
yt = yt[:t.size]/np.sum(yt[:t.size]/fs) # Normalization

# Fourier Transform:
Xs = fftpack.fft(xt)
Ys = fftpack.fft(yt)
fq = fs*fftpack.fftfreq(t.size)
Gs = (np.conj(Xs)*Ys)/(Xs*np.conj(Xs) + l**2) # Transfer function

# Inverse Transform:
gti = fftpack.ifft(Gs)
np.allclose(np.imag(gti), 0) # True
gti = np.real(gti)
gti = gti/(gti.sum()/fs) # Normalization

# Error:
ea = gti-gt
er = ea/gt

How can I justify the expression of G′(ν)?

That's easy enough. Denominator is the sum of the signal energy $|X(\omega)|^2$ and the noise energy $\lambda ^2$. If the signal energy is significantly larger, then the whole expression simplifies to $G(\omega)$. If the noise is larger, we can't do a anything useful with the information and the $1/ \lambda ^2$ term will drive the whole expression towards zero and some continuous way.

  Is there a good reference to apprehend this class of problem?

Deconvolution does work fine here. Convolution is commutative, so it doesn't really matter if you are trying to calculate the impulse response or the input signal (at least in your case).

There are a few other methods that can also work well which are typically based on a least square error optimization in either the time or the frequency domain.

Any of these methods can be greatly improved by using any a-priori knowledge of the system. For example, if you that it's a relatively well behaved physical system with a some known degrees of freedom, you can simply do a least square error fit trying to determine a small number of poles and zeros.

Each application has it's own specific requirement and the best method depends a lot on what you have and what exactly you need.

  How is called the oscillating effect at the ends of the signal?

Not sure. I think your problem here may be that your noise estimate is not a function of frequency and that somehow a single or a few discrete frequencies get overal amplified in the error signal.


What is the mixed numerology in OFDM In order to meet the user requirements in the next
generation radio access technology, the ongoing standardization within 3rd  Generation  Partnership  Project  (3GPP)  offers flexibility in terms of employing mixed numerology.

I have read many papers about it, such that: This link.

I think it's a new technique but it's difficult for me to understand it well. It's like to use different carrier spacing instead of one fixed of 15 KHz in LTE, but that's not sure.

Could you please explain what is that mean in an easy way, and what's the benefit of it?

In general, it means several waveforms will coexist in 5G ecosystem.

In the context of LTE-based 5G (up to now), several subcarrier spacings (hence CP) coexist. The first and foremost motivation for this is to support not only mobile broadband but also other use-cases such as mMTC and URLLC.

Some of benefits of mixed numerology are flexible scheduling and better spectral efficiency.

Anyway, I do agree to other comments that the question is a little too broad.


What is the state-of-the-art method for dealing with illumination issues in background subtraction? I have been struggling to find a well-cited approach to this.  Essentially, I am experimenting with different background subtraction algorithms, but none seem to perform well when there are significant lighting changes (say, from a window in the scene getting significant glare at specific times).

The question is simple (relatively): what are the state-of-the-art methods for background subtraction with sudden illumination changes?  My scenario is surveillance cameras.

The standard approach is to use a mixture of Gaussians to model the various laws that command the intensity of a pixel. In this approach, a pixel can be assigned several normal laws (with different mean and variance), each of them will model its value under different conditions.
Of course, only 1 Gaussian needs to be active at a time.

While the maximum number of Gaussians in the mixture is a (fixed) parameter, the parameters of each Gaussian are learned online.
You can find an example implementation in the OpenCV library.

A few years ago, someone (sorry, i can't find the name anymore) proposed an interesting alternative approach: using variations of the image gradient orientation instead of variations of the light intensities.
The gradient orientation has the advantage of being contrast invariant, making it more robust to illumination changes.
Intuitively, it works because the gradient orientation is tied to the shapes in the image, not to their color or luminosity.


IFFT of a complex, conjugate-symmetric transfer function - Value at Nyquist frequency Using the MATLAB 'symmetric' flag to compute the IFFT of a complex-valued, length-N (with N even), transfer function, the imaginary part of the complex value at the index ((N/2)+1) corresponding to the Nyquist frequency appears to be ignored. Why is this done instead of taking the magnitude of the complex value at the Nyquist frequency? How important are the values at the Nyquist frequency and DC for audio applications?

If you have an exactly conjugate-symmetric frequency domain vector, the values  at DC and at Nyquist must be real-valued. The 'symmetric' flag of Matlab's ifft command makes sure that the result of the inverse FFT is real-valued, implying that the values at DC and Nyquist must also be real-valued.

Let $X[k]$ be the length $N$ DFT of a real-valued sequence $x[n]$:

$$X[k]=\sum_{n=0}^{N-1}x[k]e^{-j2\pi nk/N}\tag{1}$$

The value at DC is

$$X[0]=\sum_{n=0}^{N-1}x[k]\tag{2}$$

which is clearly real-valued if $x[k]$ is real-valued. Similarly, assuming that $N$ is even, the value at Nyquist is

$$X[N/2]=\sum_{n=0}^{N-1}x[k]e^{-j\pi n}=\sum_{n=0}^{N-1}x[k](-1)^n\tag{3}$$

which is also real-valued for real-valued $x[n]$.

Conjugate symmetry means

$$X[k]=X^*[N-k]\tag{4}$$

which for $k=N/2$ (i.e., at Nyquist) gives

$$X[N/2]=X^*[N/2]\tag{5}$$

Of course, $(5)$ just shows in another way that the bin at Nyquist is real-valued for real-valued $x[k]$.


Finding similar pictures in a database I have large database of pictures, about a few million.

I need an algorithm to find pictures which are similar to those in the database, up to a scale factor or a crop and shift (I have only these types of variations).

I know about histogram matching and perceptual hash, but is there a faster algorithm that exist for these purposes?

For this you can use local features like SIFT, SURF, FREAK, etc.


How does BRIEF produce consistent keypoints with random sampling? I'm reading through BRIEF and one point isn't clear to me. If we are randomly sampling the pixel-vs-pixel tests, how is it that we obtain a reproducible descriptor?
EDIT
This is a duplicate. Find question with answer here

It looks to me as if they came up with a random pixel patch vs pixel patch set of points and they used this to train their system.  The system then recognizes different transformations by looking at the database of transformations and seeing how the transformations would affect the match.  It then picks the transformation with the best score (lowest Hamming distance).  The random pixel locations are only random once during the system design.


Why do we need a Homomorphic Filter? Please, explain in layman's terms,

Why do we need a Homomorphic Filter?
In what kind of situations should we use Homomorphic filtering?
What kind of images need Homomorphic filtering?

Please, tell what is not enough "layman".

Because standard linear smoothing (like averaging) tends to filter noise but reduces the contrast. Image composition is not often linear: cast shadows, occlusions. So linear methods tend to be suboptimal.
When you want to preserve sharp contrasts, for instance with multiplicative noise (whose intensity varies with image intensity), when you have a slow variation of the global illumination across the image (like progressive shadowing effect), when you want to both denoise and enhance
Images taken under varying illumination conditions, with noise varying with intensities, for instance related to temperature (infrared, visible of hot objects).


is there a use for a changing kernel in image convolution I understand the principle of image convolution. Maybe. I think it is this:

You use a small matrix with numbers in it,

then take a sample of an image with as many pixels as you have numbers in your matrix,

then do a dot product, then add up all the terms and

then use the result as the new value for the pixel located in the center of your image sample.

Now, my question is this:

Does it make sense to think about an image kernel for convolution whose values actually depend on the image? For example, consider a 'sharpen' kernel:

0 1 0
1 4 1
0 1 0

I want to change the value in the center to be a number between 1 and 5 depending on the brightness of the pixel which is two places to the left.

The reason I ask this is because I've noticed that when I look at an image, I actually look at 'neighborhoods'. not just individual pixels. While bright spots might catch my attention. I recognize parts of an image before it 'clicks'. So if I am using some parts of an image to inform me about the 'neighborhood' that part lives in, I am wondering if this sort of 'dynamic kernel' could produce a matrix that informs us more meaningfully about the image. Obviously I am naive and know very little about image processing but a lot of papers on image convolution get technical very fast and then I'm lost.

Yes it's true; there are image filter kernels whose coefficients may change with location of the processing and further based on image data. They're known as shift-varying filters. And data dependency also makes them nonlinear.

One prime application of shift-varying nonlinear filtering is on the edge adaptive, noise reduction where the kernel coefficients are determined based on a statistical analysis of a local group of pixels; according to whether it's an edge dominant or noise dominant block. There are different applications too.


ADC performance simulation: How to calculate SINAD from FFT? While working on this problem, I started to have doubts that my initial definition of

$$SINAD = 10 \log_{10} \left( \frac{p_f} {\sum_i{(p_i)} - p_0 - p_f} \right)dB$$

is correct. In this equation, $p_x$ is the power the FFT bin at frequency $x$, $p_f$ is the power of the frequency bin containing the signal frequency $f$ and $p_0$ is the DC component. The sum over $i$ accumulates all frequency components, before removing the DC component $p_0$ and the the signal frequency $p_f$.

More specifically, I am unsure about the $\sum_i(p_i)$ part, which I interpreted from the Wikipedia description

  The ratio of (a) the power of original modulating audio signal, i.e., from a modulated radio frequency carrier to (b) the residual audio power, i.e., noise-plus-distortion powers remaining after the original modulating audio signal is removed. With this definition, it is possible to have a SINAD level less than one.

Comparing with the equation, the "original modulating audio signal" is at frequency $f$, which is accounted for in the $p_f$ term from the FFT. The $p_0$ term I got from the following paper, which says to remove the DC component:

In the "Understand SINAD, ENOB, SNR, THD, THD + N, and SFDR" document, it says

  Signal-to-Noise-and Distortion (SINAD, or S/(N + D) is the ratio of the rms signal amplitude to the mean value of the root-sum-square (rss) of all other spectral components, including harmonics, but excluding dc

Looking at these definitions, I can think of another possible definitions of SINAD, namely

$$SINAD = 10 \log_{10} \left( \frac{p_f} { \sqrt{ \sum_i{(p_i^2)} } - p_0 - p_f} \right)dB$$

which uses the RSS (root-sum-square) of the noise and distortion bins of the FFT result. But then, what exactly is meant by "mean value" in that document?

If SINAD can be determined from only expected value and variance then it is possible to determine how SINAD transforms. Variance $\sigma^2$ is preserved whereas expected value $\mu$ grows as $ \sqrt{N} $ where N is the sampling set size. Noise + distorsion is assumed to have variance $\sigma^2$.

Thus the SINAD value would then be determined to become
$$
{\mathrm {SINAD}}={\frac {P_{{\mathrm {signal}}}+P_{{\mathrm {noise}}}+P_{{\mathrm {distortion}}}}{P_{{\mathrm {noise}}}+P_{{\mathrm {distortion}}}}} =
{\frac {N |\mu|^2+\sigma^2}{\sigma^2}} = N \cdot SNR + 1
$$

I can explain this in more detail if needed.


Problem with Covariance matrix using diagonal loading involved in calculation of eigenvalues I have the following problem : I'm calculating the sample covariance matrix in the frequency domain ( $y_{k}$ is the FFT of a time domain $k_{th}$ symbol  vector signal , basically a simulated received signal) as follows:

$$\mathbf{R}=\frac{1}{N_{f}}\sum_{k=0}^{N_{f}-1}y_{k}y_{k}^{H}$$

Well , the next step on my algorithm is to solve a optimization process,  as an essential part of it I need to compute the eigenvalues by doing SVD ,method of powers, etc... in MATLAB of the following equation:

$$(\mathbf{R}^{-1}\mathbf{A})$$

To not get into too many details because I believe my issue comes from a numerical problem and many insights of the actual algortihm / context aren't needed .Let's assume, $A$ is simply a predefined matrix that I compute.

The REAL ISSUE appears now because $\mathbf R$ is ill-conditioned as MATLAB tells me so. So the inverse procedure seems to be failing and the eigenvalues I'm obtaining are incredibly small due to this issue (In fact I only get 1 eigenvalue different than zero). The dimensions of $\mathbf R$ are typically large ( since they are compressed depends on the actual compression ratio I'm using but let's say $32\times 32$ for example).

One approach to solve this problem I found is to use diagonal loading as:

$$(\left(\mathbf R+\sigma\mathbf I\right)^{-1}\mathbf A)\quad\text{with}\quad\sigma > 0$$

This seems to be solving the problem, the eigenvalues are now scaled due to this background "noise". My question is how can I obtain the truly well scaled eigenvalues because , later on,  in my algorithm these eigenvalues will serve as weights since I'm considering them as an actual power estimate.

Note : I've been playing with the cond() function in MATLAB for $\sigma= 0.05$ ,cond(R+sigma*I)= 2 , which is not bad I believe.

Feel free to ask more questions about the problem. But I think my question relates to a purely numerical issue involving eigenvalues and covariances matrices.

You can write

$$ R=YY^H $$
where $Y$ is a matrix of size $N\times N_f$ and $N$ is the dimension of $y_k$. $Y$ contains all the measured $y_k$ as its columns.

Then, the rank of $R$ is upper bounded by $N_f$. In particular, if $N_f<N$, $R$ will always be a singular matrix. So, if you have too few measurements, you will likely run into the problem you described.

For your second part, how do the Eigenvalues of $R$ and $R+\sigma I$ relate:
Consider the eigenvalue decomposition of $R$ by
$$ R=PVP^{-1}$$
where $P$ contains the eigenvectors and V is the diagonal matrix with the eigenvalues (they can be zero). Then you have

$$ R+\sigma I=P V P^{-1} + \sigma I = P(V+\sigma I)P^{-1} $$

This means that the eigenvalues are all increased by $\sigma$.


Sound recording anonymization I want to study stress based on sound recordings. However, there are many concerns about keeping the identity of the recorded persons anonymous.
I want to keep intensity and spectral contents intact but I am not interested in the message or the identity of speakers.
I thought about time reversing the recordings, but this could be easily undone by someone.
Do you have any suggestions for this purpose?

I don't really think there's one:
You seem to be in the process of figuring out what about a piece of audio signifies stress.
As long as you don't know that, you can't erase any identifying information from it without standing a high probability of also erasing your "stress" information.


Blackman-Tukey periodogram? I have been reading about the non-parametric estimation of the spectral density, and I'm particularly interested in the Blackman-Tukey method. If we take the rectangular window we will have the truncated periodogram, computed on a reduced covariance sequence length $M$, where $M<N$:
$$\hat{P}(\omega)=\sum\limits_{k=0}^{M-1}\hat{r}(k)e^{-jk\omega}$$
where $\hat{r}(k)$ is the estimate of the covariance of the sample and $\omega$ is the frequency of interest. What I'm having trouble is: first, how do we choose our new sample of length $M$? could someone please explain how it does actually work without being too technical, like this explanation on how does the Bartlett method works

The resolution of the Blackman-Tukey periodogram estimator is $\mathcal{O}(1/M)$, whereas its variance is on the $\mathcal{O}(M/N)$. This compromise between resolution and variance is to be considered when choosing the window’s length.

I have paraphrased Peter Stoica's beautiful explanation from here (check Section 2.5)


Why samples per symbols $\geq 2$ in the GNURadio Constellation Modulator block? Let's say I want to emulate the design of a modulator I am building "in real life".  I'd like the components, input data rate, and output bandwidth to be the same as in the real system.  If I use the GNURadio Constellation Modulator block for this, it imposes a lower limit on samples/second.  This results in a different bandwidth for the GNURadio graph, than the real system.
Why does this block impose this restriction?  And, is there a canonical (or preferred) way around it, when emulating real systems?

The GNURadio Constellation modulator enforces a minimum of 2 Samples Per Symbol mostly due to simplicity and some practical reasons. Theoretically, if you're just transmitting PSK without any pulse shaping, you could do 1 sample per symbol (note these are complex symbols) . But typically we like filtered PSK, and the block does this with the RRC filter. However - a key point is that when filtering the signal, you have bandwidth expansion. So if you use a 35% excess bandwidth, your sample rate must be at least 1.35x larger , or another way to say it is that you need at least 1.35 samples per symbol to satisfy nyquist. Similarly for 50% bandwidth you'd need 1.5 samples/symbol. The easiest way to deal with this is to just upsample by 2 (i.e. oversample) and filter with an RRC designed for 2 samples per symbol and then you can support the full range of excess BW values that the RRC can do.
There is one other key idea here though - if you were to try doing a rational upsampling by 1.35x for a 35% excess BW RRC, your resulting signal would be fully occupying your bandwidth (i.e. you have no oversampling). This becomes a problem in multiple areas - for one, when you upsample by 2 later on, the image that must be filtered will show up right next to the one you want to keep, making your digital filtering requirements harder. Additionally, if you instead plan to send this out via a D/A, the image in your 2nd nyquist zone will be sitting right next to your desired image in zone 1, and again filtering with your analog reconstruction filter will be very difficult (e.g. a complex filter with steep transition band).
Also take a look at this very similar question asked by myself where Dan Boschen gave a great answer: Sample rates, Samples per Symbol, and Digital Pulse Shaping


Newbie looking for an interesting introductory book to telecommunication and signals/data processing I am currently in my first year of Bsc Applied Physics. Over the past two years, I have realized I'm interested in majoring in telecommunication after I finish my bachelor. Also, for one of my classes I need to give a presentation; I chose my topic to be about Data Transmission. While working on my abstract for my project, I came across things I never even thought of! This, however, sparked my interest in this field even more!

So is there any interesting introductory book to signals and telecommunication that you would suggest and that is interesting and accessible for a bachelor student?

Thank you very much!

Yours,
WBk

"Telecommunications" is too broad a subject: you could spend your entire life learning about just a small part of the field. However, you mention data transmission specifically; I'll recommend these books:

"Digital Signal Transmission" by C. C. Bissell and D. A. Chapman, Cambridge. This is an undergraduate level book that has very lucid coverage of both the physical layer and larger telecom networks.

"Principles of Digital Communicatios", by B. Rimoldi, Cambridge. This is an advanced undergrad textbook that focuses on the math and theory behind digital communications. A great book if you like that kind of thing.

"Data and Computer Communications", by W. Stallings, Pearson. I'm not the biggest fan of Stallings, but his books are very accessible and cover all the bases. In fact, it might be a good idea to start with this book and then continue with more rigorous books like the first two I listed.


Correlation of discrete time periodic sequences So I was reading this textbook -  Digital Signal Processing by Proakis and Manolakis. The authors define the cross-correlation operation for 2 energy signals $x(n)$ and $y(n)$ in section 2.6.1 as:

$$ r_{xy}(l)=\sum_{n=-\infty}^{\infty}x(n)y(n-l) $$

In section 2.6.3, they define the cross-correlation operation of two power signals as:

$$ r_{xy}(l)=\lim_{M\to\infty}\frac{1}{2M+1}\sum_{n=-M}^{M}x(n)y(n-l) $$

My question is, why do we need to have a separate definition for power signals? Are these two definitions the same? If they aren't same, what is the intuition behind defining them like that?

Energy signals (such as an exponentially decaying pulse) are bounded in total energy as time goes to infinity in which case you can use the first definition. Power signals (such as a sine wave) have finite energy over finite time intervals but infinite energy for all time (power is energy/time) so the first equation will not converge in that case.


Decimator effect on wide sense stationary input I've seen that the output of a decimator when a WSS process is passed through remains WSS. I am not able to immediately see why this is. What is a good explanation of why the signal maintains stationarity?

The decimator by integer $M$ can be shown to be the following block:

$$ x[n] \longrightarrow \boxed{ \downarrow M } \longrightarrow y[n] = x[Mn] $$

Assuming that the input $x[n]$ is WSS it has the ACS as
$$ r_x[k] = E\{ x[n] x^*[n+k] \} $$

Then the output autocorrelation can be defined as:
$$ r_y[n,n+k] = E\{ y[n] y^*[n+k] \} = E\{ x[Mn] x^*[Mn+Mk] \} = r_x[Mk]$$

as can be seen, the output auto-correlation $r_y[n,n+k] = r_y[k] = r_x[Mk]$ also depends on the lag $k$ and therefore we can conclude that $y[n]$ is also WSS.

(of course the mean of $y[n]$ is also independent of time $n$.)


Complex signal transform I must strengthen my knowledge in the mathematical properties of signals, and would like to know the following thing:

if I have a complex signal, written in magnitude and phase, it's like if I have 1 signal that holds 2 signal himself? I have this question because I have to do the inverse transformation of a signal in the frequency domain , that is complex and I have the magnitude and phase plot. Is it right that my phase, when I do the inverse- transform retards my signal ( the inverse transform of the magnitude) in the time domain?
Thank you so much.. I have not on focus how the phase works yet

Thank you for the answer, I know what you have written... my question is:
If I have one complex signal in the magnitude and phase form, do I have to consider magnitude and phase as two different signals? Look

The signal is written in  magnitude with a complex exponential with the phase... when I do the inverse transform the exponential is in the time domain a delay ... and this it feels strange that one signal bring with him a delay ... it's like if phase is another signal , that retards the inverse transform of the magnitude.
I hope I have explained my question in the right way

A complex signal can be written as $s(t)=a(t)+jb(t)$, which indeed corresponds to two real signals ($a(t)$ and $b(t)$) "embedded" in the single signal $s(t)$. What makes this interesting is that $a(t)$ and $b(t)$ can be easily recovered from $s(t)$. This is true whether you write the complex signal in rectangular or polar coordinates.

Now, regarding how to interpret the phase of a signal. Let's say that the Fourier transform of $s(t)$ is $S(f)=|S(f)|e^{j\phi(f)}$. Imagine that you have an infinitely narrow band-pass filter with center frequency $f_0$. Then, you connect the signal $s(t)$ to the filter's input and see the output in an oscilloscope. You'll see a sinusoidal wave, with frequency $f_0$, amplitude $|S(f_0)|$ and phase $\phi(f_0)$. Here, "phase" indicates the difference, in radians, between a zero-phase cosine $\cos(2\pi f_0t)$ and the sinusoid you see in the scope.

In that sense, the phase $\phi(f_0)$ can be thought of as a "delay"; the sinusoid in the scope is just a cosine with a delay, $\cos(2\pi f_0t+\phi(f_0))$.

As an application, consider a signal that goes through an LTI system. The system doesn't touch $|S(f)|$ but it modifies the phase $\phi(t)$ -- that is, it changes the phase of each individual sinusoid. If this results in the same delay for each and every sinusoid, then the system output will be a delayed version of $s(t)$. If some of the sinusoids suffer different delays, then the signal at the system's output will be a distorted version of $s(t)$.


How to smooth (or interpolate) phase of FFT and reduce data points I am writing some code for audio analysis, and have currently got two signals with FFT performed on them.
I get the phase of my complex array by using:

ang = Math.atan(complex.im/complex.re)*180.0/Math.PI;

I then 'format' for wrapped phase:

if(complex.re < 0.0 && complex.im == 0.0){ang = 180.0;}
            else if(complex.re < 0.0 && complex.im == -0.0){
                ang = -180.0;}
            else if(complex.re < 0.0 && complex.im > 0.0){
                ang += 180.0;}
            else if(complex.re < 0.0 && complex.im < 0.0){
                ang += -180.0;}
            phase.add(ang);
        }

As I have converted this to degrees already, I then subtract channel 1 phase from channel 2 phase to give me the phase difference:

double phaseDiff = (((double)chan1phase.get(i)) - ((double)chan2phase.get(i)));

Finally, I add a method for keeping wrapped phase whilst subtracting the two phases:

if (phaseDiff < -180){
            phaseDiff += 360;
        }
        if (phaseDiff > 180){
            phaseDiff -= 360;
        }

This works really well currently,  however there are too many data points for my application to respond quickly. I am working with 16k or 32k FFT size, and that's a lot of data to plot and redraw quickly when I manipulate my data. I also occasionally get one or two data points throughout the whole plot which are noisy or differ from the expected path of travel, which I would like to get rid off.

I have already used a smoothing function to display my magnitude vs frequency, and this responds really quickly.  However I am not confused on how best to smooth phase.  I need to reduce the data points to approximately 1000 points over the audible frequency range (20Hz - 20kHz).

However I cannot just use the same method as I used for magnitude because of the wrapped phase.  Would I need to average the data points within an octave (or preset band) prior to wrapping the phase?

Or am I better disregarding some points and interpolating the data?

Any guidance would be appreciated. Apologies if there is too much actual code here, I realise this is a DSP site and not programming site, but thought it useful to see exactly what I am doing (or where I have gone wrong!).

To downsample an FFT result for both magnitude and phase, it may work best to do an FFTShift before the FFT, then downsample the real component vector and then the imaginary component vector separately, before taking the atan2() of them to estimate downsampled phase.


Finding periodicity in discrete events I have a program that detects events in a large amount of measurement data. When it detects an event, it writes a timestamp. I have thousands of event timestamps. What I wish to do is detect if there is periodicity in the timestamps I have.

Pictures may aid my explanation: If I have a bunch of events on a timeline, as in the figure below, the events all seem to be random but there may be some kind of periodic component to the events.

What I wish to do is detect if any of the seemingly random events are following a strict repeating interval. An illustration is given below, where we see that in the seemingly random events above there are some data points that are repeating with a fixed frequency.

I am not certain what kind of method to apply. I have looked into power spectral density, fourier transform and ARIMA. I applied fourier transform (as in the answer here) but was not able to make anything meaningful out of it.

Properties of the applied method should possibly include:

A quantitative measure of how strict or fixed the intervals are, or how
certain we can be that we have detected a fixed cycle
The ability to detect periodicity on different timescales (e.g. events occurring multiple times within the same hour or multiple times during a week with a fixed pattern)

Anyone has suggestions how this could be done?

See "ePeriodicity: Mining Event Periodicity from Incomplete Observations" by Zhenhui Li, Jingjing Wang, and Jiawei Han (2013); preprint at https://faculty.ist.psu.edu/jessieli/Publications/tkde14.pdf


Analytical expression for step response of digital bessel filter If I apply a digital Bessel filter to a perfect step function, I get something that looks like the following:

The green line is the input step response data sampled at 500kHz, the red line is obtained using the scipy.signal.lfitler routine with a 8-pole Bessel low-pass filter at 10kHz, and the blue line is the result of the same filter but using the scipy.signal.filtfilt.

Given the signal (either red or blue) is there an analytical form for the response that would allow me to extract the cutoff frequency using a nonlinear fit? Both functions look like a sigmoidal of some kind - is there an analytical formula for it?

The form

seems to fit both cases quite closely, but not perfectly, and it's not clear what relation the time constant bears to the cutoff frequency.

If I fit that form to the filtfilt data, I get a linear relationship between the time constant for the fit, and the time constant for the filter (1/fc), with a slope of about 1/35. I imagine that slope will be a function of the filter order. Can anyone suggest a better analytical form for the fit function?

Yes, there is an analytical response, though it might not look like you would like to have it:

Your Bessel Filter can be written in Z-Domain as

$$H(z)=K\frac{\prod_{i=0}^7(z-z_i)}{\prod_{i=0}^7(z-p_i)}$$

where $z_i$ are the zeros and $p_i$ are the poles of the z-Transform. Now, to get the step response, you "simply" need to calculate

$$h[n]=\mathcal{Z}^{-1}\left(\frac{z}{z-1}H(z)\right).$$

i.e. you multiply the Z-Transform of the filter with the Z-Transform of a unit-step (which is $z/(z-1)$) and calculate the inverse Z-Transform. If you have access to Mathematica etc. it can do it for you. If not, you would need to perform partial fraction decomposition of the expression and then perform inverse Z-Transform of the several partial fractions.


Could not construct original matrix using SVD I need to

read a .wav file
make a 2D matrix from it using the STFT
generate its singular value decomposition
do some watermarking

but I can't build the original matrix again via matrix multiplication:

 x = audioread('0.wav');
 stft = spectrogram(x);
 [U,S,V] = svd(stft);
 stftb=U*S*V.';

The problem is that STFT and STFTb are not the same.

Don't have enough reputation to comment but you need to use the conjugate transpose in your formula for the result to be correct. So try stftb=U*S*V'; in the last line of code. Note that I removed the . which makes a difference since the matrices you are working with are complex.


Zoom an image with respect to a given image I posted this question in stack overflow, but I think it is more appropriate here:

I have a reference image - 'image_01.png' of size 200x200 and resolution 72 px/in.

I have another image 'image_02.png' (with the same size and resolution) with several levels of zoom, of which, one is close to 'image_01.png', in the sense that the objects in that image are of "similar" size. They are NOT the same images, btw.

Here is an example with three levels of zoom, and clearly zoom_02 is "closer" to image_01 than zoom_01 and zoom_02. My question is, what mathematical metric (ratio, or any measure) do we choose to determine what zoom is appropriate to ensure such similarity? Let's say I even have a screen ruler to measure the distance in the images, in pixels.

 zoom_03  zoom_01

 zoom_02  image_01

Since you want the image size and resolution to remain constant, the term optical magnification seems apt there.

You could then say

  Please use 20x optical magnification to identify said feature.

Please be careful that the term magnification is different from zoom.

Two different cameras with the same 5x optical zoom might produce images with different optical magnification depending on their optics.

Zoom refers to the ratio of the max to min focal length of your camera, whereas,
Magnification refers to the ratio of projection size to the real subject size.

Given a particular magnification, you could use your camera focal lengths and apertures to come up with at zoom metric for that particular camera to produce the desired results.

See this for a better explanation.


Design a filter which passes all frequencies except $\omega=\pm\frac{\pi}{2}$ and plot its pole-zero diagram Also draw its normalized frequency response.

What is the ROC?

This has to be done in z-plane so there must be two poles at $+i$ and $-i$ since they cannot be included in region of convergence. Is my assumption correct?

And what is meant by normalized frequency?

The kind of filter you are looking for is a notch filter. Using filter design toolbox in Matlab you can get it as I'm showing you in the following picture:

  This has to be done in z-plane so there must be two poles at +i and −i since they cannot be included in region of convergence. Is my assumption correct?

Yes! This is correct.

Assuming you want to create a causal filter (a right sided sequence with non-zero samples only for n > 0) then your Region of Convergence ROC is :

If you want to understand this better I would recommend you take a look at example 3.1 of Alan V. Oppenheim Discrete-Time Signal Processing 3rd edition.

Normalized frequency response means that instead of dealing with a sampling frequency you make a design independent of it in which you only care about digital frequency. For instance if I say that $w = \pi/2$ then if sampling frequency is 48 kHz (this is $2*\pi$), the maximum frequency that you may get is $\pi$ or 24 kHz and the notch will be at 12 KHz which is $\pi/2$. But on that very same filter if sampling frequency is later 96 kHz then the notch will be at 24 kHz.

The way in which you can do this filter in Matlab is this:

If you select a broader bandwidth (for instance 0.1 instead of 0.001) you will notice how the zero and the pole maintaining the same angle start to separate from each other in the diagram

with sampling frequency 48 kHz:

with sampling frequency 96 kHz:


Matlab : Help in getting the correct constellation diagram for QAM modulation I would like to see the constellation diagram after equalization as shown in the website Correct constellation for QAM. But I am getting a different plot. Can somebody please tell me where I am going wrong?

You are getting the correct constellation. You need to remember that you added noise to your symbols in the channel, So, you won't get a perfect constellation as you have in the first figure. Try increasing SNR (EbNo) and you will see better the constellation.

This is what I got:


A question regarding z transform and its magnitude response My teacher of signals and systems gave us a review problem as following:

given a DT rightsided LTI system with transfer function

$$\frac{1-a^*z}{z-a}, \left | a \right |<1
$$
show that the system's magnitude response is in unity for all frequencies using graphical method.

Can anyone give me any suggestion of solving this problem? I've been going back and forth between the textbook and course video for hours, but still very clueless.
I am particularly confused about the indication of using graphical method, because it seems no matter how I manipulate the fraction, there's always vector multiplication instead of simple addition.

The frequency response of a given filter is calculated evaluating $H(z)$ in the unit circle, i.e. $z=e^{j\omega}$.

Hint 1:

Note that you can express your transfer functions as follows:
$$H(e^{j\omega})=\frac{1-a^*e^{j\omega}}{e^{j\omega}-a}=\frac{1}{e^{j\omega}}\cdot\frac{1-a^*e^{j\omega}}{1-ae^{-j\omega}}$$

Hint 2:

There is a property that states that given any complex number $z$, the following equality holds:

$$|z|=|z^{*}|$$

Can you solve it now?


Is it possible to generate Hamming code with minimum distance $d$ which can correct more than 1 bit error? In literature, all I am seeing is Hamming code with minimum Hamming distance $3$. Theoretically a Hamming code with minimum distance d can detect $d-1$ errors and can correct $(d-1)/2$ error. So minimum distance of $3$ can detect $2$ errors and correct $1$ error.

So here is my question: Is it possible to generate a code with minimum Hamming distance of (let's say) $5$ which can correct $2$-bit errors in each block of coded bits?

If yes, could anyone please give me suggestion about that?

If not, then please give me an explanation.

You cannot detect 2-bit errors (reliably) with a Hamming code distance of 3 if you are correcting 1-bit errors.

At any rate, you can design Hamming codes correcting two errors.  Let's assume 32bit of data and n Hamming code bits.  The states we then need to be able to encode is
0 bit error = 1 state
1 bit error = 32+n states
2 bit error = (32+n)(31+n)/2 states
For n=10, this amounts to 1+32+10+21*41 = 904 states that can be encoded with 10 bits of redundancy.  Depending on how complex you want your correction machinery to get, you might want to use more bits in order to make multiple error correction less complex.


How to add the noise covariance matrix of my measurements to tmy 1D kalman filter? So I am new to kalman filters. This is my first experience with it.

I have some measurements but sometimes due to external factors those measurements become totally wrong (hence the peak when you run my matlab code). If I understand the principles of kalman filtering correctly, I can compensate this incorrect measurements by setting the measurement covariance matrix R as a variable, where R represents the accuracy of the measurement devide. But so far I am not sure about where I should add or how I should include this measurement covariance matrix in my measurements.

Could anybody help me out?

I have tried as follows: https://pastebin.com/KQEx5T10 But am totally not sure about my results.To me it looks like both of the outputs of my kalman filter are too far away from the actual measurements once the peak happened... I would have expected the output curves to be between the "toothshaped" curve and not above it.

If you are so new to Kalman filtering, the best first suggestion should be the most readable book on the topic: Fundamentals of Kalman Filtering, from Zarchan.

In the mean time, the measurement noise covariance matrix $R$ is adjusted according to the a priori known (either empirically or theoretically) characteristics of the measurement device (the sensors). When these characteristics change during the operation time of the filter, then you would expectedly reflect that in an update of the measurement noise matrix $R$ according to the new settings.

You will include this measurement noise covariance matrix into the Kalman update equations. Which previously using fixed $R$ will now be using variable $R$ matrix, as a matter of fact.


Advice on autocorrelation via FFT I have been creating an application in Java that transforms an audio signal and writes it to a midi file.

At first I tried using autocorrelation to find the fundamental frequency. However, I have been having difficulties implementing it. I resorted to finding the index of the peaks of the FFT output, finding their consecutive differences and deriving the fundamental frequency from the mode of the differences. This has given me pretty good accuracy. However, for lower notes, it is not working that great.

I was hoping that reviewing it would have helped, but no avail.
Following is the output from a sample FFT and autocorrelation.

Questions are:

Is the autocorrelation method correct? a) Double the input and pad with zeros. b) Forward FFT. c) Take the absolute square. d) Inverse FFT.
Am I looking for the right peak in the graph? I am assuming it should be the first peak.
Do the indexes of the autocorrelation graph relate to the frequency just as the indexes of the FFT? (see computeFrequency())
Is there any kind of preprocessing that I am missing, etc?
Is autocorrelation dependent on multiple harmonics being present?
Are there any additional techniques that I can apply to better implement autocorrelation? Or preferred fundamental frequency finding algorithms? (I have looked a little at cepstrum and was having similar problems.)

Thanks for any help received!

This is the output of an FFT. The peak is at index 47. I am positive I have implemented the FFT correctly. Index 47 matches up to E4 with a sampling frequency of 16,384 Hz and a padded FFT length of 1024.

I know, in this case, that, choosing the highest peak will suffice. However, there are FFT's that are not this kind.

Here is the output of the Autocorrelation. The first peak is at index 25 and the second at index 50.

Here is the code for the autocorrelation.

/**
 * complexData is overlapped 50%.
 * @param complexData
 * @return
 */
private Double[] autoCorrelation(final Complex[] complexData) {
    Complex[] toFFT = doubleAndPad(complexData);

    Double[] autoCorrelationAbsolute = new Double[toFFT.length];

    FFT.fftForward(toFFT);

    for(int j = 0; j < toFFT.length; j++) {
        //Same as Complex.mult(toFFT[j], toFFT[j].conjugate()) but simpler
        double square = toFFT[j].absoluteSquare();
        toFFT[j] = new Complex(square);
    }

    //Effective inverse FFT
    //forward FFT == inverse FFT for real numbers.
    FFT.fftForward(toFFT);

    for(int i = 0; i < toFFT.length; i++) {
        autoCorrelationAbsolute[i] = toFFT[i].absolute();
    }
    return autoCorrelationAbsolute;
}

Here is the code for doubling and padding

/**
 * Returns a copy of data with the size doubled and the
 * second half set to zero.
 * @param data
 * @return
 */
private Complex[] doubleAndPad(Complex[] data) {
    Complex[] doubledData = Arrays.copyOf(data, data.length * 2);
    //Set the second half to zero
    for(int i = data.length; i < doubledData.length; i++) {
        doubledData[i] = new Complex(0);
    }
    return doubledData;
}

This is what I use to calculate the frequency.

/**
 * Computes the frequency based off of the padded FFT in autoCorrelation
 * @param bin
 * @param data
 * @return
 */
public static double computeFrequency(int bin, AudioData data) {
    return bin * data.getFormat().getSampleRate() / data.getFftLength();
}

Didn't check the code, but the steps you describe in 1. seems correct.
The units of the autocorrelation are time and not frequency, and it is centred around $\tau = 0$. That means that in your $-1024 \leq\tau< 1023$. So you are looking at the wrong peaks.
Now as you can see the autocorrelation doesn't produce a sharp peak which means that it is not a good measure for your case. In general, the autocorrelation doesn't work well with narrow band signals.

I would suggest you try either the MPM or the YIN methods which are based on the autocorrelation, are extremely simple to implement and from my experience yield goo results.


Image Registration and optical flow Image registration algorithms can be classified into intensity-based and feature-based algorithms, I want to know When and Why should i choose intensity or feature based algorithm?

another question Is there a difference between image registration and optical flow?

Second question is easy: optical flow, more specifically dense optical flow, is an algorithm that takes two consecutive video frames and returns a vector field. For every pixel in frame 1 you get a vector showing where it moved to in frame 2.  You can also have sparse optical flow, which only computes the motion vectors for certain pixels, such as the feature-based approaches you've mentioned.

On the other hand, image registration is the general problem of aligning two images.  The problem of image registration can be solved using an optical flow algorithm. However, it can also be solved by other methods, such as template matching or phase correlation. Optical flow, on the other hand, can be used to solve other problems, such as detecting and tracking moving objects.

The first question is less clear cut, because there are situations when either approach will work. However, generally intensity-based methods are not applicable in the presence of a large scale change or a lot of motion between the frames.  Also, intensity-based methods tend to be slower.


Periodicity of signals $$x(t) = \cos(2 \pi t) \cdot u(t)$$
$$y(t) = x(t) + x(-t)$$Is $y(t)$ periodic. If so, what is the $T$?
$$x(t) = \sin(2 \pi t) \cdot u(t)$$
$$y(t)= x(t) + x(-t)$$
Is $y(t)$ periodic?

Yes it is.
Y(t), X(t) have the same periodicity periode for the sin signal.
The cos signal is not periodic.

That is of course assuming that u(t) is the unit step function with u(0) =1.


FFT of two sets of samples vs FFT of sum of samples If I have two sets of samples $a$ and $b$ and I have already calculated $A = FFT(a)$ and $B = FFT(b)$ is it possible to calculate $FFT(a + b)$ without performing another FFT?

Fourier Transform is linear, hence if $\mathcal F[x(t)]=X(f)$ and $a$, $b$ are complex numbers, then:

$$\mathcal{F}[ax(t)+by(t)] = a X(f) + bY(f) $$

So in your case, simply sum the results together.


Bit reversal on twiddle factors on inverse FFT I am currently working on a application which requires an inverse Fourier transform (IFFT) and am using a function library provided by Texas Instruments, in the datasheet they give an example on how to use one of their IFFT functions (shown below).

My question is: whilst I understand the need to perform bit reversal either on the input (decimation in time for an FFT) or the output (decimation in frequency for an FFT), I can't understand what's happening in the code given in the example. No bit reversal is performed on the input or the outputs only on the twiddle factors, is there some maths I am missing? Maybe something to do with the why the code accesses individual twiddle factors?

ps. The code is original from page 4-34 of this document: http://www.ti.com/lit/ug/spru657c/spru657c.pdf

void main(void)
{
 gen_w_r2(w, N); // Generate coefficient table
 bit_rev(w, N>>1); // Bit−reverse coefficient table
 DSPF_sp_cfftr2_dit(x, w, N);
 // radix−2 DIT forward FFT
 // input in normal order, output in
 // order bit−reversed
 // coefficient table in bit−reversed
 // order
 DSPF_sp_icfftr2_dif(x, w, N);
 // Inverse radix 2 FFT
 // input in bit−reversed order,
 // order output in normal
 // coefficient table in bit−reversed
 // order
 divide(x, N); // scale inverse FFT output
 // result is the same as original
 // input
}

The answer to your question can be found just in reading the comments in the C code that you quoted. The process looks like this:

The input vector x is passed through a function DSPF_sp_cfftr2_dit(), which is described as a radix-2 decimation-in-time forward FFT. It expects its input in normal order and writes its output in bit-reversed order. It expects its twiddle factors to be in bit-reversed order (which they are, as shown in the first two lines of the function).
The forward FFT result is passed through another function DSPF_sp_icfftr2_dif(), which is described as a radix-2 decimation-in-frequency inverse FFT. It expects its input in bit-reversed order (which it is, as the input is the output of the forward FFT above, which was in bit-reversed order). It expects its twiddle factors to be in bit-reversed order (which they are). It writes its output in normal order.

So, what you have is an input vector x that starts in normal order and runs through an FFT/IFFT round trip, resulting in a vector that is also in normal order (and should be approximately identical to the original vector x).


MIT 6.003 HW#8 Problem 4 - Fourier Coefficients of Triangle Wave In the mentioned homework, part of the solution involves finding the Fourier coefficients of the triangle wave.

The solution mentions that we can express this function as follows:

What does that multiplication signal means?
I approached the problem from a completely different angle of viewing the triangle wave as the integral of a piecewise constant function defined as follows:
$
f(x)=
\begin{cases}
 \frac{1}{2}&\text{if}\, -\frac{1}{4}\leq x\leq0\\
 -\frac{1}{2}&\text{if}\, 0\leq x\leq\frac{1}{4}
\end{cases}
$
Taking the integral similarly and then multiplying by $\frac{1}{2\pi ik}$ gives me a completely different answer of $-\frac{1}{4\pi^2k^2}[\cos(\frac{\pi}{2}k)-1]$.
What is wrong with my approach?

The symbol $\underset{T}{*}$ in the solution refers to ($T$-) periodic convolution, which is defined as
$$y(t)=\int_{0}^Tx_1(\tau)x_2(t-\tau)d\tau\tag{1}$$
where $x_1(t)$ and $x_2(t)$ are $T$-periodic signals. If $a_k$ and $b_k$ are the Fourier coefficients of $x_1(t)$ and $x_2(t)$, respectively, then the Fourier coefficients $c_k$ of $y(t)$ are given by
$$c_k=Ta_kb_k\tag{2}$$
Your approach to the problem is valid for $k\neq 0$, and your result is actually pretty close because
$$\frac{1-\cos\big(\frac{\pi k}{2}\big)}{4\pi^2 k^2}=\frac{2\sin^2\big(\frac{\pi k}{4}\big)}{4\pi^2 k^2}=\frac{\sin^2\big(\frac{\pi k}{4}\big)}{2\pi^2 k^2}\tag{3}$$
so there's just some scaling problem that I'm sure you'll be able to fix (check your definition of $f(x)$!). Note that you need to compute the coefficient for $k=0$ separately.


Equivalence between noise covariance of real and complex models of received signal I have a question on the equivalence of two models (complex vs. real) of a received signal in ISI channel. Yet its simplicity, I am a bit confused.

Model #1: Complex
Assume we transmit a complex symbol $\overline{a}$ drawn from a high-order QAM constellation. On a vector form, we can write the receive signal as
$$\overline{y} = \overline{G} \: \overline{a} + \overline{w},$$
where $\overline{G}$ is the channel matrix (it is real matrix of dimension $N \times N$) and $\overline{w}$ is the complex noise samples with zero mean and cross-correlation matrix $E\{\overline{w} \overline{w}^H\} =\sigma^2 \overline{G}$.

I would like to work with real-valued components instead of complex-valued components. That said, we can write the following real model of the received signal.

Model #2: Real
$${y} = {G} \: {a} + {w},$$
where $y = \begin{bmatrix}
\Re\{{\overline{y}}\}\\
\Im\{{\overline{y}}\}
\end{bmatrix},$ $G = \begin{bmatrix}
\Re\{{\overline{G}}\} & - \Im\{{\overline{G}}\}\\
\Im\{{\overline{G}}\} & \Re\{{\overline{G}}\}
\end{bmatrix},$ $a = \begin{bmatrix}
\Re\{{\overline{a}}\}\\
\Im\{{\overline{a}}\}
\end{bmatrix},$ $w = \begin{bmatrix}
\Re\{{\overline{w}}\}\\
\Im\{{\overline{w}}\}
\end{bmatrix}$, $\Re$ and $\Im$ are the real and imaginary parts of complex numbers. In this particular question, $\Re\{\overline{G}\} = G$ and $\Im\{\overline{G}\} = 0$.

As can be seen, both models are equivalent and can be used at the receiver for detection. My question is what is the cross-correlation of the noise in the 2nd model in terms of $G$, i.e. what is $E\{w w^T\}$ in terms of $G$?

My trial is as follows: $$E\{w w^T\} = \begin{bmatrix}
E\{\Re\{\overline{w}\}\Re\{\overline{w}^T\}\} & E\{\Re\{\overline{w}\}\Im\{\overline{w}^T\}\}\\
E\{\Im\{\overline{w}\}\Re\{\overline{w}^T\}\} & E\{\Im\{\overline{w}\}\Im\{\overline{w}^T\}\}
\end{bmatrix}$$
In general, how can I relate this to $G$? My expectation is/was $E\{w w^T\} = \sigma^2 G$!!!

Even though I doubt your noise model, here's how you can do it:

Given the complex noise vector $\overline w$, with the covariance matrix $E[\overline w\overline w^H]=\sigma^2\overline{G}$.

Let's separate $\overline w$ into its real and imaginary part (I use a different notation, because I think this is faster to write and read): $\overline w=w_R+jw_I$, where $w_R, w_I$ denote the real and imaginary parts of $\overline w$. Insert this into the covariance definition:

$$\begin{align}E[\overline w\overline w^H]&=E[(w_R+jw_I)(w_R^T-jw_I^T)]\\&=E[w_Rw_R^T+w_Iw_I^T]+2jE[w_Rw_I^T]=\overline{G}\end{align}$$

Let's compare this to the expression we get, when explicitely separating the real and imaginary parts:
$$\begin{align}E[ww^T]&=\begin{pmatrix}E[w_Rw_R^T] & E[w_Rw_I^T]\\E[w_Iw_R^T] & E[w_Iw_I^T]\end{pmatrix}=X\end{align}.$$

Hence, in order to define $X$ (which is the matrix you are interested in), we need to know $E[w_Rw_R^T], E[w_Rw_I^T]=E[w_Iw_R^T], E[w_Iw_I^T]$.

Let's do the comparison to $E[\overline w \overline w^H]$, which we have calculated before. We know that $G_I=0$ (i.e. the imaginary part is zero). Hence, $E[w_Rw_I^T]=0$, i.e. real and imaginary part are uncorrelated. What we dont know yet are $E[w_Rw_R^T]$ and $E[w_Iw_I^T]$. We only know that $E[w_Rw_R^T]+E[w_Iw_I^T]=\overline G$.

How to solve it? There are two possibilities:

1) (The more common assumption) We assume that real and imaginary parts of the noise have equal variance. This is often reasonble, e.g. when the noise is baseband noise and was obtained from downconversion a passband signal. Then, we have $E[w_Rw_R^T]=E[w_Iw_I^T]=\frac{1}{2}\overline{G}$.

2) Real and imaginary parts are not i.i.d. Then, the noise is not proper and the pseudocorrelation $T$ of the noise does not vanish. It is defined by $T=E[\overline w \overline w^T]$ and we calculate it to be (assuming $E[w_Iw_R^T]=0$)
$$T=E[\overline w \overline w^T]=E[w_Rw_R^T]-E[w_Iw_I^T]$$
In order to accurately describe a non-proper noise, the value of $T$ must be known. Then, we have two unknowns ($E[w_Rw_R^T]$ and $E[w_Iw_I^T]$) and two equations $$\begin{align}\overline{G}&=E[w_Rw_R^T]+E[w_Iw_I^T]\\T&=E[w_Rw_R^T]-E[w_Iw_I^T]\end{align}$$
from which you can get the unknown via $E[w_Rw_R^T]=\frac{1}{2}(\overline{G}+T)$ and $E[w_Iw_I^T]=\frac{1}{2}(\overline G-T)$.

Given this derivation, I still doubt your noise model. To me it does not make any sense, that in the equation $y=\overline G \overline a+\overline{w}$, the noise should have variance $\sigma^2\overline{G}$ due to ISI. Can you ensure that $\overline G$ is Hermitian?


What is the rms value of a signal and how to use it In some papers/books I found infos about a signal given as RMS value.

What is the RMS value and what information does it contain?

How can I compute it?

Moreover, if they say that the signal is a white noise with RMS $100$ band limited onto the interval $[0, 2000] Hz$, how can I use the RMS information to produce the signal they are talking about?
And if they say that the noise contributes to 5% of the signal RMS value?

RMS is the root-mean-square value of a signal.  For a digitised signal, you can calculate it by squaring each value, finding the arithmetic mean of those squared values, and taking the square root of the result.

Loosely speaking, it represents the average "power" of a signal.  In electrical terms, all signals with the same RMS voltage, when applied across a resistor, will produce the same amount of heat.

Note that Matlab has an RMS function Y = rms(X) in the signal processing toolbox.  if you don't have that toolbox, it shouldn't be hard to write your own function.

Since the RMS value is only related to the signal amplitude, knowing it doesn't help you to create the waveform.  However, if you have a function that will generate the desired waveform, you could measure the RMS value of the result, and multiply by an appropriate scaling factor to give you the amplitude you want.

If "the noise contributes to 5% of the signal RMS value", then presumably the signal would only have an RMS value of 95% of the total before you added the noise.  But beware that RMS values don't simply add.  Adding noise of RMS amplitude 1 to a signal of amplitude 10 will not necessarily produce a new signal of RMS amplitude 11.


How do you determine the properties of a differential equation? If I was given a differential equation of the form $$y'(t) + a(t)y(t) = x(t)$$ how would I be able to decide its linearity, time-invariance, and causality?

If a system is represented by a differential equation then it must be LINEAR.
If the coefficients of differential equation are function of time then it is time variant otherwise time invariance.
For causality and stability we will need its transfer function for any comment.


Interpolated FIR filter (from Oppenheim and Schafer's Discrete-Time Signal Processing, 3rd ed) [from Discrete-time Signal Processing by Oppenheim and Schafer, 3rd ed., p.196]
Two questions:

In this context, the filter with system function represented by Eq. (103) is called an interpolated FIR filter. This is because the corresponding impulse response can be seen to be the convolution of $h_1[n]$ with the second impulse response expanded by $M_1$.

This explanation as to why   Eq. (103) is called an
interpolated FIR filter is still not clear for me.
Why is the impulse response given by Eq (104) FIR?

The reason why "interpolated FIR filter" is an appropriate name is because the filter $h_2[n]$ has $M_1-1$ zeros between its (non-zero) filter taps, and the convolution with the filter $h_1[n]$ interpolates those $M_1-1$ values.

Eq. $(104)$ generally only represents an FIR filter if both $h_1[n]$ and $h_2[n]$ are FIR filters.