Signification of alternating signs of imaginary and real part of an FFT Let me start by saying that this question is motivated by data analysis that I am doing. It's part of a much larger set of signal processing operations, but in order to keep the question contained I'll ask a general question and hope, someone will be able to help me:

Is there any physical or mathematical characteristic that is known for a signal (containing only the odd Harmonics) whose FFT-imaginary and FFT-real parts have alternating signs ?

As a matter of fact, this alternation happens when I shifted my signal in the time domain. This signal is measured and contains many periods. As I don't know the phase shift between the excitation signal and the measured signal, I am not really sure where the correct signal begins because I don't have the real phase shift. So when I shifted the time signal with (3/4)*Number_of_bins_per_Period, I got a plausible time signal (that suits other operations I am doing on its basis). So I wanted to know if I can find the reason why this shift : (3/4)*Number_of_bins_per_Period, is the correct shift by analyzing its effect on the real and imaginary part of the signal. Before the shifting, the FFT of the signal has these real and imaginary part:
EDIT Please note that I am not asking about the fourier Transformation or why I am getting this behavior. This is clear. As I said, the shift I am doing (3 *Tr/4) I did it intuitively and it gave me afterwards the results I was looking for. I was hoping to find the answer, as WHY exactly this shift is the right one. And I thought maybe the hint would be the alternating signs. So why should be the sogns of the real and imaginary part alternating so that my shift is right . Is this a sort of finger-printing for a specific physical phenomena ?

This is exactly the expected result for a time delay: A delay in time is a linear phase in frequency; so you shifted the amount that would cause 180 degree phase shifts between your frequencies (Note that at your first frequency sample location you have a positive real and positive imaginary, which rotates to positive real and negative imaginary which is a -90° rotation, then follow the same process for your second frequency and you will see it is a +90° rotation; the phase change between the frequencies given is 180°). It is repeating at 180° due to the modulo property of phase, but you are really seeing the evidence of a linear increase in phase that is rotating past 180 degrees at each of your frequency samples where energy exists.

Consider a signal at given frequency going through a fixed length of cable; it will be delayed a certain amount in time, and given the frequency that time will relate to a certain phase shift between the signal at the input and output of that cable. Now consider a frequency that is doubled, the time delay is the same but there will be twice the phase shift between the signal at the input and output in this case.  This is the result of the Fourier Transform of a delay, the magnitude is one but the phase in linearly increasing at a negative rate.

So in your case we are seeing the resulting rotation of 180 degrees for the same distance between your frequencies tones (evidence of a linearly increasing phase).

If you plotted each frequency as a complex vector on an I-Q (real vs imaginary) diagram, as you change the delay of your signal in time, you will see the vector rotate accordingly with no change in magnitude. The higher frequencies will rotate proportionally faster.

Here is another view in case that was not clear describing the same effect and very specific to your case.  The top plot shows sine-waves in time for a first and third harmonic, and the delay at a fixed time offset such that the first harmonic has shifted -90° and the third harmonic as expected has shifted three times as much or -270°.

These last plots demonstrate your case exactly and shows the same thing on a complex plane which is consistent with your complex frequency outputs that have real and imaginary components. Note that I use the common notation I and Q for real and imaginary where "I" is "In Phase" and "Q" is "Quadrature Phase".

Here we see the same thing in that each frequency tone is a vector that rotates around the complex plane at a fixed rate; so the third harmonic is rotating three times as fast as the first harmonic (and the fifth is rotating five times as fast etc). After a given time delay, we see the diagram where the first harmonic has rotated 90° and the third has rotated 270°. To try to keep this clear, I first show the case where the first and third harmonic are in phase, starting at an angle of +45° (Such that I and Q are both positive as in your plots) and we see the exact condition where after a fixed delay the first harmonic has rotated -90° such that I is postive and Q is negative, while the third harmonic in that same time duration has I negative and Q positive. (-270° rotation.)

If we instead have a starting condition identical to your plots, where the first harmonic has both I and Q positive, but the third harmonic is 180° out of phase such that both I and Q are negative, we end up with the diagram below. We see the same thing where the first harmonic shifts -90° and the third shifts -270° after a given time delay. This is exactly the result from your plots: The first harmonic goes from I+, Q+ to I+, Q-, while the third harmonic goes from I-,Q- to I+,Q-.


What is the relationship of S-parameters to frequency response? I have some S-parameters. I'm reading about frequency response and the impulse response and step response.
If I have a set of S-parameters, then is the frequency response for a port (say output port 2, input port 1) simply $S_{21}$ for all the frequencies in my S-parameter file?

Let's use an example of a Wilkinson power divider that is used to take an input and split the power in half between two branches.
The ideal S-parameter matrix then looks like this
$$S = \begin{bmatrix}
s_{11} & s_{12} & s_{13} \\
s_{21} & s_{22} & s_{23} \\
s_{31} & s_{32} & s_{33} \\
\end{bmatrix} =
\begin{bmatrix}
0 & \frac{-j}{\sqrt{2}} & \frac{-j}{\sqrt{2}} \\
\frac{-j}{\sqrt{2}} & 0 & 0 \\
\frac{-j}{\sqrt{2}} & 0 & 0 \\
\end{bmatrix}
$$
Practically speaking, you are interested if the device is doing its job. In this case, you want to see if the power divider is indeed splitting the input power in half for the frequencies of interest. In a real divider, there will be deviations on how it performs over frequency with the hope that the best performance is achieved for the bandwidth of interest. To see this, you can plot the S-parameters which are a type of frequency response.
Take this example of a Wilkinson power divider designed to operate at 2.5 GHz

"Frequency response" is a general term for a particular output you want to see over some frequency range. As you suggest, you can plot the S-parameter values over their frequencies, which is a type of "frequency response".


Frequency of a periodic signal with distortions I would like to evaluate frequency of some unstable periodic signal coming from a detector:

The signal is registered continuously and may or may not be present (i.e. be periodic). The frequency and amplitude of the periodic signal can change but stays in some certain range. The shape of the periodic signal is more or less the same as shown on the second image. It may sometimes have distortions of higher amplitude and frequency.

I need to check 1) if there is a periodic signal and 2) what frequency it has.

Is there any good way to do it?

Edit 1:

I am not sure if FFT is reliable in this case. What do you think?

Frequency spectrum for all data points:

Spectrum with last 512 points

The spectra were calculated using Octave (online) based on example from Matlab site

N = 512; % length(all_data)
x = all_data(1:N);
Fs = N;
xdft = fft(x);
xdft = xdft(1:N/2+1);
psdx = (1/(Fs*N)).*abs(xdft).^2;
psdx(2:end-1) = 2*psdx(2:end-1);
freq = 0:Fs/length(x):Fs/2;
plot(freq,10*log10(psdx));

The signal looks like narrow band noise and may not be periodic at all.  The FFT will generally reveal a periodic signal that is in the presence of some noise.

In the FFT results, if you have a periodic signal, you should see peaks at the fundamental and its harmonics.  If the "periodic" signal is unstable (modulated in some way) the peaks will become less distinct.

You can try different window functions on your data to see if more defined harmonic peaks become visible in the FFT data.  For example compare results using a number of windows such as Gaussian and Blackman-Harris windows.

Also the choice of the window size (length of data) used when computing the FFT will affect your frequency resolution (bin spacing).  If you can roughly estimate the period of the signal, choose a data length that is a multiple of the estimated period so the bins line up with the signal harmonics.  In general a longer data set will provide better frequency resolution than a shorter data set, so if you don't know the period of the signal, go with a longer data set.


How to prepare and plot unequally spaced, irregular data to a contour plot or similar with MATLAB I've got a data set of hot-wire measurement velocity amplitudes at a given frequency bin (time data that has already been transformed to the frequency domain and I am just considering data for a given frequency bin) that have been sequentially recorded on specific locations in the xz-plane with a traverse-system. The problem is that the data does not represent a full uniform grid, but is sparsly sampled. Actually, the traversing process was performed in a way to get several locations within a given, more or less rectangular area, while avoiding hitting a circular obstacle within this area. Additionally, the regions close to the obstacle were sampled more densely as compared to regions farther away from the obstacle.

As a result the data set to start my problem with looks like the following  matrix (3 columns and e.g. 400 rows):

x z values (3 columns)

x1 z1 value1
x2 z2 value2
x3 z3 value3
.   .    .
.   .    .
.   .    .

(each value representing the velocity amplitude, e.g. at 500 Hz)

But the problem is:

The data was collected in some wierd order, e.g. the first 6 coordinates (x-component) are [-5; -5; -3; 10; 10; 0.8; ...] and the first 6 coordinates (z-component) are [0.2; 0.3; 5; 3; 2; -1; ...] ...
The data does not represent a uniform and equally spaced grid: all data (for example) lies within the range [-5 <= x <= 10] and [-2 <= z <= 6], but it is possible that there are 8 z-positions for one x-position, while at another x-position there are only 3 z-positions and at several other x-positions there are no z-positions at all (as within the obstacle region).
There is a greatest common divisor with respect to the grid resolution, e.g.    0.1   so all x and z-coordinates are positive or negative multiples of 0.1. This means it would be possible to generate a (probably) huge matrix (xz-plane) that includes all possible locations. I think this could help somehow.

What I basically want:

A simple contourf plot of the data within the xz-plane and areas that have not been recorded (obstacle region) are either interpolated or even better filled with NaNs (or Zeros) or something like that. Actually, I know the coordinates of the obstacle and might just draw it as overlay later.

Do you know how to get the desired result in a convenient manner?

I've been trying for ages now, but I just don't find a solution that solves my problem. I am pretty sure that it can be solved with some combination of reshaping, ndgrid/mesh, sorting or gridded interpolation with built-in (image  processing?) functions but I just don't get it to work.

I found a possible solution with the help of some collegue.

When gridis the matrix containing an x and a z and a value column (each having the same number of elements):

[X,Y] = ndgrid(linspace(min(grid(:,1)),max(grid(:,1)),150),linspace(min(grid(:,2)),max(grid(:,2)),150));
Z = griddata(grid(:,1),grid(:,2),grid(:,3),X,Y,'cubic');
contourf(X,Y,Z)
axis equal
scatter(grid(:,1),grid(:,2),10,grid(:,3),'wx');

So basically

a grid is generated with a relativly fine resolution (150x150 points) on basis of the min/max limits of the original data
the original data is interpolated onto this grid with a cubic interpolation method
contourf plot is generated and axis is set equal
the original data points are scattered onto the contourf plot to make sure that the interpretation of the data is somehow more reliable and one can see at once that the areas are interpolated.

I realized in between that the real difficulty is the obstacle area, as is gets wierdly interpolated although there should be no data, as none of the interpolation methods is used to such boundaries or optimized to those.


How to plot and normalize a spectrum in matlab I am given the coefficients of an FIR filter. I have created its impulse response through a convolution with a unit pulse. Now I want to plot the frequency spectrum. I perform zero-padding and windowing with a blackman-window. This does not seem to have a big effect on my result though. Still, I wonder: For some reason my amplitude of the resulting spectrum is not normalized correctly. Can somebody see what I am doing wrong?

fs = 44100;%sampling frequency
Ts = 1/fs; %time step
t = 0:1/fs:2-Ts; %signal duration

unit_pulse = zeros(length(t),1);
unit_pulse(1) = 1;

b = [1/8 1/8 1/8 1/8 1/8 1/8 1/8 1/8];
ir_filter = filter(b,1,unit_pulse);

nfft = length(ir_filter); %length of the time domain signal
%in order to obtain a good frequency resolution, pad the signal with zeros
%such that the length of the signal is a power of 2
no_zeros = (2^nextpow2(nfft)-nfft)/2;
padding = zeros(no_zeros,1);
signal = vertcat(padding, ir_filter, padding);
windowed_signal = signal .* blackman(length(signal));

spectrum_2side = real(fftshift(fft(windowed_signal))) ./ nfft;
f_2side = linspace(-fs/2, fs/2, length(windowed_signal));
spectrum_1side = 2*spectrum_2side(end/2+1:end);
f_1side = f_2side(end/2+1:end);

fig2 = figure();
hold on
grid on
plot(f_1side, spectrum_1side, 'b');
plot(f_2side, spectrum_2side, 'r');

This is a figure of what the impulse response looks like:

Short answer::

Use

freqz(coeff)

Long answer:

To plot the frequency response for a filter (which is a continuous function) you would use the Discrete Time Fourier Transform (DTFT), not the Discrete Fourier Transform (DFT) which is a discrete function. The FFT is an algorithm that implements the DFT. The difference between the DTFT and the DFT is that the DTFT is a summation from $-\infty% to $+\infty% while the DFT is a summation only over N samples. This is given in the formulas below.

DTFT:
$$X(e^{j\omega}) = \sum_{n=-\infty}^{\infty}x[n]e^{-j\omega n}$$

DFT:
$$X(k) = \sum_{n=0}^{N-1}x[n]e^{-j\omega_o n}$$

with k = 0,1,2 ... N-1 and $\omega_o = \frac{2\pi}{N}$

The DFT is a sampled version of the DTFT (the samples of the DFT are on the DTFT). You can approximate the DTFT using the DFT with zero padding (as that emulates extending the time domain to $\pm \infty$); the more zeros you add, the more samples will be interpolated in the frequency domain, and all samples will be samples on the DTFT. (So not just the next power of 2 as the OP has done, which is done to use the fft algorithm more efficiently, but enough to visually interpolate the entire spectrum. Note this does NOT increase frequency resolution but plots more samples on the underlying DTFT).

The freqz in Matlab/Octave and Python command also given above returns the frequency response and is the simplest approach.


Mid/Side encoding/decoding in frequency domain instead of time domain In the time domain let's say we have a signal $L_t$ and $R_t$. Then we can do mid/side encoding/decoding like so:

$$M_t = \frac{L_t + R_t}{2}$$
$$S_t = \frac{L_t - R_t}{2}$$

$$L_t = M_t + S_t$$
$$R_t = M_t - S_t$$

Let's say the we convert a block of $N$ values from $L_t$ and $R_t$ to the frequency domain by zero padding and taking their FFT to yield $L_f$ and $R_f$. Is it possible to do a similar encoding for mis/side in the frequency domain i.e. convert $L_f$ and $R_f$ to $M_f$ and $S_f$? So that the inverse FFT of $M_f$ and $S_f$ would yield $M_t$ and $S_t$?

Yes, it is possible. The discrete Fourier transform (DFT, of which the FFT is an efficient implementation) is a linear transformation, i.e.,

$$
a\cdot\text{DFT}\{x\}+b\cdot\text{DFT}\{y\}=\text{DFT}\{a\cdot x+b\cdot y\}.
$$

The same holds for the inverse transform.

As a consequence, you can directly convert the equations for the mid and side signals from the time domain to the frequency domain. Simply replace the subscript $t$ with $f$:

$$
M_f = \frac{L_f + R_f}{2}\\
S_f = \frac{L_f - R_f}{2}\\
\vdots
$$

From $M_f$, $S_f$ you can then compute the inverse discrete Fourier transform to obtain $M_t$, $S_t$.


How to get the Fourier series using Python's $\tt fft$ I Would like to be able to reconstruct every individual sinusoid that makes up a Discrete signal. I Have the following signal: (I am working in Python)

The signal is essentially an array with about 400 elements that varies with time. In Python after calling the fft function on the data

 dft= rfft(dat)/len(dat)  #real fft

I receive the figure below:

I am aware that I can use the result of the fft to obtain the individual Fourier series components, but I am unsure exactly how. The code I have yields me an array of magnitudes for each harmonic of the signal, how do I use that to get each individual component of the Fourier series?

I am no expert in this topic, but have some useful examples to share.

1 component example

To keep the i-eth Fourier component, you can zero the rest of the components:

n = len(y)
Y = numpy.fft.fft(y)
numpy.put(Y, range(0, i), 0.0)
numpy.put(Y, range(i+1, n), 0.0)
# Now Y holds 1 imaginary coefficient corresponding with the i-eth Fourier component

If you're looking for the Real signal corresponding to the chosen i-eth component, you can afterwards perform the inverse transform through ifft:

ifft = numpy.fft.ifft(Y) # the sinusoid you're looking for

Multiple components example

The more Fourier components you keep, the closer you'll mimic the original signal. This example shows what happens when you keep 10, 20, ...up tp n components.

Assuming x and y are your data vectors.

import numpy
from matplotlib import pyplot as plt

n = len(y)
COMPONENTS = [10, 20, n]

for c in COMPONENTS:
    colors = numpy.linspace(start=100, stop=255, num=c)
    for i in range(c):
        Y = numpy.fft.fft(y)
        numpy.put(Y, range(i+1, n), 0.0)
        ifft = numpy.fft.ifft(Y)
        plt.plot(x, ifft, color=plt.cm.Reds(int(colors[i])), alpha=.70)

    plt.title("First {c} fourier components".format(c=c))
    plt.plot(x,y, label="Original dataset")
    plt.grid(linestyle='dashed')
    plt.legend()
    plt.show()

...see how it mimics the original signal?


How to read Phase response graph? I'm a noob on DSP in general. But I'd like to understand this graph. It is the phase response of a band-stop filter, from Wikipedia:

I get the magnitude graph above: at some frequencies, the gain go down, than up, after the cut. But what about the phase?

I don't understand...

Why at the cutoff point, the frequencies (at that point) are drifted -90°/90° moved? What does it means? Is it the same phase as sin(90°) = 1?

I really don't understand that graph (visually). Is there any tool online maybe to show what is happening?

The frequency response of a filter (or, in general: of a linear time-invariant system) is a complex-valued function

$$H(j\omega)=M(\omega)e^{j\phi(\omega)}\tag{1}$$

where $M(\omega)$ is the magnitude response, and $\phi(\omega)$ is the phase response. The most straightforward way to interpret those two functions is to look at the way such a system responds to a sinusoidal input signal. For an input $x(t)=A\cos(\omega_0t)$, the response is given by

$$y(t)=AM(\omega_0)\cos(\omega_0t+\phi(\omega_0))\tag{1}$$

So, as you've correctly pointed out, the magnitude defines the filter's frequency dependent gain. The phase response determines the phase shift experienced by a sinusoidal input signal.

Assuming $\omega_0\neq 0$, Eq. $(1)$ can be rewritten as

$$y(t)=AM(\omega_0)\cos\left(\omega_0\left(t+\frac{\phi(\omega_0)}{\omega_0}\right)\right)=AM(\omega_0)\cos\left(\omega_0\left(t-\tau_p(\omega_0)\right)\right)\tag{2}$$

with the phase delay $\tau_p(\omega)$ defined by

$$\tau_p(\omega)=-\frac{\phi(\omega)}{\omega},\quad\omega\neq 0\tag{3}$$

From $(2)$ we see that a sinusoidal signal of frequency $\omega_0$ experiences a delay of $\tau_p(\omega_0)$. Note that the phase delay $\tau_p(\omega)$ is different from the group delay (cf. this wikipedia article).


Fundamental matrix rank Is there any way one can prove that $\text{rank}\big(F+F^{T}\big)=3$? where F is a fundamental matrix. If you have any idea, do let me know.
I am trying to solve this problem.

I am just trying to understand the first part.

This is not a proof, but maybe an intuition why this conjecture can be true for the points in general position. From the properties of rank we know that:
$$
\mathrm{rank}(F) = \mathrm{rank}(F^\top) = 2.
$$
Hence:
$$
\begin{align}
\mathrm{rank}(F+ F^\top) &\leq \mathrm{rank}(F)+\mathrm{rank}(F^\top)\leq 4.
\end{align}
$$
Because $(F+F^\top)\in \mathbb{R}^{3\times 3}$, $\mathrm{rank}(F+ F^\top)\leq 3$. To have any rank deficiency here (e.g. $\mathrm{rank}(F+ F^\top)< 3$), we should have an additional rank constraint on $F$, which we donot seem to have. The only additional constraint is that $F$ has 7 degrees of freedom and that $\det(F)=0$. These non-linear constraints are not sufficient to reduce the rank.


Face Classification. Is it OK to only use geometric features? I am trying to teach myself the basics of facial recognition. I see that some resources use just distances between some points on the face (e.g., distance between 2 eyes, eyes to nose, etc). Some others use more sophisticated methods that (in my opinion) take into consideration every pixel on the face (e.g., PCA). What is the difference between these approaches. Is it enough for me to just use the distances between points on the face?

What lies at the heart of pattern recognition and pattern classification is the selection of the correct features that is used in decisions. And the most important properties of correct features are 1-ease of measurement , 2-power of discrimination. So not every feature is the same. And finding the ones which yields the best performance requires research, creativity and expertise.

In your problem of face recognition, you should first find those features which most critically depend on particularity of the member of their class, but should also be easy to measure. Geometric features provide both powerful and easy to measure examples. Therefore you can rely on them.

But depending on your aim, you may consider using other features as well, if you want a better discrimination of members. So therefore it depends on your requirements of discrimination fault tolerance.


Matched filter for variable time scaling? Lets say I have some arbitrary waveform shape I want to look for. We'll use this as an example:

If my input is some mildly noisy signal, and I see something like this:

I can run that through a normalized matched filter and I should be able to easily identify this signal as my target waveform. However, lets say I see something like this:

Same general shape, but its at 1/2 the timescale. But I want to be able to detect this signal as well, and any other with the same shape but some unknown time scale(within some discrete limit).

I would like to be able to do this in real time, so simply decimating or interpolating the sampled date at different rates to see if that produces a match wouldn't be feasible. Is there any algorithm that could be used to quickly determine if any subset of my vector of samples contains a matching waveform to this target profile?

One first quick and dirty method could be to first store a sort of scale-space library of your original waveform $w(t)$, computed  at different scales (discretized), and then perform a "filter bank" of matched filters in parallel on the same observed signal. This proposition is a kind of  dual of the non-feasible solution you evoked. This way, you can save time from having the waveform scaled offline, and do not have to it on the signal. Then, you would probably have to detect the best match at a given time. Ambiguities may appear by superposition: at a given time $t_0$, the signal could look like $a w(t-t_0) + b w\left(\frac{t-t_0}{s}\right)$, where $s$ is a scale parameter.

To see how much better could an offline solution be, you may have a look at a scale-invariant matched filter. Thanks to you question, I just discovered the recent paper:  The Mellin Matched Filter, A. Monakov, 2015:

  In propagation channels signals are undergone translation and dilation
  changes. The Fourier and Mellin transforms are natural foundations of
  the analysis of wideband signals. Three interconnected problems, where
  the Mellin transform plays a key role, are considered in the paper: i)
  synthesis of a linear matched filter that is invariant to signal
  scale; ii) estimation of the signal scale; and iii) the wideband
  ambiguity function and its properties. The filter that is invariant to
  the signal scale and maximizes the output SNR is synthesized in the
  paper using the Mellin transform. The filter is named the Mellin
  matched filter. Introduction of the Mellin matched filter allows to
  give a new interpretation of the signal scale estimation and the
  wideband ambiguity function. Concept of the Mellin matched filter
  facilitates to produce a physical interpretation of the resolution of
  wideband signals and the time-frequency-scale uncertainties.

I have no idea about how it performs in practice.


What are the advantages of higher order Kalman Filters like EKF, UKF? Kalman Filter provides the optimal estimate of the states of a stochastic dynamical system if the system is linear, the measurements are also linear functions of states and the errors in system modeling and the measurements are Gaussian white noise. However, for a nonlinear system, the use of Extended kalman Filter (EKF) or Unscented Kalman Filter(UKF) provides a sub-optimal estimate. In general, for a nonlinear system, EKF gives a less accurate measure of covariance than UKF. In order to correct the covariance, higher order EKF have been proposed. Similarly, higher order UKF also leads to a more accurate covariance. My queries are:

Apart from improvement in the covariance, do higher order filters provide any other advantage?
Can we capture skewness or kurtosis using these high order filters?

They should improve the convergence properties of the state estimate, i.e. if you have an initial value far from the trajectory, they could improve the region of attraction. They should also reduce the bias in the estimate.
In the EKF and UKF, there are no direct ways of including higher order moments of distributions, in both cases the state distributions are assumed to be Gaussian. Presumably, you can approximate non-Gaussian distributions with additional modeling. For (some) higher order approximations, it should be possible to include higher order moments. [However, as far as I understand, the additional computational complexity is usually not worth the performance gain (if any), since there are so many sources of uncertainty when modeling a system.] The go-to solution for determining the distributions of state-estimates is Monte-Carlo simulation.


Converting raw I/Q to dB I am getting I/Q data from a software-defined radio.  I want to do some stuff on signals in the data, but only if it exceeds a certain range.  What is the general procedure to get dB (dBm, or anything) from this kind of data?  Programs such as SDR# do it, but I don't know what exactly they do so that I can imitate them.

Basics

The amplitude of an IQ signal is just the vector magnitude, $\sqrt{I^2 + Q^2}$.

The power of an IQ signal is the squared magnitude, $I^2 + Q^2$.

When you see a logarithmic (dB) meter, it is usually measuring the log of the power, i.e. $10 \log_{10}(I^2 + Q^2)$. (This can also be calculated as $20 \log_{10}$ of the amplitude, but unless you already have the amplitude that wastes a square root operation.)

Units

Remember, dB is a relative figure. If you just take $10 \log_{10}(I^2 + Q^2)$, then 0 dB corresponds to an amplitude of exactly 1. If your hardware driver takes the usual floating-point convention of the absolute extreme sample values being from −1 to +1, then you can say that your dB power values are dBFS — decibels relative to full scale. Any signals stronger than that level will be clipped, distorting the signal.

dBm, decibels relative to one milliwatt of power, just uses a different reference level. You can convert to dBm just by adding or subtracting the proper calibration value from the dBFS value — but you need to know that calibration for your hardware at the frequency of interest, such as by measuring it (using a signal source of known power output); it is impossible to perform that calibration purely in software since the digital samples are just numbers with no inherent units.

(A mistake I've seen is to refer to the sample values, or parameters which scale according to them such as an amplitude threshold, as being in “volts”; this is complete nonsense unless your ADC (and other hardware) is actually calibrated to one volt. This is unreasonably large for a radio receiver.)

Practical application

If you are just looking to ignore signals that aren't sufficiently strong (this is known as carrier squelch or power squelch), it doesn't matter what units you use, or even if they're logarithmic or linear, because you're just doing a greater-than comparison. The only critical component is that you start with $I^2 + Q^2$ (as opposed to, say, $I + Q$, which would be just plain wrong).

Note on bandwidth you probably don't need to read

It may also be relevant to note that if you filter a signal, you are by definition removing some of the signal power, so the measurement will be smaller.

In particular, a FFT (such as is the primary visual display in tools like SDR#) can loosely be thought of as a large collection of extremely sharp filters; each output “bin” collects some fraction of the input power. Accordingly, the power in each bin depends on the width of the bin. If you divide by the width of the bin in hertz (that value being $\text{sample rate}/\text{FFT length}$) before taking the logarithm, then instead of dB power, you measure dB power spectral density, which has the advantage of being independent of the FFT bin width if the features you care about are wider than one bin (e.g. a wideband modulated signal); if they are narrower (e.g. pure tones) then the power value is more useful.


Moving from books to implementations So after a good while of digging and practicing, I think I've managed to gain a decent intuition for some of the very core concepts related to linear systems, filter coefficients/parameters, FIR/IIR filters, the DFT/FFT, analyzing spectra and using statistical measures for random signals (stationary,non-stationary). I have a fairly good idea of magnitude and phase, and implementing and analyzing an STFT.

BUT... I struggle quite a lot with implementing something I read in a paper, especially if it's to do with more practical applications that need adaptive filters or estimating filter parameters. I can follow along for most of the math and understand what is being said for the most part, but find myself very easily lost for weeks to no end.

I'm curious if anyone out there has had similar issues, and if there were any solutions that were particularly useful for getting more comfortable with the core concepts. I honestly feel it would help a lot to have some kind of supervisor/TA that I could go to with my questions. Perhaps there are tutoring/mentoring websites out there? StackExchange is honestly really great, but sometimes I either find myself struggling to ask the right question, or not knowing enough to understand the given answers.

In my opinion,

DSP is a tool. It is interesting in of itself but it is most interesting when used to solve real problems.

Find an application, problem, or community that interests you.

Take medical imaging as an example. Understanding the slice theorem is necessary but so is knowledge of anatomy and physiology.

If something like RADAR interests you, knowledge of RF propagation is essential.

Do you like music? You can use DSP to like it more.

Look at other books beyond  signal processing. Think about how you could use DSP in your hobbies.

Practicing DSP is about balancing function and constraints.


What's a Normalized function? I'm studing Representations of Random Processes and the book talks about Orthonormal functions, but doesn't make it clear what is it.
I was able to realize that a set of functions are orthonormal if they are normalized and orthogonal.
I'm familiar with the concept of orthogonal functions, but I couldn't find out what exactly is a normalized function in that case, since is a widely used term in math.
I'm studing using the following book: Detection, Estimation, and Linear Modulation Theory Part 1 by Harry L. Van Trees

Similar to a a vector that is normal when it's magnitude is 1, which is the inner product with itself. In the same fashion, a function is normal over a defined range $[t_1,t_2]$ when the root of the integrated squared magnitude is one (if it's normal the root is trivial but this shows the same form as derived from an inner product):
$$||x(t)|| = \sqrt{\int_{t_1}^{t_2}|x(t)|^2dt} = 1 $$
Normalization of functions is done with a simple scaling: If the above formula results in a number of than one, then the entire function is scaled by that number.


Scratch Detection On Noisy Textured Surface I wish to extract long thin scratches from the noisy background image.I have tried several ways such as stacking smoothed image with input image to remove some of the noisy background patch before applying blob detection filtering on those images to get the long thin scratch. However, overkill will happen due to the noisy background patches.

Any good suggestion or guidance on other things to try?

I tried PHOT (suggested by @applesoup) using the unofficial implementation given under https://github.com/thinkng but it did not work for these image sets. Maybe one needs to further investigate or tweak the algorithm a bit.

First, if you have sufficient amount of data, I do not believe that one could easily outperform a good deep architecture in the task of texture defect classification. The literature is full of these data driven methods. Some examples are here: [1], [2], [3], [4], [5].

However, due to limited amount of data, I tried a classical approach which seems to output somewhat acceptable results. In essence, the method involves:

Compute the multi-scale difference of Gaussian (DoG) image. Call this $D=\{D_i\}$, where $D_i$ is a single DoG:
$$
D_i = I * g_i^1 - I*g_i^2
$$
where $I$ is a single channel input image and $g_i^1$ and $g_i^2$ are two Gaussian kernels with $\sigma_2>\sigma_1$. $*$ denotes convolution. As your defects are line-like, such edge enhancement should ease the job of any subsequent analysis. The resulting DoG images look like:

You could already see that the defective areas are a bit emphasized. I show here only one of the scales, whereas multi-scale convolutions form a scale pyramid (see scale-space theory). That is what $D$ is.

Perform a curve (line) detection on this DoG pyramid $D$ by using the algorithms provided in one of my previous responses. The raw response of such a subpixel line detector on a single level of the hierarchy:

Even though the images are messy and cluttered by many false positives, the lines we are looking for seem to be there. I avoid showing all the detected lines as this would severely harm the perception.

The final stage is to filter these curvilinear structures to arrive at the desired scratches. While many filtering techniques can be used, I just let the lines that are long enough survive. The output after all these stages looks like:

This is not ideal, but probably close to what you have been looking for. I also did not spend much time to tune the parameters.

Note that this approach is solely based upon a geometric modeling of the problem. For different types of scratches these assumptions won't hold. Nevertheless, I thought it might be a good starting point.


Number of microphone limitation for Spherical Harmonics based DOA estimation I am following this book "Theory and Applications of Spherical Microphone Array Processing", as well as this paper on DOA estimation "Eigenbeam-ESPRIT for DOA-Vector Estimation".

My dataset is composed of a small number of channels (4) from the Eigenmike32. The book states "all spatial sampling schemes require at least $(L_f + 1)^2$ microphones to sample a sound field of order Lf without aliasing.", where $L_f$ is the order of the sound field (3.4 Spatial Sampling). Both the paper and the book use 32 microphones in their implementations.

First, I do not find any information on what is the order of the sound field. It sounds to me that this is a signal-dependent quantity, much like a sampling rate is lower bounded by the highest frequency of the signal you want to sample. Did I get this right? Can anyone elaborate please?

Second, does that mean that I can do very little with 4 microphones? Is there a lower limit on the number of microphones for a DOA estimation and for a number of active sources estimation?

What other, none neural network, methods can I use to estimate the number of active sources and DOA with (no necessarily by the same model)?

First, I do not find any information on what is the order of the sound field. Can anyone elaborate please?

$Y_{lm}(\Omega)$ is the spherical harmonics with mode $m$, order $l$ and frequency $\Omega=(\theta,\phi)$ (dictating angles of arrival, ex: azimuth, spherical, elevation, etc.). According to [1], we have the following formula:

$$Y_{lm}(\Omega) = \sqrt{\frac{2l + 1}{ 4\pi} \frac{(l-m)!}{(l+m)!} }P_l^m(\cos \theta)e^{j m \phi}$$
where $P_l^m(\cos \theta)$ is the associated Legendre functions of degree (mode) $m$ and order $l$. The above serves as a basis of spherical harmonics (similar to Fourier expansion).

  Second, does that mean that I can do very little with 4 microphones? Is there a lower limit on the number of microphones for a DOA estimation and for a number of active sources estimation?

The classical rule of thumb is that given $q$ sources, you need at least $q+1$ antennas to resolve uniquely your sources. However, this is really algorithmic dependent. For example, the paper in [2] (which formulate a 2D AoA DoA estimation) describes shift invariant methods with different lower bounds. So, in principle, the more antennas you have the better resolvability and lower MSE you will have, since the Cramer-Rao bound of DoA estimation is proportional to $\frac{1}{N}$ where $N$ is the number of antennas.

  What other, none neural network, methods can I use to estimate the number of active sources and DOA with (no necessarily by the same model)?

You can use subspace methods like MUSIC, or co-variance fitting. Some generic methods aim at directly optimizing the Maximum likelihood problem like Expectation-Maximization or IQML for specific antenna configurations. A 2D IQML approach was derived in [3].

[1] Rafaely, Boaz. Fundamentals of spherical array processing. Vol. 8. Berlin: Springer, 2015.

[2] Bazzi, Ahmad, Dirk TM Slock, and Lisa Meilhac. "Single snapshot joint estimation of angles and times of arrival: A 2D Matrix Pencil approach." 2016 IEEE International Conference on Communications (ICC). IEEE, 2016.

[3] Bazzi, Ahmad, Dirk TM Slock, and Lisa Meilhac. "Efficient maximum likelihood joint estimation of angles and times of arrival of multiple paths." 2015 IEEE Globecom Workshops (GC Wkshps). IEEE, 2015.


Ideal reconstruction after down sampling The signal $x_a(t) = \cos(2\pi450t)$ is sampled.

F = 450
Fs = 1000 Hz
f = F/Fs = 450/1000   // Sampling theorem is fulfilled
x(n) = cos(2*pi*(450/1000))

The signal is then down sampled with a factor 3.

fNew = f*3 = 450*3/1000 = 1.35
xNew(n) = cos(2*pi*1.35)

Now the signal is prepared. How to make an ideal reconstruction using 1000Hz?

The signal can not be reconstructed because in the new sampled signal the sampling theorem is violated.  Aliasing would occur and the signal reconstructed would  be a lower frequency than the original signal.


FFT frequency bands and filtering I don't really understand FFTs, and filtering neither. I hope I can explain my confusion in a way that is not too confusing for you to understand.

consider the following example: I have an EEG signal, and am interested in the frequency bands from 1 to 40 Hz.

if I apply a LP filter at 40, and a HP at 1, my signal of interest would be distorted around the beginning and end because of the filter's transition band. right? why does it have a transition band at all, and how does it work? is a filter a sum of sine waves in the specified frequency range? (sorry, seems like a basic question, still I couldn't find an answer to it)
now to get the FFT of different frequency bands in the 1-40Hz range, I could either first bandpass filter the specific ranges and apply FFT to each (e.g. to 1-4,4-8,8-12,12-20). however, this would lead to a distorted result in each band, because of the transition bands, right?

if on the other hand, I do FFT over the whole 1-40 range, and then just 'extract' the results from 8-12 etc, the results would be more accurate? and how would the extraction mathematically work, is it just the mean of all values at the frequencies between 8 and 12?

I hope you can help me understand these things

Regarding the transition band of the lowpass and highpass filters; an ideal filter would indeed have zero transition bandwidth but it would also have infintely many filter coefficients in its impulse response. In practice you will have a finite length approximation to that infinite length impulse response, as a result the practical filer will have non-zero transition bandwidth. In case of an FIR (finite impulse response) implementation, the longer the impulse response the narrower will be the transition bandwidth.  As a result of this nonzero transition band of the finite length FIR filter, the resulting frequency response will display deviations from the ideal response through the transition band, you should (and can) take this into account. You can also consider using IIR based recursive filters, but that will still have a transition band, in additon to the nonlinear phase distortion being introduced, which doesn't exist for a typical (symmetric) linear phase FIR filter.

Regarding your second question on the spectral analysis of the filtered signal through its FFT, you first obtain a long enough signal observation window and then take a dense enough FFT of the filter output, so that you can focus on any region withn 1 Hz to 40 Hz depending on your interest. Assuming that your sampling rate is $F_s$ samples per second, then an $N$-point FFT would provide you an equivalent analog frequency resolution of $$\Delta f = \frac{F_s}{N}$$ where each frequency bin, in Hz, would be located at $$f_k = \frac{F_s}{N} k ~~,~~ \text{for}~ k=0,1,...,N-1$$ Note that FFT of real valued input signal would have conjugate symmetric complex output, so only the first half (roughly) carries relevant information from $f=0$ to $f=F_s/2$

Note that your spectral resolution is fundamentally limited by your observation window length i.e., FFT length does not improve frequency resolution, but merely interpolates the already resolved spectrum.


Bpsk demodulator implementation For a project I'm trying to implement a BPSK demodulator on some random wave files, I'm getting a few difficulties:

How do I know my bit rate, i.e. how many bits each 'real data bit' I'm having.
How can I figure out the size of my real data.
I don't know exactly the corner frequency $f_c$, is it OK for me to estimate it for the carrier wave?

You want to take a wav file, run it through a BPSK modulator and I'm guessing (from your #3) to plot some waveform or something. So you should find a wav file like handel.mat which is built into MATLAB, and convert it to binary. This may require some decisions like how many bits you choose to represent each sample if its not already something you can work with.

All a modulator is doing is mapping bits to symbols. For BPSK, a single bit gets mapped to a single symbol (Bits 0 --> Symbol -1, Bit 1 --> Symbol +1). So now you have mapped to symbols.

Now to make a plot you can just look at the equations for BPSK modulation and plot for some $f_c$ that you choose. Or, what I think you may have meant, is to plot the baseband signal in which case you then take your symbols, upconvert them, and pulse shape them (you'll have to choose some pulse shaping filter, check out rcosdesign() function in MATLAB).


How many bytes would I require to store few seconds of fan controller PWM signal? How many bytes would I require to store few seconds of fan controller PWM signal?

Such PWM is typically around 25Hz, so I think I need to sample it at 50Hz, right?

But what about the bit depth? Since it's a pulse wave (and such takes only values 0, 1 I believe) then 2 bits does suffice, right?

So is the computation then:

50 times per second * 2 bits per sample = 100 bits per second = 12.5 bytes per second of signal?

The easiest answer is: find the maximum frequency in the signal, and double that. In the case of PWM, call $T_p$ the duration of the shortest pulse possible. The signal bandwidth is in theory infinite, but in practice, frequencies above $10T_p$ will contain very little energy. So, you can safely sample the signal at $20T_p$ samples per second.

If you need better fidelity, then sample faster, accordingly.

Now, a PWM signal typically has only two amplitudes, so you can use one bit per sample, or $20T_p$ bits per second.


Where should triangle, sawtooth and square waves start? I am reading a text on DSP, which shows generation of sin, cos, triangle, sawtooth and square waves.  Generated cos waves start at 1 and sin waves start at 0, which makes sense (cos(0)=1, sin(0)=0).  Generated triangle waves start at 1, square waves at -1, sawtooth at -1.  Wikipedia gives a number of constructions for square wave, which would lead to starts at 1, -1 or 0.

Is there a convention in the DPS/electrical engineering/oscillosope world on where triangle, sawtooth and square waves should start?

If you look at it on the surface, all periodic waves have phases associated with them such as $\sin(\omega t + \phi)$ or $\cos(2\pi f_0 t + \theta)$ which has phase angels of $\phi$ and $\theta$ radians respectively. Therefore, they do not necessarily begin at their peaks at $t=0$ always. Indeed at $t=0$ above waves assume the values of $\sin(\phi)$ and $\cos(\theta)$. And you can therefore shift them by amounts of $\phi$ and $\theta$ left or right to make their phases zero with respect to display graphics plotting systems.

Oscilloscopes make this adjustment by their trigger levels. For example, to set the peak of the signal $A \cos(2\pi f_0 t + \theta)$ at the origin ($t=0$) then you shall adjust your trigger level to about $A$, assuming that the trigger operation will start plotting the wave at that amplitude positioned at the origin of the coordinate grid.

The same should be true for other waveshapes such as sawtooth, triangle, square or many unnamed periodic possibilitied.

I cannot see if you have a deeper question here, but on the surface the the answer is as above.

Plot of $\text{sgn}( \sin( 2\pi t ) ) $ :

And plot of $\text{sgn}( \cos( 2\pi t ) ) $ :


What is the basic concept of Fourier transform I use discrete Fourier transform for digital image processing purposes , but I don't understand basic concept behind it. For example :

What information exists in frequency domain?
What is difference between spatial domain and frequency domain?

What information exists in frequency domain?

As JasonR says in the comment, The frequency domain
is a different view of the same data.

No new information is created, it just takes the "spatial" domain data (the image pixel values and their locations) and re-presents it as the coefficients of complex exponentials (sines and cosines).

  What is difference between spatial domain and frequency domain?

The spatial domain is the domain of the image: the pixel values are located at particular positions in the image --- they are spatially distributed, usually in a regular grid.

The frequency domain takes this same data and finds any underlying periodicities (sine waves and cosine waves) in the spatial data, and their amplitudes and phases (spatial offsets).

For example, suppose I have the following image (Scialab, not matlab):

N = 100;
x = [1:N];
y = [1:N];
phi = 2*%pi*0.0987298374*ones(N,1)*x + 2*%pi*0.033102974*y'*ones(1,N);
im1 = sin(phi);

Which looks like (appropriately scaled to the grey scale values):

Then the FFT of this is:

(again, with appropriate scaling).

The frequency domain version shows up the periodicities of the spatial domain as a small number of large coefficients.


Idea for Noise Level Estimation / Automatic Thresholding in the Presence of Peaks I have the fft of some signal, and want a rough estimate of the noise level in order to choose an appropriate threshold for our peak detection algorithm. In general, the fft contains mostly noise with a handful of peaks (which are usually pretty high compared to said noise). For reference, I've attached a screenshot of a pretty typical fft:

Now to my proposed algorithm for noise level estimation, which is based on the following assumptions:

What i want is the mean and standard deviation of the fft
If the fft contained only (gaussian) noise, the median would be very close to the mean
Peaks shouldn't affect the median too much
The mean and standard deviation are independent of the order of the samples
Therefore i can sort the fft and still get the same result

So the resulting algorithm looks like this:

Sort the fft (ascending)
Calculate the mean and standard deviation of the sorted fft, but stop as soon as the mean exceeds the median
Use the calculated mean and standard deviation to choose the peak detection threshold
Profit?

I've tested this algorithm on some typical signals and the results seem to be pretty decent (keeping in mind that i don't need the exact noise level, but just something robust to choose an appropriate threshold). Sorting is a little expensive of course, but i don't expect the noise level to change too much, so i don't need to run it for every fft.

That being said, the algorithm itself feels kind of "wrong", because i don't ask the actual question "What would the mean and standard deviation be without the peaks?" but rather "What would the standard deviation be if the mean was X?".

So what are your thoughts on this?
Are there better algorithms for this kind of problem? (i bet there are)

I'd do some small adjustments to your idea (You really nailed them).

Assumptions

The Signal Model - Signal + Additive White Gaussian Noise (AWGN)
Probably we could generalize it more but this is beyond the scope of this question.
The DFT of the signal contains Peaks with relatively small roll off
This is important as we're almost saying the Signal is a combination of sparse number of Harmonic signals. As the lobes width will hide the noise data. So the ratio between the observation window, signal lobes width and sampling rate means that the DFT of the data creates sharp peaks.

Estimation

Under the above assumptions we know that:

Noise is AWGN in the DFT (As Gaussian Vector stays Gaussian under linear transformations). It also suggests that the Mean Value and the Median are similar (Realization of Symmetric Distribution).
Without the Signal Data the Mean and the Median value of the beans should be similar.
Number of bins contaminated with signal (Funny to write this, ha?) is sparse.

So we can sort the bins, but sorting and evaluating the mean on sub set will create bias unless we make sure the sub set is symmetric around the Mean / Median.
So, I'd sort the data and keep bins in the 25-75 Percentile range only. I'd use Mean on those.


Finite sequence input to DTFT i'm studying the practical utility of Fourier transforms and i have some questions. I hope to receive answers in layman terms.

1) Does the DTFT take only infinite input sequences?

2) If i apply the FT to a step function i get an infinite band spectrum. If i sample the step function and apply a DTFT does i get replicas of that infinite spectrum (obtained from FT) in the DTFT spectrum ?

I doubt that Hilmar will agree (maybe he will) but I will offer this answer as a counter to the one that has already been (prematurely) accepted by the questioner.

Hilmar's answer is wrong.  It is factually wrong and analytically wrong.

    1) Does the DTFT take only infinite input sequences?

  No. There is nothing that prevents you from applying the formula to a finite sequence.

That's sorta true but it must be said that the "formula" is explicitly of an infinite sequence.  However, there is nothing in that formula that says any particular sample must take on any particular value (like non-zero) other than a finite value.  There is nothing that prevents all but a finite subset of samples having zero value.

The DTFT is:

$$ X(e^{j\omega}) = \sum\limits_{n=-\infty}^{\infty} x[n] \, e^{-j\omega n} $$

That's what the formula is and all of the $x[n]$ have a defined value, even if it's zero.

  However, a finite time sequence has an infinite spectrum, so you can't sample it without some amount of aliasing

No one is sampling it (again).  It has already been sampled.  The question is: What happens if the length of the sequence (which is, in the DTFT, infinite in length) is shortened to a finite length?

And the answer is: It depends on the manner in which the sequence is shortened to a finite length.  Or, if the sequence was never infinite in length and later applied to the DTFT, how are the infinite number of samples outside of the original finite length defined.

But if you do define all of those samples outside of your original finite length to be zero, that is saying:

$$ x[n] = 0 \qquad \text{for } x<0 \text{ and } x\ge N < \infty  $$

Then the DFT is equivalent to:

$$ X(e^{j\omega}) = \sum\limits_{n=0}^{N-1} x[n] \, e^{-j\omega n} $$

And then it starts to look like the DFT.  If $X(e^{j\omega})$ is sampled at $N$ equally-spaced values on the unit circle including $X(e^{j0})$, then it is equivalent to the DFT

$$ X[k] = X(e^{j\omega}) \Big|_{\omega = \frac{2 \pi k}{N}} = \sum\limits_{n=0}^{N-1} x[n] \, e^{-j2 \pi k n/N} $$

But something misleading in Hilmar's answer must be corrected:

Applying a finite-length sequence to the DTFT has nothing to do with aliasing.  Nothing at all.


What is the filter with the less phase shift? I have to analyze a dynamic signal but there is too much noise so I applied low pass filter but then there is too much phase shift.So what is the most reactive filter I can apply to my signal ?

Best regards

Technically, the most reactive filter is the all-pass filter with a gain of 1, this filter has no phase shift at all. But it is not a really useful filter.

Here's what you need to take in account :

1 - IIR versur FIR

An IIR filter will have less phase shift than a FIR filter for the same cut-off frequency

2 - Cut-off frequency

For an IIR low-pass filter,  a filter with 0.1 Fs cut-off frequency will have a higher phase-shift than a  with a 0.2 Fs cut-off frequency. An all-pass filter i.e. a gain of 1 has no phase shift.

3 - IIR filter order

For an IIR, the higher the order, the higher the phase shift.

4  - FIR number of taps

Assuming you a linear-phase low-pass FIR filter, the more taps you have the higher the phase shift. If you want a high cut-off frequency you can get by with less taps. If you want a really low cut-off frequency you will need more taps thus more phase shift.

5 - Minimum phase-filter

There is a class of filter called minimum-phase filter. They are the most reactive (fastest decay) for a given cut-off frequency according to this reference

https://www.dsprelated.com/freebooks/filters/Minimum_Phase_Filters.html


Filtering random noise from a signal using excel I was given an excel file, and I was asked to filter the noise so that signal would be as close as possible to the original one, using also excel.

I kinda understand that I need to do $(x(n)+z)-z$ ($z$ being the noise), but I don't understand how to find the exact formula to find the noise when the noise itself is generated using rand(). We were given option to use DFT FFT FIR IIR, so I guess the solution would be that, but my lecturer didn't really teach us about FIR and IIR (due to corona I guess, not enough time) and the lecture never really cover the part of using DFT FFT to signal with noise. So, how to solve this problem using excel?

Signal given (without noise) is as below:
$$x(n) = \sin\left(2\pi \frac{8}{128}n\right)$$
Information Frequency: 8Hz

Sampling Frequency: 128Hz

I kinda understand that I need to do $(x(n)+z)-z$ ($z$ being the noise)

There goes an engineering joke that is just "If you know you've got noise in your signal, why don't you just subtract it!".

It's a joke, because the main property of noise is that you do not know its value, so you cannot just subtract it!

  I don't understand how to find the exact formula to find the noise when the noise itself is generated using rand().

Exactly. That's because (if rand is actually random), there is no formula.

So, you cannot subtract the noise like that completely.

What you can do is filter (and all the options you've been given are filters, or filter banks, with DFT being identical to the FFT, not quite sure you're fully grasping your assignment there)!

So, you want to build a filter that lets through nothing but your signal of interest – your 8 Hz sine wave.

(In the following: italics denote terms that your lecture most probably introduced. If you don't know them – look them up! Relevant to future understanding, and quite possibly to your grade.)

As you know (or should know, sounds like that was the whole point of the lecture), filtering is applying a linear system to a signal. And I'd recommend you consider a FIR first, so that filtering is just convolution of the signal with the impulse response of that filter (i.e. the filter coefficient vector gets convolved with the signal).

Because of linearity, you can look at how the convolution with that filter affects a) the pure signal and b) the pure noise, because if you add the results, you'll get the same as if you filtered signal added to noise.

So, you're looking for an impulse response that maximizes the amplitude of the sine passing through it. Convolution has a formula, and hint: if you use the same (just ever so slightly modified) as the thing you convolve with that you want to maximize, then you get a high amplitude.

This is as far as I want to do your homework.


Is acceleration noise modelled differently in EKF and UKF Kalman Filters? In a lecture on basic Kalman filter, I came across the following assumption about acceleration noise.

Every component of the noise vector $\nu$ is itself a product of time and acceleration values. Nothing is known about the accelerations $a_x$ and $a_y$ but each component of $\nu$ is assumed to have zero mean and every pair of components $\nu_j$, $\nu_k$ has some known variance $\sigma^2_{\nu_j\nu_k}$. As a consequence the last term in each of the four equations is set to zero (not shown here). The covariance matrix Q contains time though as in the screenshot below.

In a latter lecture on Unscented Kalman Filter the assumption was changed as below.

Here the components of noise vector $\nu$ do not contain any effect of time. These time-free components are assumed to have mean zero and covariance $Q$. As a result their effect on the mean state vector $x$ cannot be set to zero. Also $Q$ is independent of time.

Is it standard practice to make different kind of assumptions?

Absolutely!

Every (extended / unscented) Kalman filter starts with a signal model.  The state update equation in the second image is a very different beast from the little you show of the first model in the first image.

Without seeing the complete model, it's hard to say whether it's all correct, but if the signal models are different, then certainly the state update and output equations for the filter will be different.

For example, this paper shows two EKFs for frequency tracking that use two different signal models.

Regarding incorporation of the $\Delta t$ terms in the acceleration, that does seem a bit odd.  Usually those terms form part of either the state transition matrix or the input matrix. See this answer for one way to do that for a simpler model.


In computer stereo vision, what does it mean to "look in"? In trying to understand this answer to my earlier question, can someone explain why it is important that the cameras are not in the same plane when the focal distance is large compared to the separation between the cameras and the object?

What is meant by "looking in"?

Does this change if the number of views/cameras increases?

"Looking-in" simply means that the cameras will need to be angled such that their principal axes cross over at some point in front of the camera (rather than being parallel as is often assumed for stereo work)

The reason for doing this when the cameras are far apart is also simple:

Stereo requires you to match features from one image to the other. Hence the cameras must both have a view of the objects you are interested in.  Parallel cameras only overlap when the object is sufficiently far away - by pointing them inwards, you can solve this. (Of course, if your object looks significantly different from the two viewpoints, you still won't be able to do a good stereo match.)


Group delay plot of a Moving average filter I am implementing the moving average filter in Matlab. Using fvtool on this filter, gave the following output curve, is this group delay curve correct?

Edit:

My confusion is regarding the values of the delay output as well as the constant nature of the curve.

I implemented the filter using 'FilterDesigner' tool in matlab
(In which, filter dimensions can be specified in the GUI, which then handles the code itself) Attaching a screenshot of the tool window with filter parameters.

After which I exported this filter as an object called 'Hd'

To plot the filter, I did -
fvtool(Hd)

Which resulted in the curve above (figure 1)

You could have added to your question why you doubt the correctness of those numbers. Anyway, the group delay of a moving average is indeed constant. It's just the parameter 'window size' that is off by one. The window size $N$ is the number of taps of the moving average filter, and the group delay is related to the window size by
$$\tau_g=\frac{N-1}{2}$$
However, in your plot, the group delay is half of the window size, so I expect that what is called 'window size' actually equals $N-1$, i.e., the number of taps minus $1$.


What's the point of using Morphological processing concerning images? I dig deeper in internet to figure out what's the point of morphological processing on images, but i'd rather to hear the answer from people with experience, can you give me a hand with this? Thank you in advance

Many operations can be performed on images for enhancement and restoration. Sometimes, they seem mere handcraft and tinkering.

To make more sense of such workflows, and provide it with solid roots, it is useful to consider that  images are classes of abstract objects, and operations as well-defined structural actions on the former. Depending on premises or assumptions on the image shape (is it binary, continuous-valued, what do we expect from filtering, what elements are meaningful, etc.), several constructions can be derived.

Standard linear filtering techniques in images assume that images are composed of linear superpositions of components. This yields convolution-based methods, and treating images as elements of a linear vector space (the underlying structure), with a strong "continuous" underlying theory. This can be related to Fourier or harmonic analysis, such a strong influence of least-squares method.

But image formation is often non-linear, for instance because of occlusion or saturation. So you have a collection of non-linear tools, with homomorphic, min, max or median filters for instance, or stack filters for bit-wise versions. Plus,  they can be discretized in both values and space.

Mathematical morphology is another comprehensive framework, based on another underlying mathematical structure, lattices/set theory and related operations (erosion, dilation), with a strong discrete background (especially on image values).

Both have been proved recently connected in some ways, through the Cramer transform, see for instance An Explanation for the Logarithmic Connection between Linear and Morphological System Theory. Another interesting reading is 1989's Fourier analysis, mathematical morphology, and vision by C. Ronse.


Normalized autocorrelation of a sum of two signals I obtained the autocorrelation function of a sum of two signals:

as

but I want to know what is the normalized autocorrelation of above equation?

That would be a scaled autocorrelation such that its maximum value, hopefully attained at $\tau=0$, has value $1$.


Is document image binarization a closed research field A lot of work has been done in the field of document image binarization. Is it a closed research field or still an interesting challenge area? I am new in the field, and just curious.

First, in science, a field is rarely closed, sometimes asleep only. Resistance to low-contrast, real-time, badly  scanned, composite documents/writers or from aging medium seem to remain challenges, in similarity with other digital data: robustness, speed, poor acquisition, source separation or error concealment (the same topics, with generic names) are common for other data source.

As many data fields, this domain is being pervaded by (deep) learning techniques threatening the old handcraft (eg Learning Document Image Binarization from Data, Yue Wu et al.,  2015).

DIBCO 2017, the ICDAR 2017 Document Image Binarization COmpetition, is still on, a good sign altogether, and a source of challenges.7

As a personal note, I have been very interested in the analysis of old written documents, and their authentification, based on:

binarization from alternative light sources (multispectral, X-ray, ultraviolet), especially to detect previously scrapped, washed off  or erased writings, like from palimpsests (Binarization of MultiSpectral Document Images, 2015)
accurate analysis of non-written elements like the texture of paper fabric


ASK and Manchester Encoding Firstly, the context: I've been given a signal in a wav file, the sampling rate and carrier frequency along with the fact it is Manchester encoded and has been ASK Modulated. Five periods of a sinusoidal wave at half frequency of the carrier denote the end of a 'row' and a parity column/row is added. From this I am to reproduce a 300x300 black and white image. My main tool I've been given (never used before) is MatLab. This IS an assignment, and I'm just looking for some answers that I can't find clearly elsewhere.

So, from what I understand I've got this y[t], which is equal to x[t]*cos(2*pi*1000*t). From this signal I need to try and get x[t] out, then decode it. One of my issues is how to go about extracting this signal.

Sorry for the long start, and if any of this seems basic. It's an introductory course that just throws us in the deep-end all of a sudden!

Question 1: As this is ASK after Manchester Encoding, I'm guessing each 'bit' will be represented by a group of values from the sample that go from either positive to negative or negative to positive? Is my guess correct?

Question 2: In the samples, there are exactly five (or six, I've forgotten at present) zeros between each group of four non-zero numbers (except at the parts which I'm guessing are the five periods denoting the end of a row). Would these zeroes be useful in the demodulation, or are they just introduced during the sampling phase and not useful for demodulation?

Question 3: So, assuming the above questions are correct. To demodulate this signal, would I simply remove all the zero values, then iterate through and denote a '1' when it is positive and a '0' when it is negative (maybe some additional value for when the end of a row is encountered), then pair it up so 10 -> 1 and 01 -> 0?

Many thanks for any help, and sorry once more about the length!

I would suggest you first demodulate your signal. Take the received signal y[t] = x[t]*cos(2*pi*1000*t) and multiply it by cos(2*pi*1000*t). Remember, $ cos(A)cos(B) = 1/2 (cos(A-B) + cos(A+B)) $. You'll be left with x[t]*0.5*(cos(0) + cos(2*pi*2000*t) = x[t]*0.5*1 + x[t]*0.5*cos(2*pi*2000*t). Filtering the remaining signal with a low pass filter at 1k Hz, leaving only x[t]. From there you can try and answer the rest of your questions.

In Matlab this will look like:

t = 0:1/Fs:END_TIME;  % Fs is your sampling rate. This makes a vector of individual time ticks
y_demod = y.*cos(2*pi*1000*t);   % This multiplies every element in Y with the cos(...)
                                 % ... element. Make sure they're the same size.
a = (1/Fs)*(Fs/4);               % A cutoff frequency of approximately 1/2 the sampling rate
x_est = filter([1 1-a],a,y_demod);

The simple 1-pole low pass filter was taken from here: https://stackoverflow.com/questions/1783633/matlab-apply-a-low-pass-or-high-pass-filter-to-an-array


Better way to mute fundamental baseband frequency? I'm working on a project to monitor the fluctuations and harmonics present on standard 120VAC mains power. I built an adapter that drops the voltage from 120VAC to 12V and then down to 0.5-0.6 VAC nominal and feeds that into sound card input to convert it to digital form where I can analyze it with gnuradio.

The dominant frequency that is present of course is the fundamental mains power frequency (60Hz here in the US). I want to mute this frequency so I Can study all the other attributes of the signal.

I came up with a way to mute the fundamental 60Hz signal but I don't think it's the best way. My method is simple - it takes the current signal and adds the inverse of the same point on the signal to it. Picture a sine wave and a cosine wave plotted together and when added, they always add up to 0.

The problem is that I suspect this method may be adding artifacts to the output. I suspect so because I created another flow graph with sound card input and output sink and performed the same "add signal to its delayed inverse self" method I mentioned above. What I hear is a very distorted version of the actual input which includes an echo of the input signal.

Aside from the artifacts it does an amazing job of muting the 60Hz fundamental. So maybe it is the best way, but I want to crowd-source this challenge and ask the community for ideas on how to mute the 60Hz fundamental while adding little or no artifacts to the output signal.

Fig. 1: Flow graph and "zoomed out" scope focused on the current signal (blue) and the inverse delayed signal (green):

Fig. 2: Zoomed in scope showing the calculated output (red line). Note that it does appear to show the precise fluctuations of the input signals while muting the fundamental as expected, I'm just suspicious that it's also a somewhat unrealistic representation due to my method of calculation.

EDIT 1:
Thank you @dan-boschen for the notch filter suggestion. To get a notch filter in my gnuradio flow graph, I used a FIR filter with a band_reject taps function (firdes.band_reject(1,samp_rate,56,64,2)). Playing with the start/stop/transition parameters doesn't change much. The fundamental is significantly attenuated but is still present and much stronger than the very minute signal deviations that are present, and can be heard using the other method (both methods are compared in the video through left and right audio channels)

I just uploaded this video to illustrate the gnuradio flow graph, and left and right speaker channel comparison between FIR notch filter along, and FIR filter plus the "add to inverse last cycle" method I came up with.

I suspect the band_reject notch filter is tuned as best it can be. Perhaps there is another level of filtering I can apply. Maybe a series of band_reject filters? Looking for any ideas.

Use a notch filter such as the one shown in the figure below, where $\omega_n =2\pi 60/f_s$, where $f_s$ is your sampling rate, and $\alpha$ is chosen based on the bandwidth of the notch and how long you can allow for settling in the time domain; the tighter the filter bandwidth the longer it will take to settle; you can use a first order approximation of 10%/90% rise time = $\frac{0.35}{BW}$ where the BW is in Hz (and the amount of precision you can allow in the filter is an important factor---insufficient precision and $\alpha$ to close to 1 can lead to instability).

See this post for details: Transfer function of second order notch filter

Note the diagram below to facilitate implementation, mapping the coefficients to the transfer function as commonly referenced in Matlab and Octave. Do note that the signal flow below is going from right to left, which is opposite the diagram above (this was because my slide was showing how the implementation came from transposing the Direct Form II implementation). To be clear, the resulting coefficients are:

$b_0 = 1$

$b_1 = -2cos(\omega_n)$

$b_2 = 1$

$a_1 = -2\alpha cos(\omega_n)$

$a_2 = -\alpha^2$

A normalizing gain of $(1+\alpha)/2$ is added to the whole filter. Since $\alpha$ is typically close to 1, this can reasonably be omitted in most cases.


BIBO stability and input at system's pole My understanding of BIBO stability is that if a bounded input is applied to a discrete-time, linear, time-invariant BIBO system, the output will also be bounded. A sufficient condition for this is that the poles of the transfer function are within the unit circle.

Suppose I have a system with the transfer function
$$
H(z)={(1−2z^{−1})(1−{\frac{1}{6}} z^{−1}) \over (1−\frac{1}{2}z^{−1})(1−\frac{1}{3}z^{−1})}
$$
The poles of this system are within the unit circle and so is a BIBO stable system.

If I apply an input of $x[n]=(\frac{1}{3})^n$, which is bounded, the output would involve $H(\frac{1}{3})=\infty$, which is unbounded - i.e. a bounded input produced an unbounded output.

Obviously, there is something amiss in my understanding here, but what?

The signal $x[n] = (\frac{1}{3})^n$ is not bounded. As $n \to -\infty$, $x[n] \to \infty$; therefore, the system's response to your candidate signal $x[n]$ has no bearing on its BIBO stability.

A bounded, causal version of the signal (and likely what you were thinking) is $x_b[n] = (\frac{1}{3})^n u[n]$, where $u[n]$ is the unit step function. The response of the system to this signal is bounded as you would expect.


Is there any new scientific capability to build non-causal filters in real world? According to this post:

  In discrete-time systems, causality is a requirement only when
  processing (filtering) signals in real time; i.e. when the
  index nn corresponds to a physical time n×Tsn×Ts. In this case, a
  non-causal system is impossible to implement since calculating the
  current output would require inputs from a future time.

  In practice, however, many systems are not real time.

  Image processing, for example, has indexes for horizontal and vertical
  coordinates. So it's OK to apply a non causal filter to an image
  (usual examples being 0-centered FIR filters for blurring or
  sharpening, or IIR filters applied in both directions to obtain 0
  phase distorsion).

  In time sequences, many times the signal is stored and then processed
  in non real time, so you may apply non-causal filters there also.

I like to know is there any new scientific theory which could yield to making none causal FIR Filters in the real world. Like quantum theory in physics or other scientific fields?

Thanks for your attention.

is there any new scientific theory which could yield to making none causal FIR Filters in the real world?

No. If anything, scientific theories seem to be challenging determinism itself, which is what could allow you to peer ahead in the future. The key here is the last phrase of that quote ("...many times the signal is stored and then processed...").

If you already had acces to the true $x[n+1]$ at the level of observation of a signal, why not use a causal filter anyway? What is so special about $x[n+1]$?

For instance, take a very simple finite difference filter. You can definitely express it as:

$$ y[n] = x[n+1] - x[n]$$

But you would not be able to apply it in real time. The objective here is to take the difference between two samples. If you wait for $x[n+1]$ to happen, you can simply:

$$ y[n] = x[n] - x[n-1]$$

If it was possible for this difference to work in the presence of $x[n+k]$, then $y$ could start producing some output before $x$ even begins to occur.

You can observe exactly what this would sound like in a time reversed reverb. Here is an example of that.

To try this yourself, load a sound in Audacity, apply "Reverse", "Reverb" and finally "Reverse" again, before hitting "Play".

In reverberation, present samples bear some resemblance to past samples. If you stand at the far end of a large room, sounds from sources with a direct line-of-sight to your ears reach you much faster than sounds that have to bounce off of the walls. The sounds that bounce off of the walls travel a longer distance and therefore reach you later. This is what is causing the time difference of different reflections in reverberation.

By reversing the sound clip, you make the current sample to bear resemblance to its future values, thus starting to hear reflections of an event that has not yet happen.

This is only possible with signals that are stored.

Hope this helps.


Recommendations for SINAD equipment I am interested in buying a device that measures the SINAD of the audio signal(FM signal) that I receive from PC.
Here is the setup:

Signal generator (generating 1kHz FM signal ) $\rightarrow$ attenuator $\rightarrow$ HackRF connected to PC with SDR#

When we receive a 1kHz FM signal generated by the signal generator at the HackRF (which is connected to PC), we get an audio tone from the PC as output. I am interested in measuring the audio quality of FM tone that we receive at hackrf (connected to PC with SDR# support).

I know that for some signal generators, there is a separate port for measuring SINAD. The signal generator that I have does not have the capability to measure the SINAD. So, if there are any other devices that can help me measure this SINAD, please let me know.

Note: Sorry, I asked the same question without registering in dsp.stackexchanege (Any recommendations for buying SINAD measuring equipment?). I was unable to edit the question. So created a new one here.

So, you've got the audio signal in digital on your PC.

Analyze that digitally; you don't need any measurement equipment to do that.

On the contrary: Measuring properties of the analog audio signal is highly undesirable, because you'd be including the sound card of your laptop in your measurement, but you probably don't care about that even remotely. If you actually wanted to know the properties of your sound card, then you'd need to separately measure that – but you'd do that with a "clean" test signal generated by your PC, not something received over a noise system.


Different type of padding in image? I want to implement a gradient operator on an image in Matlab. Should I have to pad the image before its implementation? How do I decide which padding to use e.g. whether to use zero-padding or Symmetric padding? How are they differ from each other? I am new to this field of image processing, Can you provide some references to learn this stuff?

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

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

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

A couple of pointers:

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

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

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

Image padding, REFLECTION padding


Fourier Synthesis With Varying Frequencies Over Time I'm doing some fourier synthesis by making bins of complex numbers, filling them out with the amplitude and phase information and using IFFT to turn that into time domain samples.

I'm able to make sine, saw, triangle and square waves, and other more exotic periodic wave forms.

My question is, how do you do fourier synth for sounds that have frequency amplitudes that change over time?

Do you do it in small windows and append them together? If so, how do you make sure the phase's line up for each frequency from window to window?

Thanks so much for any info!

You might want to have a look at the overlap-add (http://en.wikipedia.org/wiki/Overlap%E2%80%93add_method) and overlap-save (http://en.wikipedia.org/wiki/Overlap%E2%80%93save_method) methods. But, if all you are trying to do is additive synthesis, you don't absolutely have to use the IFFT to generate your signal.

You can set up a bank of $M \in {[1 \ldots N_{Osc}]}$ elementary oscilators of the form $x_m(n)=\alpha_m \times sin(\frac{2 \times \pi \times f_m \times n}{Fs}+\phi_m) $ where $m \in M, N_{Osc} \in \mathbb{N}$, $\alpha_m, \phi_m, F_s, f_m \in \mathbb{R}$ and final output $y(n) = \sum_{m \in M}{x_m(n)}$.

In this case, $\alpha_m, \phi_m$ are amplitude and phase coefficients for each oscillator at frequency $f_m$ and to make them variable in time, all you have to do is turn them to $\alpha_m(n), \phi_m(n)$, i.e. to make them dependent on time sample $n$. This can be an absolute definition, by actually defining a varying waveform or by using some form of interpolation to provide initial and final values for these quantities and let the intermediate values be calculated automatically.

This would work both for generating one sample at a time or a frame of samples at a time (by iteration) and would not have any problems with continuity.


Practical Applications of upsampling and downsampling I'm a complete novice to digital signal processing, and I've started reading a few books on it.  In particular, I've run across the concept of upsampling and downsampling. As I understand it, given a discrete time signal $x(n)$, upsampling transforms the signal by (essentially) inserting $L-1$ zeros between each sample.  Similarly, one can downsample by omiting ever $L-1$ samples.  I'm not sure I understand why one would want to do so.  Why insert 'data' into the signal that wasn't there originally (upsampling)?  Wouldn't removing data (downsampling) remove features that are essential to its analysis?  I'm curious to know under what circumstances is this particularly useful?

Simplest reason for downsampling is simply to reduce the amount of data you have. Say you have a audio raw stream sampled at 44.1 kHz. You may simply want to reduce your data rate to 22.05 kHz because your algorithm can work with that sampling rate. Apply the same principle anywhere else. Just because you have more data does not mean you always need it. It may sound untrue but it is like saying do I need something upto 10 decimal points of accuracy or will 2 decimal points do? [Not exactly the same but you get the idea]

A simple reason to upsample to to interpolate. You have data which is sampled at a particular rate. You want to upsample it to get intermediate points between the actual data.  [Think of curve fitting. you want more points on the curve.]

There are implications to downsampling and upsampling. Not every downsampling is the same and not every upsampling give accurate interpolation. You will learn more as you go along.

Edit: Thought I'd give another simpler real life example
Image scaling: if you have a small image. You upscale it. This is upsampling.
If you want a smaller version you resize to a smaller size. This is downsampling.

Ofcourse practical scalers will not just add 0's and upsample and downsample. [You could do that too, but the output will not look good]. But it should give you another idea of where upsampling and downsampling is done.


How to calculate the phase shift AND time delay of non-periodic signals Say I have two non-periodic signals, $f_1(t)$ and $f_2(t)$, with Fourier transforms $F_1(\omega)$ and $F_2(\omega)$. Basically, I need to line up $f_1(t)$ and $f_2(t)$ as close as possible, and I am allowed to shift the signals in time and multiply their Fourier transforms by a constant phase.

What is the most efficient method to find the values of the phase shift $\Delta \phi$ and the time delay $\Delta t$ such that, when the Fourier transform of $f_2$ is modified into

$$ F_2(\omega) e^{i (\Delta \phi + \omega \Delta t)}, $$

the modified $f_2$ best "lines up" with $f_1$ according to their cross-correlation (or some other measurement)?

I am interested in simple and computationally efficient methods, even if they are not necessarily perfect.

A pure time delay could be determined by looking for a peak in the cross correlation. But in your case $f2$ might also have an overall phase offset.

You could try to compute two cross correlations:

$$
\begin{align}
x &= cross(f1,f2) \\
y &= cross(f1,hilbert(f2)) \\
\end{align}
$$

where $hilbert(f2)$ refers to an overall 90° phase shifted version of $f2$. If you combine those two like this

$$ z = \sqrt{x^2 + y^2} $$

you should get something that is independent of the phase shift and shows you a peak at the correct time delay $\Delta t$. The "phase" at that peak, $atan2(y,x)$ should give you the phase offset $\Delta\phi$.

I don't know if such a problem is usually solved this way and I have not tried it myself. But it might work.


I can't figure out why spectral features are "displaced" on a spectrogram I'm trying to create a spectrogram of a complex signal 'y', but I'm running into a problem where the spectrogram that I create appears to have its spectral features "displaced" in frequency space, and I think that it has something to do with not having pre-processed the signal correctly.

I know the signal's:
sampling rate (Fs = 400 KHz)
bits/sample (Nbits = 16)
frequency that the signal is centered on (Fc = 5.60 MHz)
total time that the signal elapses (dT ~ 145 sec).

I tried to use matlab's built-in spectrogram function with the following arguments:
[S, F, T] = spectrogram(y, num_pts_win, num_pts_overlap, num_pts_fft, Fs);
where:
num_pts_win = num_pts_fft = 2^14
and:
num_pts_overlap = 1/2 num_pts_win = 2^13

After subtracting Fs/2 from F to get a frequency axis from -Fs/2 -> Fs/2 rather than from 0->Fs, I plotted the spectrogram of the log(S) and obtained:

while I was expecting to get something like this:

I believe that this displacement of the spectral features has resulted from not accounting for the central/carrier frequency (Fc) anywhere, but as I'm sort of new at this, I don't really understand exactly how I should account for it.

Edit 1:

To clarify, I suspect that the spectral features I'm seeing in the first graph are actually aliases of the true spectral features, and I think this aliasing is happening because I haven't accounted for the central/carrier frequency.

Edit 2:
After FFTSHIFT'ing over the frequency dimension I think I've cleaned things up nicely!

After subtracting Fs/2 from F to get a frequency axis from -Fs/2 -> Fs/2 rather than from 0->Fs

You can't center the frequency by changing the labels on the x-axis. Use 'fftshift' instead.


How to combine bins of my DFT I have a time series and apply the FFT to get a spectrum.

Let's assume that my sampling frequency and the length of the time sample are chosen such that I end up with a $\Delta f = 0.1$ Hz.

As this is still rather noisy, I'd like to change $\Delta f$ to 0.5 Hz a posteriori.

How do I combine the fourier components of the adjacent bins? Am I right to just sum them? Or do I take the square root of the sum of squares?
Is there a term for what I am trying to achieve, so that I can find some references on my own?

There are different methods for doing what you want.
But first i have to tell you power spectrum density (PSD) is real and positive and you have to obtain it from squared FFT of signal or from FFT of signal's auto-correlation (FFT of auto-correlations is always positive).
To reduce noise you could apply some damping window on your signal. This method is called modified periodogram.
Also you could divide your signal to smaller segments and find PSD of each segment then find average of PSD of each segment. This method is called Bartlet method.
Also your suggested method, summation of adjacent bins or applying moving average works too (you must apply it on PSD).


How exactly does one apply a Hanning Window for a spectral estimate? PREAMBLE

I am a graduate student researching evolutionary ecology. My supervisor and I have been trying to learn how to simulate colored noise, discretely, by the $\frac{1}{f^\alpha}$ power law for $\alpha \in \mathbb{R}$, specifically on the interval [0,2].

We found an article that seems to cover a method in great detail.

Link here: Discrete simulation of colored noise and stochastic processes and $1/f^\alpha$ power law noise generation (Kasdin, 1995)

On page 806, we were able to prove to ourselves that by a symmetric autocorrelation function, the formula for the sampled spectrum of an arbitrary segment of a process is indeed,

$$\hat{S}(\omega)\triangleq E\{\tilde{S}(\omega)\}=\int_{-T}^{0}\left(1+\dfrac{\tau}{T}\right) \Big[ \dfrac{1}{T} \int_{t_o+\tau/2}^{t_o+T+\tau/2}  R(t, \tau)\Big] e^{-j\omega\tau}d\tau + \int_{0}^{T}\left(1-\dfrac{\tau}{T}\right) \Big[ \dfrac{1}{T} \int_{t_o+\tau/2}^{t_o+T+\tau/2}R(t, \tau) \Big] e^{-j\omega\tau}d\tau$$

$$=\int_{-T}^{T}\left( 1-\dfrac{|\tau|}{T}\right)R(t,\tau)e^{-j\omega\tau}d\tau$$

Where,

$$R(t,\tau)\triangleq \{x(t+\frac{\tau}{2})x(t-\frac{\tau}{2})$$

Next, the paper computed the spectral estimate for Brownian motion. That estimate in the paper is given as,

$$\hat{S}_B(\omega)=2(1+\dfrac{t_o}{T})\dfrac{1}{\omega^2} - 2(\dfrac{t_o}{T})\dfrac{cos\omega T}{\omega^2} - \dfrac{2}{T\omega^3}sin\omega T$$

Where for $T$ much larger than $t_o$, this reduces to $\frac{2}{\omega^2}$. We were able to reproduce this result.

QUESTION AND PROBLEM

The paper now makes the claim that we have unknowingly been using a rectangular window. We are told we can improve the spectral estimate by using a Hanning window with normalization $8/3$ (After some reading, I've been able to grasp the concept of spectral leakage).

The new spectral estimate for Brownian Motion, with $t_o$ set to 0 is,

$$\hat{S}_B=\dfrac{1}{3}\{  \dfrac{4}{\omega^2} - \dfrac{4sin \omega T}{T \omega^3} + \dfrac{80\pi^4 - 24\pi^2T^2\omega^2 + T^4\omega^4}{T^4(4\pi^2/T^2-\omega^2)^3)} + \dfrac{48\pi^2T\omega sin \omega T - 4T^3\omega^3sin\omega T}{T^4(4\pi^2/T^2-\omega^2)^3)}  \}$$

Try as we might, we have not yet been able to reproduce this result. We tried the following,

$$\int_{-T}^{0}\left(1+\dfrac{\tau}{T}\right) \Big[ \dfrac{1}{T} \int_{t_o+\tau/2}^{t_o+T+\tau/2}  R(t, \tau)\Big] e^{-j\omega\tau}d\tau \cdot HW(x=t)\cdot HW(x=(t+\tau)) + \int_{0}^{T}\left(1-\dfrac{\tau}{T}\right) \Big[ \dfrac{1}{T} \int_{t_o+\tau/2}^{t_o+T+\tau/2}R(t, \tau) \Big] e^{-j\omega\tau}d\tau \cdot HW(x=t)\cdot HW(x=(t+\tau))$$

Where we defined the Hanning Window (HW) as,

$$1-cos((\pi x)/T)^2$$

Done with the following code in Sage:

Q, T, t0, w, tau, t, x = var('Q T t0 w tau t x')
hanning_window = (sqrt(8/3))*(1 - cos((pi*x)/T)^2)
rectangular_window = 1
assume(T>0)

window_1 = hanning_window
window_2 = hanning_window

fneghann = Q*(t+tau/2)*exp(-I*w*tau)*window_1.substitute(x=(t))*window_2.substitute(x=(t+tau))
fposhann = Q*(t-tau/2)*exp(-I*w*tau)*window_1.substitute(x=(t))*window_2.substitute(x=(t+tau))

result = (1/T)*((fneghann.integrate(t,t0-tau/2,t0+T+tau/2)).integrate(tau,-T,0) + (fposhann.integrate(t,t0+tau/2,t0+T-tau/2)).integrate(tau,0,T))
f = (result).expand().full_simplify()
view(f)

Computing the limit as $T\rightarrow \infty$ and setting $t_o = 0$ DOES produce the spectral estimate of $\frac{1}{\omega^2}$ as indicated in the linked paper, but we got something dissimilar to the initial spectral estimate submitted in the article, leading us to believe the result of our limit was luck. I've chosen to not post the result, $f$, because of it's size.

Where did we make a mistake, if any? Did we apply the Hanning Window incorrectly? Is there a mistake in our code? Your help is greatly appreciated.

Winston and I believe that the observation window is from 0 to T in two variables, $t$ and $\beta$.  The integral equation Winston cites is after a transformation to $\tau = \beta-t$ (see page 147 in "Random Data" by Bendat and Piersol).  One of the things which makes us suspicious is that the result of the integral in Sage has terms with $i$ in them, which cancelled in the simpler, rectangular window case.  Thus we believe we are doing something wrong.


Choosing the right overlap for a window function On Wikipedia there is quite a list of window functions.  However, when you are doing a lot of STFTs, you often have to have some overlap to get better results.

Generally an overlap factor of 2 (i.e., 50% overlap) works fine (I'm doing analysis & resynthesis), but for some windows it isn't good (I currently have this experience with the Blackman window).

I tried my best but I couldn't find a resource that would list the "good" overlap factors for common and less common windows.

Please refer to this paper on the optimum overlap percentage for the Blackman-Harris window, which is derived to be 66.1%. It has a lot of useful information on spectral analysis and windows.


Discrete time sine wave generation near nyquist I'm trying to wrap my head around how to generate sine-waves out of a DAC near the Nyquist frequency (or determining how close I can get for reliable results). So if I want to generate a 499 Hz sine wave using a 1 kSPS DAC, I will calculate each sample:
$$y(t) = \sin(2 \pi 499 \ t) $$
Converting this to discrete time:
$$ y[k] = y\left( \tfrac{k}{1000} \right) = \sin \left(2 \pi 499 \frac{k}{1000} \right)$$
I end up with something like this:

The red line shows the 499 Hz signal, while the blue dots show the DAC sample. I get these large amplitude fluctuations. I know it's common to use a reconstruction filter to get rid of DAC images, which I always thought of as a low-pass filter from DC to Nyquist, which wouldn't really help get rid of the amplitude fluctuations...or would it?
Is there a better approach for generating sine waves? Or is there just simply a limit I cannot exceed (and what is the limit)?

There's nothing wrong about this – it's absolutely what you should be seeing.
The sine wave is intact, if you applied a reconstruction filter – there's no amplitude fluctuation of the wave, as soon as you limit its bandwidth to the nyquist bandwidth, the analog signal "has" to follow the full sines.


Adding a constant to a bandpass-signal in transmission I do understand AM-DSB-SC and AM-DSB-TC, but can't see the point in adding a constant DC-offset of 2 in the upper part of this Simulink-model?

Source

You're right, from a pure information point of view, the constant just wastes transmit power and doesn't add anything to the signal that you'd care about – the receiver will filter it out.

However, you've got to realize that all AM audio transmission schemes are positively ancient. So, receiver's abilities to demodulate a signal were very limited.

Simple (analog) AM receivers are nothing but a frequency selective circuit (e.g. a tuned antenna) and a rectifier, and a low-pass filter – done. If you're off a bit when tuning the circuit, it doesn't really matter.

Now, in your stereo AM (never seen that, wouldn't do that, seems strange and wasteful in spectrum) transmitter,  it's very important to tune to the right frequency for the first and second channel to be correctly received.

Now, that requires something a simple receiver could "lock" on. And that's exactly what the constant is: It adds a constant peak in the spectrum. If you have something like a PLL, you can lock on to that (make the PLL's bandwidth small, so that it's not affected by the AM).

Modern receivers could be much smarter – it's not very hard to guess the "center" of a double-sideband AM spectrum without a pilot tone.