Why does the standard waveform from audio out get distorted when I reduce the amplitude? In my attempt to generate a standard sine wave out of a sound card, I started using matlab which has a function soundsc for scaling the output

 t = 0:1/48000:10;
 y = sin(1000*2*pi*t);
 soundsc(y,48000,[a b]);

It's ok when my a and b are -1 and 1 respectively. When I try to make them even -0.9 and 0.9 respectively, the standard sine wave of 1Khz gets distorted, why?

Are our sound cards not able to handle standard sounds?

Primarily because the documentation is poorly written. If you want to play at half amplitude you need to use

soundsc(y,48000,[-2 2]);

The range argument is the input amplitude range that gets mapped to the maximum available range of the sound card. So if the max of your input exceeds the max in the range argument, you will get clipping.


How to compare multiple polynomials in MATLAB? I am trying to fit a second order polynomial to various curves on a binary image using MATLAB polyfit function.There are many such curves on the image, and the coordinates of the pixels forming the curve are known.

How do I compare the co-efficients, so as to arrive at one approximation to all such curves?

Here are some co-efficients of such curves:

 0.0010126     -1.5981       1074.8
 0.0024871     -3.3327       1866.3
 0.002662      -3.7067       1631.4
 0.0025574     -4.5129       2148.4
 0.002114      -3.3457       2195.6
 0.0015383     -2.6023       1682.9

Based on what you share, this is the best answer I can give.
I see two methods:

Method #1,

Instead of making N polynomials (if N is the number of ribs), you could fit a curve for all the ribs at once. Of course, for that, you will need to have a common (0,0) coordinate, in other word, stack them one over the other on your image. Since you are able to make a distinction between the ribs, that means that you can already isolate them, so that shouldn't be an issue.

Or if you prefer, you can generate a new cloud of points with the polynomials you already found.

coeffs = [0.0010126, -1.5981, 1074.8; 0.0024871, -3.3327, 1866.3];
x = linspace(0,1000,1000); % fictives points
merged_x = [x x];
merged_y = [ polyval(coeffs(1,:), x), polyval(coeffs(2,:),x)];
coeffs_final = polyfit(merged_x, merged_y, 2);
figure
hold on
plot(polyval(coeffs(1,:), x), 'red')
plot(polyval(coeffs(2,:), x), 'blue')
plot(polyval(coeffs_final, x), 'black')

Method #2

You could fit 3 polynomials to find the values of the coefficient you already have found. If we take the first polynomial you provided, then
$a_2=0.0010126$ $a_1=-1.5981$ $a_0=1074.8$

You could fit a curve that would follow all $a_2$, another for all $a_1$, and same for all $a_0$. This would yield 9 new coefficients (3 polynomials) that will allow you to find your polynomials based on the rib position.

Here's an example

coeffs = [ 0.0010126     -1.5981       1074.8;
    0.0024871     -3.3327       1866.3;
 0.002662      -3.7067       1631.4;
 0.0025574     -4.5129       2148.4;
 0.002114      -3.3457       2195.6;
 0.0015383     -2.6023       1682.9 ];

x = linspace(0,1000,1000); %fictive value
ribs_count = size(coeffs, 1);
a2 = polyfit(1:ribs_count, coeffs(:,1)', 2); % careful to the ' symbol that transpose the column into a row
a1 = polyfit(1:ribs_count, coeffs(:,2)', 2);
a0 = polyfit(1:ribs_count, coeffs(:,3)', 2);

figure
for rib=1:ribs_count
    subplot(6,1,rib);
    hold on
    new_coeffs = [polyval(a2,rib), polyval(a1,rib), polyval(a0,rib)];
    plot(x, polyval(new_coeffs, x), 'red')
    plot(x, polyval(coeffs(rib,:), x), 'blue')
end


Difference Between Iteratively Reweighted Least Squares (IRLS) and Sequential Quadratic Programming? Part of my work is concerned with applications in Sparse Bayesian Learning and therefore I occasionally stumble over interesting papers in the field of compressed sensing.

I recently read Iteratively Reweighted Algorithms for Compressive Sensing by Rick Chartrand and Wotao Yin (Available also on DocDroid).

The paper describes how using $ {L}_{p} $ -norms with $ p < 1 $ can be used to recover signals with fewer measurements than with the LASSO ($ {L}_{1} $ Regularization).
There is even a Wikipedia entry on Iteratively Reweighted Least Squares (IRLS).

However, I can't wrap my head around the difference between IRLS and Sequential Quadratic Programming (SQP). Is there any difference? The papers I have found on IRLS never mention SQP..

Many thanks in advance!

SQP is a method for solving smooth (objective and constraint functions are at least twice differentiable) constrained nonlinear optimization problems.  It solves a series of quadratic programming problems to converge to a solution to the Karush-Kuhn-Tucker conditions for the constrained optimization problem.

IRLS is a method for solving unconstrained minimization problems by solving a sequence of least squares problems which are obtained from the original problem by a rescaling at each iteration.


Narrowband Signal Decomposition Why do the in-phase and quadrature components of a narrowband signal each have half the bandwidth of the signal?

The in-phase and quadrature components of a bandpass signal are real-valued baseband signals. So if the spectrum of a bandpass signal with bandwidth $B$ extends from $f_c-B/2$ to $f_c+B/2$ (where $f_c>B/2$ is the carrier frequency), the spectra of the in-phase and quadrature components extend from $-B/2$ to $B/2$. Since bandwidth is measured at positive frequencies, the baseband bandwidth is $B/2$.


Generate signals with a particular variance and SNR Consider a system model of the form: $y_n = ax_n + v_n$ where $x_n$ is the input that is corrupted by $v_n$ which is an Additive White Gaussian Noise of zero-mean and variance 1 for $n = 1,2,...,N$ samples. $a$ is an unknown parameter representing the channel coefficient. How do I generate $v_n$ which is White but of a different variance other than $1$ and a particular SNR?

QUESTION: I am not sure whether the following way creates a signal of non-unit variance and how to determine the variance? If I create a signal z_10 using awgn() with a particular SNR say 10 dB, then would its variance be different from another signal, z_20 created using SNR = 20 dB? What is the proper way to create signal of a particular variance and know its SNR?
y_wnoise = y + sqrt(variance)*randn(size(y))  but how do I include the SNR value?

The way I have implemented in MATLAB is as follows. Did I do it correctly?

%generate data:
N = 50; %number of data points
s = randn(1,N);
a = 0.6;
for n = 1:N
   y(n) = a*s(n);
end
SNR = [10,15]

%generate noisy signal of different variance
z_10 = awgn(y,10,'measured');
z_15 = awgn(y,15,'measured');
OR
z1 = y + sqrt(0.6)*randn(1,N);

It depends on what you mean by SNR. It's a common joke in the DSP community to spell it out as "something to noise ratio", referring to the fact that there is no unique definition of SNR, so the term by itself means nothing.

Define it yourself and use it appropriately.

What's common is to define it as ${\rm SNR} = \frac{P_{\rm s}}{P_{\rm n}}$ where $P_{\rm s}$ is the power (variance) of the signal samples ($x_n$ in your notation) and $P_{\rm n}$ is the power (variance) of the noise samples.

Hence:

P_s = 1;   % target signal power
SNR = 15;  % target SNR in dB
P_n = P_s / 10^(SNR/10); % calculated noise power

s = randn(1,N)*sqrt(P_s);
v = randn(1,N)*sqrt(P_n);
y = a*s + v;

You don't need the for loop there by the way. You can use the awgn function, it should do something very similar. Again, there is no clear need for it though.

Note that this definition of the SNR is independent of a. Depending on your application this may be desirable or undesirable. In the latter case, you can include it in your SNR definition. Feel free to do so. Just define everything very clearly.

*edit: In response to your question: the SNR is really a ratio of signal power to noise power. So to change the SNR, you can change both. If you want the SNR to be 10 dB higher, you can decrease the noise variance by 10 or alternatively increase the signal variance by 10.

If you want to change the SNR while keeping the noise variance fixed (say, 1), you're welcome to do so. Almost the same approach:

P_n = 1;   % target noise power
SNR = 15;  % target SNR in dB
P_s = P_n / 10^(-SNR/10); % calculated signal power

s = randn(1,N)*sqrt(P_s);
v = randn(1,N)*sqrt(P_n);
y = a*s + v;

This also means looking at the noise variance itself doesn't tell you about the SNR conditions. If you assume it's 1 this can mean a low SNR (if the signal power is, say, below 1) or a high SNR (if it is much higher than 1).


Do Causal Discrete-time systems have proper transfer functions? In the case of continuous-time systems, if the system is causal, its Laplace transfer function is strictly proper (the degree of the numerator is less than the degree of the denominator).

Is this true also in case of discrete-time systems? That is, if the system is causal, is its Z transfer function strictly proper?

Note that a stable and causal continuous-time transfer function does not need to be strictly proper but only proper, i.e. the degree of the numerator does not exceed the degree of the denominator, but numerator and denominator degree can be equal. E.g.

$$H(s)=\frac{as^2+bs+c}{s^2+ds+e}$$

can represent a causal and stable system, as long as its poles are in the left half of the complex $s$-plane.

For discrete-time systems the same is true. A transfer function for which the degree of the numerator is greater than the degree of the denominator has at least one pole at infinity. Since causal and stable discrete-time systems must have all their poles inside the unit circle, such a system can't be causal and stable. However, as is the case with continuous-time systems, equal degrees of numerator and denominator are possible.

Assume a causal and stable discrete-time system is described by the following linear difference equation with constant coefficients:

$$y[n]+a_1y[n-1]+\ldots+a_Ny[n-N]=b_0x[n]+b_1x[n-1]+\ldots+b_Mx[n-M]\tag{1}$$

The corresponding transfer function is

$$H(z)=\frac{\sum_{m=0}^{M}b_mz^{-m}}{1+\sum_{n=1}^Na_nz^{-n}}\tag{2}$$

Note that there are of course no restrictions on the choices of $M$ and $N$. Nevertheless, the transfer function $H(z)$ is always proper. Assuming $M>N$ gives

$$H(z)=\frac{b_0z^M+b_1z^{M-1}+\ldots+b_M}{z^M+a_1z^{M-1}+\ldots+a_Nz^{M-N}}\tag{3}$$

Obviously, the transfer function given by (3) is proper. $M>N$ just means that the last few coefficients of the denominator polynomial are zero. A very similar argument holds for $M<N$.


Kalman Filter Process Noise - Model Where It Vanished I am trying to use the Kalman filter (the scalar version) to estimate the steady state of a set measurements which is a random process.
I have used a constant dynamic model as the state equation,

$$ x(t+1) = x(t) + w(t) $$

Where $ w(t) \sim N(0, Q) $.

The measurement equation

$$ y(t) = x(t) + v(t) $$

Where $ v(t) \sim N (0, R) $.

The measurement noise $ v(t) $ is taken as the uncertainty of the instrument and R is a known value with a assumed Gaussain distribution.
The process error, $w(t) = x(t+1) - x(t)$ has high values as a constant variance (Calculated using the total error variance) or a time varying (Taking last 2, 5 readings as the set f sample).
However, the state reaches the steady state when $ Q $ is very small or  $ Q = 0 $.

I have seen several other theoretical examples uses $ Q = 0 $ and reaches the steady state. But I am working with actual data. Also I have seen some other posts explaining the similar behavior. Can some explain this behavior with reasons, please?

The process noise $w(t)$, which is typically assumed to be a zero mean, white Gaussian noise with (power) variance $\sigma^2$, is used to account for any mismatches between the assumed dynamic model of the states and the actual truth.

If your assumed constant state (DC) model perfectly matches with the underlying nature of the observed signal, then you should set the process noise to zero; i.e., your model is exact.

Then you only have measurement noise due to the sensor's intrinsic noise.

If you do not set your $\sigma^2$ to zero (adding non-zero process noise to a constant dynamic model) then your state estimation will never reach steady-state, but wander around it; as you are telling the Kalman filter that the true signal is not a constant...


Play IQ data in the soundcard I have short-wave digital receiver that send RP packets with IQ components over network depend of the RF frequency and passband. I need to work with the telephone passband (0.3 - 3.4 kHz). The task is:

Using ordinary notebook receive RTP packets in the real-time, extract IQ-data, build FFT, visualize it and then I need to playback this data on the ordinary soundcard.

Please, help me out. How can I do this quickly (in the short time)? Maybe prepared solution exists.

If the spectrum of your I/Q samples is centered at zero then you'll have to perform either AM or FM demodulation before routing any real-valued audio samples to a sound card. For AM demodulation you'll need to implement a complex-input "envelope detector" which produces a real-valued audio signal riding on a DC bias. (In a few days check the web page:

www.dsprelated.com

for a blog of mine describing different envelope detectors.) That amplitude-biased audio signal then needs to be passed through a "DC blocking" filter (a highpass filter whose cutoff frequency is very close to zero Hz) and the output of that filter goes to your sound card.  From Chapter 13 of my DSP book, the most common DC blocking filters look like the following:

where alpha is a number close to, but less than, one.

Now if your I/Q signal (whose spectrum is centered at zero Hz) is an FM modulated signal there are two common ways to perform FM demodulation to produce a real-valued audio signal sequence.  One way is unpleasant because it requires the computation of arctangents and the other way is unpleasant because it requires you to perform divisions. Oleg, search the Internet for "digital FM demodulation" and you will be able to learn about those two different FM demodulation methods.


Why are the Continuous Wavelet Transforms of the same signal drastically different? I'm currently studying wavelets and am running into confusion with regards to CWT coefficients. Ideally, I want a CWT algorithm that produces outputs similar to that of a STFT - i.e. produces coefficients that correspond to frequency magnitudes of the input signal. However, I'm finding conflicting results.
From the CWT Wikipedia page, I'm given the illustration below. This illustration shows that the wavelet coefficients are essentially just sampling different subbands of the signal. This is not what I'm looking for.

However, when I look at a MATLAB web page discussing the CWT and DWT, I find a different illustration of the CWT when operating on the same input signal. This diagram looks more similar to what I was hoping for where each frequency has a constant coefficient.

What is the difference between CWTs?

There are differences in axis-scales. However, I suspect that the main difference is that the first picture relates to either a real CWT (or time-frequency transform) or the real-part of a complex CWT (or time-frequency transform, again), while the second seems to be a magnitude.
In the first  image, the apparently wiggling part could "planified" by combining it to some imaginary part. Then changing the $y$-scale from linear to logarithmic could make the plots look closer.
Apart from representation issues (scales, magnitudes), with admissible continuous wavelets, there exists a wide variety thereof.  Factors that affect scalograms in practice are also the global shape of wavelets (symmetry and anti-symmetry for instance), their regularity, and the ways it is discretized (number and sampling of voices).


Audio Bit Depth "mapping" from DAW to Audio card: is this correct? Let say I'm generating a 1 kHz sine wave within my DAW, using some synth (Sytrus in my case). I normalize it, so it play at max "digital" amplitude 0db.

Thus, any sample of the signal would gets a Floating Point values (in 32 or 80 bit of precision, using FL Studio 32bit; but that's irrelevant right now).

This means that all values (since the signal is normalized) goes from -1.0 to 1.0 max.

Now, the signal go through ASIO Drive and reach my sound card (an M-Audio FireWire Solo), 24-bit (fixed-point representation). Due to ENOB, it will use 21-bit, thus the max value would be 1048576 (which is 100000000000000000000 in binary).

Is correct to state that the value 1.0 is mapped to 100000000000000000000? And -1.0 to 000000000000000000000?

I'm really not sure about this. Does this mapping happen? Is it correct the range I've used (-1.0/1.0 => 000000000000000000000/100000000000000000000)?

Can you help me to understand this "mapping", if any?

If your sound card has input range 2Vp-p (2 volts peak to peak = -1V to +1v) and has 24 bits of resolution (ENOB doesn't play any role in this, you still get 24 bits for each sample for each channel, but the higher ENOB the smaller noise you get), here is what voltage represents the output code of ADC:

00000000 00000000 00000000 (lowest possible result from ADC) => -1V
10000000 00000000 00000000 (midpoint) => 0V
11111111 11111111 11111111 (highest possible value) => +1V

Or with 3 bit ADC:

000 (0 dec) => -1V
001 (1 dec)
010 (2 dec)
011 (3 dec)
100 (4 dec) => (midpoint approximately) 0V
101 (5 dec)
110 (6 dec)
111 (7 dec) => +1V

So the voltage is proportionally distributed in the full ADC range, to calculate the voltage correspondig to specific code, you can use following formula:

U = (code/(2^bits-1)) * (Umax - Umin) + Umin

where Umin = -1V, Umax = 1V, bits = 24
For conversion between binary and decimal numbers refer to this document


Interpolate DFT Coefficient of a Frequency That Is Not in the DFT Bin I'm using the Jacobsen interpolation to get a more precise frequency of my signal. To get the corresponding DFT coefficient I'm doing:

$$X_{f} = \sum_{n=0}^{N}{x_{n}  e^{-2\pi ifn}}$$

Where $x_{n}$ is the samples of my signal and $f$ is the frequency found by the interpolation (that is not one of those in the DFT bins).
This is a quite heavy computation for large signals.

Is there a method to interpolate the coefficient corresponding to that frequency without having to compute the sum?

If you need very specific frequency you should either use the Discrete Fourier Transform (DFT) definition or the Goertzel Algorithm.

Yet, if you need higher resolution Grid, you can use many off the shelf algorithms to interpolate.
The native Interpolation Kernel would be the Dirichlet Kernel since it is the Sampling Kernel for the DFT.
An easy way to interpolate the grid in the Frequency Domain using the Dirichlet Kernel is by Zero Padding of the signal in the Time Domain.
As this is equivalent of higher rate sampling of the Discrete Time Fourier Transform (DTFT).


How to calculate the energy levels of frequency intervals of an audio signal? Using the Fourier transformation I can translate a signal from the time-domain to the frequency-domain. If I divide the frequency axis into intervals, how can I calculate the audios signals energy for each of these intervals?

You may research for spectrogram. The above comments from pichennets and B Z are correct, but you have to pay attention to the effect of just take some subset of coefficients and calculate the power. When you take this way you're applying a retangular window in frequency, and some peaks in the frequency domain power. You must apply a window in time domain before take the FFT and calculate the power.

I want to add this as a comment, since it's not a real answer but I'm not abble to do this.


Getting a more accurate frequency read from autocorrelation and peak-detection algorithm For a project I am attempting to create an automatic tuner for a guitar, which reads the audio from the guitar jack, determines the frequency and adjusts the string by a motor.

Using http://www.instructables.com/id/Reliable-Frequency-Detection-Using-DSP-Techniques/ the algorithm is quite accurate.

void setup() {
      sum = 0;
      pd_state = 0;
      int period = 0;
      // Autocorrelation
      for(i=0; i < samples; i++)
      {
        sum_old = sum;
        sum = 0;
        for(k=0; k < samples-i; k++) sum += (values[k])*(values[k+i]);

    // Peak Detect State Machine
    if (pd_state == 2 && (sum-sum_old) <=0)
    {
      period = i;
      pd_state = 3;
    }
    if (pd_state == 1 && (sum > thresh) && (sum-sum_old) > 0) pd_state = 2;
    if (!i) {
      thresh = sum * 0.5;
      pd_state = 1;
    }
  }
  // Frequency identified in Hz
  freq_per = realFreq/period;
  lcd.cls();
  lcd.printf("%f",freq_per);
  }

At the end the frequency is determined by the by dividing the realFreq, which is the precise reading frequency with the period determined by the peak detection portion.

The frequency is quite accurate between 40-200 Hz normally about 0.5 Hz out which is accurate enough for my needs. But above 200 Hz the accuracy falls to 2-3Hz.

Is there a technique to improve the accuracy of this function. The sampling frequency is 4410 and I wanted to reduce this to improve the speed but keep the accuracy as high as possible.

Thank you.

If you were using an omnidirectional microphone to pick up the instrument sound you wanted to determine the pitch for, then using the autocorrelation method of determining the period of a waveform could have some merit in reducing the impact of random uncorrelated noise.

But, the way it is attempted to be used here seems to be an incomplete match to the problem's specifications. The signal is not expected to have high levels of noise, the hardware platform is "weak" (limited wordlength and processing power) and all that the selected method is doing is a "fancy" way of deriving a (periodic) signal's period. In addition to this, at some point later on, the same platform is going to have to control the motors too.

Since the signal will be coming directly from the guitar's pickup, it will be loud enough and clean (assuming no problems with the instrument's wiring). In this case, the zero-crossing method is probably the better choice for a guitar tuner. All that the processor would be doing is count how much time passes between setting and resetting a counter as the signal crosses a "zero voltage" point within limits. This method would return the pitch of a guitar string with much higher accuracy given the clock speeds that an Arduino is capable of running at and also leave plenty of time to the CPU to perform other tasks with ease (such as controlling the motors).

For more information on pitch detection in DSP, please see this link.

For more information on zero crossing, please see this link (linked from the same Instructables article that was posted along with the question).

For more information on the basic technique required to create the circuit that will be converting the guitar's signal to the digital form that Arduino expects, please see this link.


doesn't OFDM inherently deals with ISI? why cyclic prefix is needed? As far as I can understand, the motivation to OFDM and multi carrier communication in general was alot due to the ISI caused by multi-path which is equivalent to frequency selective channel.
Now, by using OFDM, we send our symbols over sub-carriers, thus each symbol in the frequency domain is affected by approximately a constant channel response.
Then, why is there an aditional need for cyclic prefix? doesn't using sub-carriers enough for dealing with the ISI?

OFDM (sub)carriers are only exactly orthogonal over a certain IFFT block length (and integer multiples).  Multipath echos will mix a subcarrier signal from the end of one symbol block into the beginning of the next, rendering the symbol in that carrier to be ambiguous over the initial one IFFT block length.  A Cyclic Prefix extends the length of each symbol beyond just that of one IFFT length so that the IFFT can be (hopefully) done over a portion of the signal past where all the orthogonal subcarriers have been messed up with multipath echoes.  Skipping past the messed up portion allows the (constant or otherwise) channel response to be mostly ignored by the receiver, while maintaining the orthogonality of all the subcarriers from each other.


Zero-order interpolation problem Let $x_c(t)=\cos(\omega_0t)$. This signal is sampled with $\omega_s$, which is greater than the Nyquist rate. It is then interpolated with a zero-order interpolator. The signal obtained is $y_c(t)$. What should be the relationship between $\omega_0$ and $\omega_s$ so that $y_c(t)$ is periodic? Also, find the relationship between the amplitudes of the first two non-zero Fourier coefficients of $y_c(t)$.

I guess that, for $y_c(t)$ being periodic, $Y(j\omega )$ should consist of just deltas. The problem is that I don't see how this could happen. Mathematically,

$$X_p(j\omega)=\frac{\omega_s}{2}\sum\limits_{k=-\infty}^\infty [\delta(\omega -k\omega_s- \omega_0)+\delta(\omega -k\omega_s+ \omega_0)]$$
$$H_0(j\omega )=\frac{2\pi}{\omega_s}\text{sinc}\left(\frac{\omega}{\omega_s}\right)e^{-j\frac{\pi\omega}{\omega_s}} $$

With $x_p(t)$ being the sampled signal and $h_0(t)$ the impulse response of the interpolator. So

$$Y_c(j\omega)=X_p(j\omega)\cdot H_0(j\omega )$$

You're thinking too complicated.

For a digital signal to be periodic, it simply has to repeat every $P\in\mathbb N$ samples.

Now, you have a continuous signal which is $\frac{2\pi}{\omega_0}$-periodic.

Obviously, if a period happens to be exactly "worth" a natural number $N\in\mathbb N$ sample times, the resulting sample from sampling it will repeat with a period of $N$.

More generally, if $M\in\mathbb N$ periods of the continuous signal are as long as $N$ sampling times, the digital signal will also be periodic.

So, as a general rule, one period of the input signal must be worth $\frac{N}{M}$ sampling times, or:

$$\begin{align*}
\frac{2\pi}{\omega_0} &\overset != \frac NM \frac {1}{f_\text{sampl}}&\quad N,M\in \mathbb N\\
&= R\frac {1}{f_\text{sampl}} &\quad R\in \mathbb Q\text,\\
\text{which means that}\\
\frac{f_\text{sample}}{\omega_0} &= \frac {R}{2\pi}\text;\\
\text{with}\\
f_{sample}&=\frac{\omega_s}{2\pi}\\
\text{it follows that}\\
\frac{\omega_s}{\omega_0} &= R\text,
\end{align*}$$

i.e. the relationship between signal frequency and sampling rate needs to be rational.

The fact that you have a zero-order interpolator afterwards doesn't change this consideration -- the interpolator's output will be periodic iff its input is.


Transmit diversity vs multi-path? Are transmit diversity signal and the multi-path signal same? Multi-path contains LOS+delayed and attenuated/distorted versions of it and transmit diversity signal contains LOS1+delayed LOS2+ +delayed and attenuated/distorted versions of the two.

Kindly, elucidate in terms of the impact of each type of signal on the equalization used at the receiver.

Thanks, JK

Multi-path is a property of the communication channel (mostly the wireless channel). The received signal $r(t)$ is the sum of delayed and weighted versions of the transmit signal $x(t)$:
$$
r(t) = \sum_{i} a_i x(t-\tau_i)
$$

Transmit diversity is a property of the transmitter. It describes a technique using multiple transmit antennas, thus creating multiple communication channels that are different because the antennas are placed at different locations in space. The idea is to transmit the same information at the same time at different antennas.

For example, let's assume one receive antenna and two transmit antennas (2x1 multiple-input single-output (MISO)) the received signal is a superposition of both multi-path channels:
$$
r(t) = \sum_{i} a_i x(t-\tau_i)+\sum_{j} b_j x(t-\tau_j)
$$
The hope is that if transmission on one channel fails you will still receive a useful signal through the other channel thereby increasing the robustness of the system.

There is also the concept of receive diversity where multiple receive antennas are used (SIMO). And of course you can use multiple antennas on both receiver and transmitter side (MIMO).

When using the scheme of transmitting the same data at the same time, an FIR filter can be used no matter if transmit diversity is applied or not, because to the receiver it will just appear like an ordinary multi-path channel. Advanced equalizer structures include the fractionally spaced equalizer and the decision-feedback equalizer. If you apply some space-time coding a the transmitter an according space-time decoder is required at the receiver.


The Scale Space Theory Understanding In the scale-space theory the scale-space representation of the signal $f(x), x = (x_1, ..., x_d)$, (in case of image $d = 2$) is given as:
$L(x, y; t) = g(x, y; t) * f(x, y)$ where $g(x, y; t)$ is a gaussian kernel with parameter $t$ and $*$ is a convolution. By changing the $t$ parameter we receive a more or less smoothed image. As the result coarser representation (parameter $t$) will not contain small objects or noise.

The main point is to find a way of scale-invariant feature detection, right? So that for some image at it's reduced in size copy the features like keypoints will be detected correctly, even if size is different, without finding other noise keypoints.

In the paper they are using the $\gamma$-normalized derivatives.
$\delta_{\xi, \gamma-norm} = t^{\gamma / 2} \delta_x$.
What is the meaning of using the $\gamma$-normalized derivative, how does it helps in scale-invariancy?
From this image we can see that at near the same positions the different keypoints found (different in size). How is that possible?

If you can explain the step-by-step algorithm of scale-invariant feature detection, this would be great. What is actually done? The derivatives can be taken by $x, y$ or $t$. Blob can be detected by taking the the derivative of $L$ by $(x, y)$ variables. How is the derivative by $t$ is helping here?

The paper I was reading is: Feature detection with automatic scale selection

It really has been a long time since I have read Lindeberg's papers, so the
notation looks a bit strange.  As a result, my initial answer was wrong.
$\gamma$ is not a scale level.  It seems to be a parameter of some sort that can be tuned.  It is true that you need to multiply the derivative by the
appropriate power of $t$. $t$ itself corresponds to a scale level, and the power  depends on the order of the derivative.
You can find keypoints at multiple scales in the same location.  That is because you look for the local maxima over scales.  Here's the intuition: think of an image of a face.  At a fine scale you get a blob corresponding to the nose.  At a course scale you get a blob corresponding to the whole face.  The two blobs are centered at the same point, but have different scales.
Here is the whole algorithm:

Decide which image features you are interested in (e. g. blobs, corners, edges)
Define a corresponding "detector function" in terms of derivatives, e. g. a Laplacian for blobs.
Compute derivatives that you need for your detector function at a range of scales.
Multiply the derivative responses by $t^{m \gamma / 2}$, where $m$ is the order of the derivative, to compensate for magnitude decrease.
Compute the detector function over the whole scale space.
Find local maxima of the detector function over $x, y, t$.
These are your interest points, or keypoints.

Edit:

Lindeberg proves in the paper that $t^{\gamma / 2}$ is the appropriate factor for normalization of derivatives.  I don't think I can reproduce the proof here.
You do not take derivatives with respect to $t$.  You only compute derivatives with respect to $x$ and $y$, but you compute them at a range of scales. One way to think about this is you generate a Gaussian scale space first, by repeatedly blurring the image with a Gaussian filter of some variance $t$.  Then compute derivatives with respect to $x$ and $y$ at each scale level.
You want to find local maxima over scales because you may have image features of different size at the same location.  Think of an image of concentric circles, like a bulls-eye.  It will give you high responses of a Laplacian at several scales.  Or think of an image of a real human eye filtered by a Laplacian at a range of scales.  You will get a high response at a fine scale for the pupil, a high response an some medium scale for the iris, and a high response at a coarse scale for the whole eye.

The whole point is that you do not know at what scale the features of interest might be ahead of time.  So you look at all scales.


What are the properties of continuous-time improper systems? I am trying to better understand the properties of improper systems $H(s) = \frac{b(s)}{a(s)}$, for which the order of the numerator $b(s)$ is greater than the order of the denominator $a(s)$ (in the broader context of continuous-time linear time-invariant systems). I have trouble reconciling the application of region of convergence (ROC) rules and my intuition.

Let us take the simplest example that is the pure differentiator, i.e.,  $H(s)=s$.

Application of ROC rules

Causality. The system has a pole at $s=\infty$. Thus, the ROC is not a right-sided plane. So the system is not causal.
BIBO stability. The ROC includes the imaginary axis $s=j\omega$. The system is BIBO stable.

My intuition

Causality. The differentiator can be expressed as the limit of a forward difference or as the limit of a backward difference. So depending on how we define the differentiator, it can be a causal system or an anti-causal system.
BIBO stability. I have two intuitions here.

(1) The derivative of a bounded input like the step function can be unbounded as it gives the Dirac function. So the system is not BIBO stable.

(2) The impulse response of the differentiator is $h(t) = \delta'(t)$, i.e., the derivative of the Dirac function. This is not an absolutely integrable function as $\int_{-\infty}^{\infty} |\delta'(t)| dt = \infty$. See
https://en.wikipedia.org/wiki/Unit_doublet. So the system is not BIBO stable.

How can I reconcile the ROC analysis and my intuition? Where am I wrong?

Definitions

Causality. A LTI system with impulse response $h(t)$ is said to be causal, that is, it has the property that the value of the output at time $t_0$ depends on the values of the input and output for all $t$ up to time $t_0$ but no further, i.e., only for $t \leq t_0$, if and only if $h(t)=0$ for all $t < 0$.

BIBO stability. A LTI system with impulse response $h(t)$ is said to be BIBO stable, that is, it has the property that its output $y$ is bounded whenever its input $u$ is bounded, if and only if
$$\int_{-\infty}^{\infty} h(t) dt \triangleq \| h \|_1 < \infty $$
in which case we can bound the peak of the output by
$$ \| y  \|_\infty \leq \| h  \|_1 \| u  \|_\infty $$
where $\|u\|_\infty \triangleq \mathrm{sup}_{t\in\mathbb{R}} u(t)$.

Note that the imaginary axis is not inside the ROC because $|H(j\omega)|\to\infty$ for $\omega\to\infty$. So the system is definitely unstable, regardless of causality, which confirms your intuition.

Concerning causality, the ROC is neither a right half-plane nor a left half-plane (both because of the pole at $s\to\infty$), which is consistent with a (generalized) function concentrated at $t=0$.


How to calculate the bandwidth of a continuous signal? Can someone please explain how to compute the bandwidth of a signal? Like why the bandwidth of $5\sin(2t)$ is $0$ and of $\sin(2t) + \sin(3t)$ is $1$?
Thanks.

I suppose the correct expressions you meant are:

a) $5\sin(2t)$

b) $\sin(2t) + \sin(3t)$

Bandwidth is defined as the difference between the highest and lowest frequencies of a given signal ou system.

With this in mind, signal a) has one single frequency of 2 rad/s and so its bandwidth is 2-2=0 rad/s.
Similarly, signal b) has 2 frequencies: 2 rad/s and 3 rad/s. So, bandwidth is 3-2=1 rad/s.


$f \star (g \,* \,h) = (f \star g) * h$? In this paper that I'm reading, I see this equation:
$$S_0 \star (S_0 * N) = N,$$
where $\star$ denotes cross-correlation, $*$ denotes convolution and $S_0 \star S_0 = \delta$. Is it true in general that $f \star (g*h) = (f \star g)*h$ and if not, why is it true here?

It's easiest to look at this in the frequency domain. We have
$$\mathcal{F} (x*y) = X \cdot Y$$
$$\mathcal{F} (x \star y) = X' \cdot Y$$
where ' is the complex conjugate
So we get
$$\mathcal{F} (f \star (g * h)) = F' \cdot G \cdot H $$
$$\mathcal{F} ((f \star g) * h) = (F \cdot G)' \cdot H = F' \cdot G' \cdot H  $$
These are only identical if G is real.
In your case you let's simplify the notation to $f = S_0, g = S_0, h = H$ $S_0$ has a unit impulse as autocorrelation, which means that the spectrum of the autocorrelation is just unity, i.e. $1$
$$\mathcal{F} (f \star f) = \mathcal{F} (\delta) = 1 = F' \cdot F$$
We use that property and can show that
$$\mathcal{F} (f \star (f * h)) = F' \cdot F \cdot H = H$$
since $F' \cdot F = 1$


AR model order selection for half second EEG fragments I am using MATLAB to evaluate power spectral density estimates of half second EEG signals, using modified covariance method. Can anyone suggest me how to select the AR model order for this process? Is there any function in MATLAB which can be used for selecting the best order using any model order selection techniques?

Thanks a lot for your help.

Do you want to select the model order adaptively based on the content of the input EEG, or choose the AR order as a design parameter for a system that will always be the same?. In the second case (not adaptive), maybe if you have some sample EEGs you could try a a lot of models from order = 0, to some number N (large enough but not so that running all the models takes forever). And Graph the Error between the EEG signal and white noise filtered by the model. You should arribe at a graph or Error vs order. At some value of N the error only marginally decreases. That is your model order.

EDIT

Matlab offers AIC and others.


Why is `fftfilt` (i.e. `fft` of both inputs, then element-wise multiplication, then `ifft`) faster than direct convolution? I have the following matlab code

N = 500;
M = 32;
K = N + M - 1; % Length of full convolution

x = randn(N);
y = ones(M);

% first one: convolution
r = conv(x, y, 'same');

% second one: multiplication
z = ifft(fft(x, K, K).*fft(y, K, K));
idx = round((K-N)/2+1):round((K+N)/2); % Indices of the central portion
z = z(idx); % Obtaining the central part of the result

May I know which opertation is faster? getting $r$ or $z$? Why is this so? Thanks in advance

Whether the direct convolution or the FFT/IFFT method is faster depends on the length of the impulse response, $N_\mathrm{i}$ and the signal length $N_\mathrm{s}$. With the formulas taken from here I've created a small Matlab script that calculates the required number of real multiplications and additions for the direct convolution and FFT/IFFT method, respectively. I've assumed real-valued impulse response and signal. Taking your example of $N_\mathrm{i}=32$ (or the other way round, this makes no difference), the direct convolution always requires less multiplications:

But if we increase $N_\mathrm{i}$ to 50, for example, the direct convolution requires less multiplications up to a signal length of about 60. For signal length greater than 60 the FFT/IFFT method requires less multiplications:

The actual computation time depends on some other parameters, especially the CPU architecture but the number of multplications is usually a good indicator.

The MATLAB code for the above figures for reference:

Ns = 1:1000; % length of signal
Ni = 50; % length of impulse response

K = Ni + Ns - 1;

% number of real multiplications for convolution:
M_R_conv = Ns*Ni;
% number of real multiplications for convolution:
A_R_conv = Ns*(Ni-1);

% number of complex multiplications for freq domain conv:
M_C_dft = 3/4 * K .* (log2(K) + 1) + K;
% number of complex additions for freq domain conv:
A_C_dft = 3/2 * K .* log2(K);

% number of real multiplications for freq domain conv:
M_R_dft = 4*M_C_dft;
% number of real additions for freq domain conv:
A_C_dft = 2 * A_C_dft + 2 * M_R_dft;

figure('Position', [0 0 400 300]);
plot(M_R_conv, 'b');
hold on;
plot(M_R_dft, 'r');
hold off;
legend('direct', 'fft', 'location', 'Southeast');
title(['number of real multiplications. Ni: ' num2str(Ni)]);
xlabel('length of signal');

A sidenote to your code: as it stands it is not running as you're creating NxN and MxM matrices for x and y. What you probably want to do is

x = rand(1,N);
y = ones(1, M)
% ...
z = ifft(fft(x, K).*fft(y, K));

p.s. for those of you interested in Python code instead of Matlab, here is the same code in Python: https://gitlab.com/snippets/1789085


Equivalence of Maximum Likelihood (ML) and Discrete Fourier Transfrom (DFT) Peak Finding for Single Tone Estimation My understanding is Maximum Likelihood and DFT Peak Finding for a single tone produce the same results assuming the ML is restricted to the same frequencies as the DFT.

I was wondering if there was an easy way to prove this?
My understanding is that ML will be the same doing Least Squares regression so I suppose an equivalent statement would be to prove that the DFT peak gives the best Least Squares estimate of single tone sinusoidal.

Maximum Likelihood under the assumption of Additive White Gaussian Noise (AWGN) is always equivalent to finding the hypothesis with the minimum distance to given data.

Since minimizing distance is equivalent (In the euclidean Space) of maximizing the correlation you can always build the idea of Match Filter for parameter estimation in the settings of ML under AWGN.

Matched filter means that you create a set of hypothesis parameterized on a grid and goes one by one until you find the one which maximizes correlation (Minimizes distance) to data.

This is exactly what's the DFT is doing.
The grid is based on the different tones (Frequencies) and the correlation level is given by the bin value.


How to filter noise from continuous path in image Problem statement

In the dataset above, you can notice a continuous path and then a noise-band at around 320-400 on the x-axis. I want to be able to apply a filter that amplifies the continuous path while suppressing the noise-band. The data inside the noise-band is useless.

Attempt 1

I am new to signal processing, but I took a shot at designing my own filter. I tried designing a minimum filter. I defined it to make the output signal $ S $ the minimum of the next $ w $ terms in the input signal $ F $: $ S[n] = \min \{F[x] | x \in \{n, n+1, n+2, ..., n+w\}\} $

for i in range(len(original_signal) - window_len):
    filtered_signal[i:i+window_len] = min(original_signal[i:i+window_len])

It represents each pixel by the lowest pixel in the area which should make the noise go away because the noise band will have a lot of low values in between a lot of high values. It takes out too much of the signal as well.

This is an image showing what the filter took away (the difference between the original signal and the filtered signal)

Attempt 2

A median filter will not work, because the median of the noise-band is relatively high (1e4):

whereas the median of empty space is not that high (1e3):

So the median filter would not reduce the noise to be even with empty space.

Do you know any filters I could try?

Additional Data

Histogram of desired signal:

All histograms normalized and put on the same axes:

So entropy filter sounded good, but it didn't really work because the actual signal is also considered high entropy, even after adjusting the window.

What I did was take a derivative of all of the data in the time axis (left to right). Then I capped the maximum value. Then I integrated the derivative.

deriv = data[i, 1:] - data[i, :-1]
deriv[deriv > max_] = max_

The derivitave is $ \frac{dy}{dx }$. In my mind dy is a vertical line segment, dx is a horizontal line segment. Putting them together makes a right triangle, where the function is the hypotenuse is f, the horizontal leg is dx and the vertical leg is dy. Thus $ tan^{-1} {\frac{dy}{dx}} $ is the angle of elevation for f. In my mind, putting a maximum on the derivative is like putting a maximum angle on f. It is like saying, "Cap the angle of elevation for f. If it goes above 150 degrees from the horizon, make it smaller."

It seems to work. It limits the maximum angle of elevation for the signal. You can see the original signal (in red) and the signal after derivative filter (in green).

Then I take original signal minus derivative filter applied to original signal. This leaves me with just the noise band, which I can amplify, blur, and then subtract out.


Find response of discrete time LTI system given input and impulse response For this question the guide says to use $Y(e^{j\omega}) = H(e^{j\omega}) X(e^{j\omega})$.
I have been able to find the discrete time Fourier transform of the impulse function as $$H(e^{j\omega}) = \frac{1}{1-\frac12 e^{-j\omega}}$$
Finding the same for the impulse function is what I do not understand.

Note that the given definition of $x[n]$ is undefined for $n=0$. But even if you changed it by replacing the unit step $u[n]$ by the shifted step $u[n-1]$ to exclude the value $n=0$, then the problem would still be very hard to solve. So it's more than likely that the correct definition of $x[n]$ has no $n$ in the denominator, and the current form is just a consequence of a typo or of a misinterpretation.


How MATLAB calculates matched filter gain Here is a bunch of code from MATLAB documentation: Designing a Basic Monostatic Pulse Radar:

wav = phased.RectangularWaveform(...
    'PulseWidth',1/pulseBandwidth,...
    'PRF',PRF,...
    'SampleRate',Fs);
matchingcoeff = getMatchedFilter(wav);

hmf = phased.MatchedFilter(...
    'Coefficients',matchingcoeff,...
    'GainOutputPort',true);

[rx_pulses, mfgain] = step(hmf,rx_pulses);

How MATLAB calculates the mfgain (matched filter gain)? Any formula for that?

Matched filter gain is calculated as (in dB):

$G_{dB} = 10 \cdot \log_{10}(L)$,

where $L$ - is filter length. It's maximum possible SNR improvement that filter can provide.

Try this formula and compare result with MATLAB one. I suppose it will be the same.


DFT from Fourier transform I always studied the DFT starting from his formula, but for some reasons I need to do comparison between the FT and the DFT. I found the pdf in this link very useful http://www.robots.ox.ac.uk/~sjrob/Teaching/SP/l7.pdf because it explains how the DFT can be obtained starting from the FT.
After the first formula it says " We could regard each sample f[k] as an impulse having area f[k]". What I don't understand is how it possible that the sample, being an impulse, has an area f[k]. Shouldn't t have an area f[k] multiplied by a very low number? This passage for me is crucial because I cannot understand why, passing from the FT to the DFT, the time increment disappears.

The area under an impulse function is one. So if you multiply a continuous time waveform with an impulse, the result will be a weighted impulse with the area equal to whatever the waveform was at that moment in time.


Calculating phase response of maximum phase filter using Hilbert Transform Given only a magnitude response $A(\omega)$ of a minimum phase filter, one can calculate the phase response using the Hilbert Transform:

$$θ(ω) = -\mathcal{H}\{\ln(A(\omega)\}$$

This paper suggests that there is another Hilbert Transform relationship between magnitude and phase for maximum phase filters, but I can't find anything beyond that in the cited references and elsewhere.

Does anybody know what the exact Hilbert Transform relationship is between magnitude and phase for maximum phase filters?

Note that a Hilbert transform relationship for maximum phase systems can only exist if a maximum phase system is defined to be an anti-causal stable system with an anti-causal stable inverse (cf. these MIT course notes, p.4). Considering only systems with rational transfer functions, this means that not only the system's zeros but also its poles must be in the right half-plane (for continuous-time systems) or outside the unit circle (for discrete-time systems), respectively. The reason for this is that for the Hilbert transform relationships to hold, $\ln(H(s))$ (or $\ln(H(z))$) must be analytic, either for $\text{Re}\{s\}\ge 0$ (or $|z|\ge 1$) (minimum-phase), or for $\text{Re}\{s\}\le 0$ (or $|z|\le 1$) (maximum-phase). Since both poles and zeros of $H(s)$ cause singularities in $\ln(H(s))$ we need to make sure that for maximum phase systems there are neither poles nor zeros in the left half-plane $\text{Re}\{s\}\le 0$ (or inside the unit circle for discrete-time systems).

Note that the above definition of a maximum phase system is different from another common definition which says that a maximum phase system is a causal and stable system with all its zeros in the right half-plane (or outside the unit circle), but - since it's causal and stable - all its poles in the left half-plane (inside the unit circle).

According to the first definition above, a maximum phase system is obtained from a minimum phase system by time-reversal, which leaves the magnitude of the frequency response unchanged but inverts the sign of the phase. Consequently, for a maximum phase system according to above definition, we simply get the following Hilbert transform relationship between phase and magnitude:

$$\theta(\omega)=\mathcal{H}\{\ln(A(\omega))\}\tag{1}$$

Note that the sign inversion in the Hilbert transform relationship between magnitude and phase can also be seen from the fact that a minimum phase system has a causal cepstrum (from which the Hilbert transform relationship follows), and a maximum phase system has an anti-causal cepstrum, resulting in the same Hilbert transform relationship between magnitude and phase but with an inverted sign.


IIR filter design fitting step response I have access to the step response of a system and I want to find its poles and zeros without knowing the order of the system.

Consider an example of a step response of the system shown in the figure below.  Here is the link to the raw data if someone wants to try it (https://pastebin.com/BpT4d9RW).

I want to build an IIR filter (with arbitrary number of poles and zeros) that has the same step response (or pretty close to it).

I have tried using the prny() function in MATLAB for that purpose, and the closest result I could get was the following. I am actually not sure what I am doing wrong here and if there are any considerations for using the prony() function that I am not aware of.

I would appreciate if someone can help me on correct use of the function OR can come up with some other approach that works.

Please also note that the approach should be general enough to accommodate any type of step responses (over-damped or under-damped).

There are two problems here that need to be solved. The first is the delay of the measured step response, i.e. how many samples from the beginning you need to discard. This is critical in your example because you chose a second-order system to model the data, and such a system is not very flexible in adding or subtracting delay. Delay could be added by adding more zeros to the system, i.e. by increasing the order of the numerator. If you don't want that then you will need to optimize the delay of the step response such that it can be optimally approximated by a second-order system.

The other problem is the fact that the designed filter doesn't reach the final value of the measured step response. The reason is that the impulse response required by Prony's method is derived by computing the difference of adjacent step response samples without taking into account that the step response was truncated. Since several samples of the measured step response are discarded, you also discard the same number of samples of the impulse response, and if you were to compute a step response from that truncated impulse response (by computing the cumulative sum), you would get a step response that is different from the measured one (and that would be close to the one implemented by the designed filter). The solution to this problem is very straightforward. Just truncate the step response and add a value of zero at the beginning. Then compute the corresponding impulse response by taking differences, and the first value of the impulse response will equal the first (non-zero) value of the step response, i.e., it will represent the accumulated past of the original impulse response (before truncation).

I did that and found as optimal starting value for the step response an index of $530$ for a filter order of $2$. The figure below shows that the step response of the designed filter closely approximates the measured step response:


Plotting Hilbert Spectrum from Hilbert Huang Transform Before anyone calls it a repeat post from past posts. I'd like to point out I've read all of those; however, I have questions on making a hard coded version of the hilbert spectrum.

I do indeed have all the ingredients aka the hilbert transform and empirical mode decompositions.

Here is my code below:

def hilbert_huang(self, imf, res):
    # Reconstruct the signal
    sig = np.sum(imf, axis=0) + res
    # Hilbert Transform of IMFs
    hil_imf = np.array([sp.fftpack.hilbert(imf[i]) for i in range(len(imf))])
    # Instantaneous Amplitude
    inst_amp = np.array([np.sqrt(imf[i]**2 + hil_imf[i]**2) for i in range(len(imf))])
    # Instaneous Phase
    inst_phase = np.array([np.arctan(hil_imf[i]/imf[i]) for i in range(len(imf))])
    # Instaneous Temporal Frequency
    inst_freq = np.array([(1/(2*np.pi))* np.gradient(inst_phase[i]) for i in range(len(imf))])

    return inst_amp, inst_phase, inst_freq, sig

However, I'm not sure how to construct the spectrum from all this, I know that the formula is given by
$$ H(\omega, t) = \Re\left(\sum_{i=1}^{n} a_{i}(t) e^{j \int \omega_{i} (t) dt}\right) $$

But I'm not sure how to fix $\omega$ , given that $\omega$ is an array with values , in the wikipedia article it says to let $\omega = \omega_{j}$ where $j$ indicates which IMF it's from but again $\omega$ is a vector, so do I just find the max? Let me know if I'm not clear or need to provide any more mathematical expressions or code.

The formula for the Hilbert spectrum in your question does not make much sense. Its right-hand side it is not a function of frequency, instead it is simply the sum of the intrinsic mode functions (IMFs), i.e., it represents the original signal. The Hilbert spectrum is the sum of the individual Hilbert spectra of each IMF:

$$H_j(\omega,t)=\begin{cases}a_j(t),&\omega=\omega_j(t)\\0,&\text{otherwise}\end{cases}$$

where $a_j(t)$ and $\omega_j(t)$ are the instantaneous amplitude and the instantaneous frequency (IF) of the $j^{\,th}$IMF, respectively.

Apart from that, note that your computation of the instantaneous phase and frequency will not work. For the phase you have to use arctan2 instead of arctan, and for computing the IF from the instantaneous phase, you need to unwrap the phase. However, it is generally better to compute the IF directly from the signal, as discussed in this answer.


How to combine a perfect signal with a limited dynamic range with a poor one with high dynamic range? I have two sensors that measure speed $v(t)$ of a moving vehicle. The first sensor produces a signal $f(t)$ which is a very accurate estimation of speed. However, it only works for slow to moderate speeds and beyond a certain limit, it simply fails to report anything. Let's assume we can model this behavior with an additional output of this sensor using signal $h(t)$, which is $1$ while $f(t)$ is valid and is $0$ when $f(t)$ is invalid.

Now, the second sensor produces signal $g(t)$ for all speeds but it is not particularly accurate. When $h(t)$ is 1, i.e. when speed is low and the first sensor produces perfect results via $f(t)$, the signal from the second sensor could be modeled as:

Equation 1: $g(t) = \beta_{0}(t) + \beta_{1}(t)v(t) + \epsilon$

The parameters for this linear relationship, $\beta_{0}$ and $\beta_{1}$ are unknown. All we know is that they vary very slowly over time.

So the problem is to estimate $v(t)$ given the 3 signals $f,h,g$.

One solution is of course to forget about $g(t)$ and just report the last value of $f(t)$ when it becomes invalid. Effectively, it means to saturate the output at a certain high speed.

Alternatively, my initial idea was to use a form of rolling linear regression to estimate the $\beta$ parameters with a fixed window size $w$: When $h(t)$ is $1$, I simply report $f(t)$ and update $\beta_{0}$ and $\beta_{1}$. When it is $0$, I report a $\hat{v}(t)$ estimated using $g(t)$ and the last calculated $\beta$ coefficients. This is obviously very slow since I have to solve a linear regression at each sample.

Questions:

(1) Does the idea of using rolling linear regression make sense? If yes, how can I do it fast?

(2) Would it be possible to something like a Kalman filter to solve this problem? I'm not sure how to express this problem in a linear state-space model.

(3) I guess the most challenging question would be what to do in the case that $f(t)$ is quantized (e.g. an integer value) while $v(t)$ and $g(t)$ are not. I guess we might be able to simulate it by assuming a uniform distribution for $\epsilon$ in equation 1.

In response to your question (1):

Thinking along the lines of your initial idea, there is a way to avoid regression and do the following approximate procedure instead:

Without loss of generality, assume that $\epsilon$ has zero mean (if not, then just absorb it into the constant $\beta_0$). Also assume it is i.i.d and has finite variance. Suppose you obtain two sets of data for $(f,g)$ pairs when $h=1$. I will denote these as $\{(f_{1,i}, g_{1,i})\}_{i=1}^N$ and $\{(f_{2,i}, g_{2,i})\}_{i=1}^N$.

Let $\bar{f_1} = \frac{1}{N} \sum_{i=1}^N f_{1,i}$ and similarly define $\bar{g}_1$, $\bar{f}_2$ and $\bar{g}_2$. If $N$ is sufficiently large, we can now claim that, at least approximately, the noise $\epsilon$ can be averaged out and we can write:
$$
\bar{g}_1 \approx \beta_0 + \beta_1 \bar{f_1}
$$
and
$$
\bar{g}_2 \approx \beta_0 + \beta_1 \bar{f_2}.
$$
Now we only need to solve two equations in two unknowns to get the scale parameter $\beta_1$ and offset $\beta_0$, provided $|\bar{f}_1 - \bar{f}_2|>0$.

Of course, whether this will work in your application will depend on whether the variance of $\epsilon$ is small enough (remember: averaging will reduce the variance by a factor of $1/N$), and whether the $f$'s are really noise free when $h=1$.

In response to question (3):

The same procedure should work in case $f$ is a quantized version of the true velocity. The argument here is that quantization introduces a small amount of noise that is approximately uniformly distributed with zero mean. It gets averaged out for large $N$, or so I hope.


Uncountable Set of Poles? It is easy to define an (ideal) LTI system that would have an infinite number of poles - for instance, if the transfer function is
$$
H(z)=\frac{1}{\cos(z)-1}
$$
However, this would only define a countably infinite set of poles.
I am curious: are there systems with uncountably infinite sets of poles?
I would guess that the very definition of a pole would prohibit a function from having a contiguous region of poles, but what about a hypothetical transfer function that would have its poles defined on something like Cantor set?

This really comes down to how you define LTI.
Far as I can tell, a time-discrete system being LTI means its in- and output are related by a countably large system of ordinary differential equations.
Since each of these transforms into a polynomial through the Laplace-, Fourier- or z-transform, the transfer function, in which ever domain it exists, is always a quotient of polynomials. Coefficients of polynomials are countable.
Since polynomials have, including their multiplicity, exactly the number of roots as their degree, and poles are just the roots of the denominator polynomial, we'll always have countably many roots – and thus, poles.

Not quite sure you'd want to relax the constraint that poles be countable, anyway: If you do that, transfer functions cease to be meromorphic, and that means the residue theorem doesn't apply, and as far as I can intuit, that leads to unsolvable integrals that you'd need to define the $z$-transform to begin with...


What does double hat mean in Spectral theory? I was confused about the double hat notation $\hat{\hat .}$ and my teacher and I thought it probably means an estimate of an estimate most of the time, but we could not be 100% sure of its meaning.
One hat $\hat .$ is an estimate but what does this double hat usually mean? Does it have many different meanings in DSP? Can it mean an estimate of unit vector sometimes? Does it have other signs whatever it means?
Example

You can see such double hat on page such as 213 in Prentice Hall - Modeling Of Dynamic Systems book. In the chapter on Periodograms: Blackman-Tukey's Method, they use this kind of spectral estimate for N length function:
$$\hat\phi_N(\omega)=\int_{-\pi}^{\pi}W_\gamma (\omega-\xi)\hat{\hat{\phi}}_N(\xi)d\xi$$

It seems this notation is used for empirical transfer function estimates (see pg. 207, 8.4 Fourier Analysis). To wit, on pg. 79 of Modeling: identification and simulation of dynamical systems we have

  The double hat denotes that the ETFE is a very rough estimate of the transfer function.


What are Local and Global Inpainting Techniques in Image Processing? Is Diffusion-based inpaiting Local or Global?
Is Pixel-based inpaiting Local or Global?
Is Patch-based inpaiting Local or Global?
Can Local-diffusion be used inside Patch-based in-painting problems of Image Processing?
Or, three of them can be applied separately to operate on the same type of problem?

Usually one could divide the methods according to:

Global Method
The function optimized is all over the image with no local data.
Namely it finds best solution involving all pixels.
Diffusion Based is usually iterative method to solve such a Global model.
Local Method
The optimization function is local for a group of patches and holes.
Most "Patch" based methods are local in their nature.

In practice, state of the art methods combine both approaches.
Namely global regulator over a local optimization function or the other way around.


Sharpening Using Laplacian Kernel - How to Handle the Output (Specifically Negative Values)? In order to sharpen an image using Laplacian kernel, it should be convolved with the original image.

It is evident that the result might have negative values. I am not sure what to do with these values:

Subtract the minimum from laplacian then map to [0, 255]
Clip negative values to zero and map to [0, 255]
Just clip negative values to zero

I know the next step is to add/subtract the laplacian to/from the original image and map to [0, 255] but I have doubts about the laplacian.

Indeed, I tested all three methods and 2, 3 got me to the right results but according to my Image Processing textbook, Laplacian should be scaled but it is not mentioned clearly that which version of laplacian is used to derive the sharpened image.

Thank you!

@SOFUser, Welcome to our community.

Since the High Pass Filter, such as the Laplacian Kernel, output has zero mean (As the DC Component is removed) it indeed can (Should) have negative values as the output oscillates around zero.

What to do with it depends on the objective:

Display Purpose
In case on wants to display the output of the Sharpening / HPF Filter process it is very common to add 50% gray to it. In case of 8 [Bit] image it means adding 127 / 128. This is what's don in Adobe Photoshop.
Processing Purpose
Usually we keep the negative values as they might be important in the next steps. For instance in Laplacian Pyramid if one won't keep the negative values the reconstruction won't be accurate. Another example would be sharpening where we add the output of the Laplacian Filter scaled by the sharpening factor to the image to have sharpened image.
Generating Image
If one wants to generate image and archive it in one of the UINT based formats (JPEG, PNG, etc...) usually one would clip data into the compatible range.


Parametric EQ vs multiple bandpass filters Will filtering a signal through multiple bandpass filters, where each filter is passed the complete signal, and then combining the output have the same result as a parametric EQ?

What if the signal is passed into a series of bandpass filters?

Are there any performance issues or optimizations that can be made when EQing a signal this way?

Update:

After reading the answer in this question which-filter-for-an-audio-equalizer it seems splitting the signal through a filter bank (in parallel) is the correct way.

The answer explains how the bands need to blend together to avoid audio issues, how is this calculated if each bandpass filter can have a varied bandwidth?

a graphic EQ need not mean cascaded shelves.  it could mean cascaded peak/cut bell filters.  maybe with shelves as bookends, maybe not.

the bell filters are symmetrical in log freq (or the pitch scale) and are naturally spaced equally in log frequency.  you can draw any reasonable shape with a sum dB of these cascaded peak/cut filters.  but nothing too wild.

you could also do this more analytically with a sum of FIRs, all time-aligned.  make a series of overlapping complementary window functions (like a Hann) but operating on log frequency instead of linear frequency.  but that frequency response can be interpolated and sampled in linear frequency and an FIR can be derived from that magnitude response with a delay of some sufficient constant value.


Visualise/compare numpy arrays from Matlab/Octave to matplotlib I'm new to python and matplotlib, and I'd like to visualise / compare 3 mfcc files stored as numpy arrays in txt format.
I have the Octave code below, and I'd like to know how it can be done using python/matplotlib.
Any help is much appreciated.
load /dir/k11.txt
load /dir/t11.txt
load /dir/a11.txt

subplot(1,2,1);imagesc(j11);axis('xy');colormap(jet);colorbar;subplot(1,2,2);imagesc(t11);axis('xy');colormap(jet);colorbar;

c=[k11(:,end),k11(:,1:end-1)];
figure(1);
Ncep=size(c,2)-1;
a=real(fft([c,zeros(size(c,1),512-Ncep*2-1),c(:,end:-1:2)]'));
imagesc(a(1:end/2,:));
axis('xy');
colormap(jet);

c=t11;
figure(2);
Ncep=size(c,2)-1;
a=real(fft([c,zeros(size(c,1),512-Ncep*2-1),c(:,end:-1:2)]'));
imagesc(a(1:end/2,:));
axis('xy');
colormap(jet);

c=a11;
figure(3);
Ncep=size(c,2)-1;
a=real(fft([c,zeros(size(c,1),512-Ncep*2-1),c(:,end:-1:2)]'));
imagesc(a(1:end/2,:));
axis('xy');
colormap(jet);
```

# in python you will need to import the modules you use
# some users with matlab/octave background prefer
# to use from numpy import *, making all the functions global
# but with the time they learn (some not) that modules are good.
# you prevent overwriting an internal function with a custom function
# for instance
import numpy as np;

# for plotting we use matplotlib
import matplotlib.pyplot as plt

# load /dir/k11.txt
# load /dir/t11.txt
# load /dir/a11.txt
## yow need to save the result in variables
k11 = np.loadtxt('/dir/k11.txt')
t11 = np.loadtxt('/dir/t11.txt')
a11 = np.loadtxt('/dir/a11.txt')

Matplotlib will provide basically the same functions
plt.subplot(1,2,1);
plt.imshow(j11); # or you can use pcolormap
plt.set_colormap('jet')
plt.colorbar()
plt.subplot(1,2,2)
plt.imshow(t11);
plt.set_colormap('jet')
plt.colorbar()

Unfortunately python does not have a way to build block matrices, if you want to do vertical concatenation you can simply build an array from a list of arrays. Also instead
of using the keyword end you simply use negative indices, indices count from 0 not from 1, and you use square brackets to differentiate indexing from function calls, and you can omit an index when using a slice ending at the end of the axis.
#c=[k11(:,end),k11(:,1:end-1)];
c = np.hstack([k11[:, -1], k11[:, :-1]])

For the above operation you could use np.roll instead.
plt.figure() # I don't use figure identifiers probably you can
Ncep=c.shape[1]-1;

# the step size in python is the third so x[::-1] is flipping x.
# by default fft operats on the last axis (axis=-1)
# if I remember correctly matlab by default
# will apply fft on the columns (axis=0)
a=np.real(np.fft.fft(np.hstack([c,np.zeros((c.shape[0],512-Ncep*2-1)),c[:, :0:-1]]), axis=0);

I think with this you already got started, in the module np.fft you will find routines like real fft (fftr, ifftr), complex fft (fft, ifft), fftshift, fftfreq. And you find documentation here


Channel Capacity - is this only for discrete/digital signals? The units of channel capacity is bits/second. Does this mean this only refers to discrete/digital signals? Is channel capacity analogous to bandwidth in case of an analog signal transmission?

The definition of channel capacity can be applied to either digital or analog cases. The meaning depends on how you calculate it. The definition is:

$$C=\max_{p_X(x)}\ I(X;Y)$$

In the digital world, $p_X(x)$ is a mass probability function and the mutual information is calculated as $H(Y)-H(Y|X)$, where $H(\cdot)$ denotes the entropy.

In the analog case, then $p_X(x)$ is a probability density function and the mutual information is calculated as $h(Y)-h(Y|X)$, where $h(\cdot)$ denotes the differential entropy.

Regarding your other question, channel capacity and bandwidth are not the same thing. Take for example a Gaussian channel, where $$C=B\log_2(1+\mathrm{SNR})$$

You can see that the capacity depends on the bandwidth, but they are clearly not the same.


How to find spectral density of a signal whose correlation depends on time? Suppose we have two uncorrelated, real white noise $\xi_1(t)$ and $\xi_2(t)$, with zero mean, and equal variance.

Now, I create the signal $\xi(t) = a\xi_1(t)\cos(\omega_0 t) + b\xi_2(t) \sin(\omega_0 t)$. ($a,b$ are real)

I want to find the spectral density of $\xi$.
The correlation $<\xi(t) \xi(t+\tau)>$ depends on both $t$ and $\tau$. Now, to calculate spectral density, should I calculate
$\int_{-\infty}^{\infty} e^{i\omega\tau} <\xi(0) \xi(\tau)> d\tau$?

Or should it be $\int_{-\infty}^{\infty} e^{i\omega\tau} (\text{average over all $t$ of} <\xi(t) \xi(t+\tau)>) d\tau$ ?

I have calculated the correlation $<\xi(t) \xi(t+\tau)>$

$=<[a\xi_1(t)\cos(\omega_0 t) + b\xi_2(t)\sin(\omega_0 t)]\times [a\xi_1(t+\tau)\cos(\omega_0 (t+\tau)) + b\xi_2(t+\tau))\sin(\omega_0 (t+\tau)]>$

$= a^2 \cos(\omega_0 t)\cos(\omega_0 (t+\tau))<\xi_1(t) \xi_1(t+\tau)>$

$+ b^2 \sin(\omega_0 t)\sin(\omega_0 (t+\tau))<\xi_2(t) \xi_2(t+\tau)>$

$+ ab \cos(\omega_0 t)\sin(\omega_0 (t+\tau))<\xi_1(t) \xi_2(t+\tau)>$

$+ ab \sin(\omega_0 t)\cos(\omega_0 (t+\tau))<\xi_2(t) \xi_1(t+\tau)>$

Since $\xi_1(t),\xi_2(t)$ are uncorrelated noise, this is equal to $= a^2 \cos(\omega_0 t)\cos(\omega_0 (t+\tau)) \cdot c\delta(\tau) + b^2 \sin(\omega_0 t)\sin(\omega_0 (t+\tau)) \cdot c\delta(\tau)$

(where $\delta(\tau)$ is the Dirac delta function, $c$ is a real number that depends on the amplitude of the noises)

Your process is not stationary. As you already correctly noted, your autocorrelation function depends on $t$ and $\tau$. Let me call it $\varphi(\tau,t)$.

There are multiple ways of dealing with such cases. One is to simply consider Fourier transforms with respect to each of the time variables, treating them independently: The transform with respect to $\tau$ gives you frequency (say, $f$), where as the transform with respect to $t$ gives you a rate of change as in how fast do your statistics change, the latter often being referred to as Doppler frequency (say $\alpha$).

Now you can define four functions:

Time-varying ACF $\varphi(\tau,t)$
Time-varying Power spectral density $\Phi(f,t)$
Delay/Doppler cross spectral density $\varphi(\tau,\alpha)$: often also called scattering function
Frequency/Doppler power spectrum $\Phi(f,\alpha)$

These are also called the second set of Bello functions, the concrete naming of each of them varies widely across sources.

Another way of attacking the problem is to go to the Wigner-Ville distribution and its variants, have a look here.


What is the phase response of SRRC/NRZ/RZ pulse shapping filters and does it matter? So I would like to know if filters usually used for pulse shaping (like SRRC) have a linear phase response, and if not, what this response looks like and does it matters in any circumstances?

FIR filters that have coefficients symmetric about their center coefficient(s) are linear phase.  Digital pulse filters commonly in use, including a Root-Raised-Cosine filter, are FIR and have this property, so they will have will have a linear phase response.

For all digital modulations, I would avoid non-linear phase response filters.


Mean of input based on acf of filtered output By looking at the ACF of the output of a filter whose impulse response I know, can I approximate the mean of the input?

For example

x = randn(1000,1);
y1=filter(ones(5,1),[1],x);
y2=filter(ones(5,1),[1],x + 0.2);
res1 = xcorr(y1,'unbiased');
res2 = xcorr(y2,'unbiased');

When I plot res1 and res2 it is clear res2 has an identical shape but shifted up. This makes me think that just by looking at the ACF I can tell which process has more 'DC' but is there a way to quantify it?

Yes. Mathematically, the limiting value of the autocorrelation function $R_X(\tau)$ of a wide-sense-stationary process is $$\lim_{\tau\to\infty}R_X(\tau) = \lim_{\tau\to-\infty}R_X(\tau) = \mu_X^2$$ where $\mu_X$ is the mean of the process. So it is not too surprising that res1 and res2 have identical shapes but res2 is shifted upwards in comparison to res1 and the above equation tells you how much the shift is.


What does PSD represent in the spatial domain? People use FFT methods to analyze topographic features in the spatial domain, such that the signal varies over (x) rather than (t). If a PSD curve characterizes amplitude versus frequency of an electrical signal, the dominant peaks of the PSD point to which frequencies contain the most power - if I understand correctly.

That being the case, what exactly does the y-axis in a PSD for a physical surface represent? Power has a very specific definition in physics. Is the geometric term Area more appropriate than Power when interpreting a spatial domain PSD curve? I'm having difficulty conceptualizing and translating the terminology and units from the temporal to the spatial domain.

DISCLAIMER: I have many questions to ask about FFT, but am new to the community and am trying to respectfully constrain the scope of my questions, so I plan to ask them separately (which I believe is the preferred format). That said, here is some background.

This is 1m-resolution multibeam sonar data, gridded to 5m, with near 100% coverage. The surface in question is coral reef seafloor, which under certain conditions forms linear features (seen below). The dimensions of these features vary with both depth and coastal aspect. My interest here is to use FFT to quantify the changes in both frequency and amplitude, so that I can assess their correlation with various environmental variables. My understanding is that if I can use FFT correctly, I can extract the dominant signal frequency/frequencies; if so, I hope to document how these dominant frequencies change from one area to the next.

The image above gives a good sense of my workflow. I start with a grid of points, which serve as neighborhood centroids for my moving-window analysis. Selecting one point (in green), I define my analysis window, calculate the average direction of features, and draw a line perpendicular to them. What you don't see is the elevation profile that's extracted, but below is the resulting PSD. In this question I'm seeking to understand the output of FFT.

I hope to execute this analysis for each grid point, and write results to file. In later questions I'll ask what metrics I can pull from FFT analysis; how exactly to tune the input data for FFT (windowing/ zero-padding/ base-2 # of samples); and whether there are best practices for determining how many spectral peaks to record.

Please see this response. Radiant Flux should take care of the physical units.

The PSD in the spatial domain would represent a "dominant" blob, whose shape is determined by the spectral content of the PSD, underlying the image. (Please note that the PSD is defined over the Fourier Transform of the autocorrelation of a signal).

EDIT:

So, your pixel size is 5mx5m and the actual length of your red line is (let's say) about 500 meters across and contains 164 samples. Let's assume a square pixel size and call it $P_s$, call the length of the line $L$ and the number of samples across the line $L_s$.

You perform an FFT on the profile line ($x_n$) and calculate the amplitude spectrum as the absolute value of the complex spectrum of $x$ (i.e numpy.abs(numpy.fft.fft(myData))). This gives you $|X_k|$.

From that, you can see one sinusoid clearly sticking out of the rest at a frequency bin of 12. What does this mean for the ridge?

The formula that links the bin to the actual spatial frequency is $f = \frac{k}{L_s} \times \frac{1}{P_s}$ where $k$ is the frequency bin.

But, what is $\frac{1}{P_s}$?

The FFT decomposes a waveform into a sum of sinusoids at different frequencies. It will take the waveform from the spatial domain (height  of ridge ($x_n$) at well defined points in space ($n$)) to the frequency domain (amplitude of sinusoids $|X_k|$ that make up the waveform over the same well defined points in space). The variable $k$ ranges from 0 (also known as DC) to $L_s$, corresponding to actual physical frequencies $f$ between DC and $F_s$.

By the way, you only show the real part of the FFT, which means that your $k \in [0 \ldots \frac{L_s}{2}]$ and your $f \in [0 \ldots \frac{F_s}{2}]$.

$F_s$ stands for Sampling Frequency and is usually expressed in Hertz (Hz). 1 Hz means 1 cycle per second. PAUSE.

Imagine a signal in the time domain, it looks like $x(t)$. For every time instance $t$ we collect a sample $x(t)$. If we sample it regularly at some sampling frequency $F_s$ then the distance (the sampling period $T_s$) between successive time instances $t,t+T_s$ is $T_s = \frac{1}{F_s}$. PLAY.

So, you know that you have sampled space once every 5m. That is, $T_s= 5m$. Then, what is $Fs$ when I know $T_s$? ... That's what $\frac{1}{P_s}$ is.

And as you can see from the above simple formulae, the longer the line, the more samples you obtain from the bottom, the more clearly you can tell spatial frequencies apart.

OK, now we know how to translate between what the FFT tells us and actual spatial frequency. For this example $f=\frac{12}{164} \times \frac{1}{5} \approx 0.014Hz$. OK...ok...What does this mean? It means that your ridge, as idealised by a single sinusoid has a period of $T = \frac{1}{0.014} \approx 71.428$ meters. That is, one period of sand raising and then dipping is about 71 meters across.

But, what about height? I read an amplitude of 42 on my single powerful harmonic, does this mean 42 meters?!?

No. Different scientific packages might "implement" slightly different formulas for the FFT. If you are within the Python ecosystem, you are likely to be using numpy. So, for this purpose, take a look at this link. You see how the forward transform is just a sum but the inverse transform has a $\frac{1}{n}$ in front of it? That's because of the way the "averaging" works. So, when you read 42, in reality it is an amplitude of $\frac{42}{L_s} \approx 0.25$ meters or about 25cm. Of course, I am assuming that your height is measured in meters here.

These are the "basics", if you like. There are a few more things that, I suppose, are going to come into play later on.

For example, is a single sinusoid good enough to represent the cross section of a ridge pattern? Personaly, I have no idea, but I am sure that there are ways for you to find out. In general though, you can order the amplitudes of the harmonics of the FFT in decreasing order and keep adding them until you have reached something like 90% of the total amplitude. This might not extend to more than 4-5 coefficients in places where the ridge is very clear. Another way to select coefficients is, of course, by estimating the goodness of fit of a set of $n$ sinusoids to explaining the total variation of the amplitude in the obtained sample, there are many different ways to carry out this step. The more coefficients you include of course, the better your representation but the more difficult the interpretation. However, you can collect these coefficients and correlate them later on with other features from your reef and see if there is something that could explain the ridge pattern better. Maybe stiff corals create strong objects which result in turbulent flow but soft seaweed or flexible flora bends and follows the current and slows it down and contributes to a more laminar flow (?).

Another thing is, how do you judge that what you "see" in front of you is good enough to be characterised as a ridge pattern. The algorithm is unsupervised so it might well go over a patch where the flow is so turbulent that it doesn't leave a well defined pattern. The mathematical pipeline will still do its work and come back with numbers but it is likely that these won't mean anything. You can easily see for example that if we were looking naively for a single powerful sinusoid from [4,4,4,4,5,4,4,4,4,4,4,], we would say 5...But that's not correct. For this reason, I would suggest that as a criterion to whether or not you go ahead with harmonic analysis of a path you first apply a Spectral Flatness metric which will tell you very quickly if the spectrum is "spikey" or close to flat. So, now, you can put a threshold on flatness and say something like "If the spectrum is too flat then don't even bother with FFT and the rest, just characterise this patch of sand as "unsuitable""...or something like that.

Hope this helps.

EDIT:

To clarify points in the ensuing discussion, I am adding here the following code with detailed comments.

import numpy
from matplotlib import pyplot as plt

#Some simple setup
signalDuration = 4; #Seconds (METERS in your case)
samplingFrequency = 100; #Samples Per Second (NUMBER OF SAMPLES ALONG YOUR SONAR PROFILE)
timeVector = numpy.linspace(0,signalDuration,samplingFrequency * signalDuration) #Discrete time (DISCRETE SPACE IN YOUR CASE)
phaseVector = 2.0 * numpy.pi * timeVector #Discrete phase. I am doing this to get rid of mul by 2*pi all the time, don't worry too much.

#Some Signal (the actual shape of it doesn't matter too much here)
Q = 7.0
A = 1.0
B = 1.0
C = 4.0

#With S here, we simulate your readings. The actual signal doesn't matter
S = Q + (A * numpy.sin(8.7*phaseVector) + B * numpy.sin(8*phaseVector)) + C * numpy.sin(2 * phaseVector)

#Let's look at our signal
plt.plot(S);plt.grid(True);plt.title("S in the time domain")

#Fourier Transform (this is the original spectrum)
FS = numpy.fft.fft(S)

#AMPLITUDE SPECTRUM. NOTICE THE SCALING FACTOR!
FS_AMP = 1/float(len(S)) * numpy.abs(FS)

#In a semilogy axis so that we can see all harmonics
plt.figure()
plt.semilogy(FS_AMP);plt.grid(True);plt.title("S in the frequency domain (AMPLITUDE SPECTRUM)")

#In a linear axis so that we can look at the amplitudes of each individual harmonic
plt.figure()
plt.plot(FS_AMP);plt.grid(True);plt.title("S in the frequency domain")

plt.show()

#The first spectral line is your C (because it is at the lowest frequency).
#The second spectral line is your B
#The third spectral line is your A
#
#Look at the spectral line that corresponds to C. We specified 4.0. Fourier Transform tells us 2.0.
#WE SPECIFIED 4.0 AND WE DO GET ABSOLUTE AMP OF 4.0 (remember the symmetry of FFT), similarly for the rest.
#Now, if you ignore Q (which is just a constant that I used to offset the measurements on purpose)
#Harmonic C is approximately 66% of the total amplitude.

Hope this helps.


Mixing two signals creates power at dc? I have a signal with components at 1200 and 1800 Hz.  When I mix this signal with a 1500 Hz signal to create components at -300 and +300 Hz (then LPF) a DC offset presents itself when the input signal is going form 1200 to 1800 Hz at a period of 1200 Hz (msk modulated signal).  This strange affect only occurs with specific time vectors.  What could cause this?  I understand how mixing two signals at different frequencies create signals at f1+f2 and f1-f2 however where is this power at DC coming from?

Below is an image to show my problem:

The first signal is the MSK signal, the second is the 1500 Hz signal I'm mixing with it, the third is the result and the fourth is the filtered version of the result to recover only the low frequency parts.

The problem I believe lies in that at certain time offsets there exists an alias of 1500 at 300 (and -300) that is destructively mixing with my msk signal and producing a dc value.  The result in the center there when you see the dc value should be a series of half sine wave pulses at 300 Hz.

Is there a way I can avoid this problem and recover my +/- 300 Hz signal from the msk signal?

If I understand correctly, you are mixing a 1500 Hz sinusoid with an MSK signal, which uses 1200 Hz to transmit 0 (or 1) and 1800 Hz to transmit 1 (or 0). The spectrum of the MSK signal, though, is quite large and includes energy at 1500 Hz. So, when mixing, you can expect to find a DC value.


STFT on time varying signal with good time and frequency resolution I am trying to determine the main frequency of a noisy signal that varies in frequency over time. Ideally I want to detect changes in the frequency as rapidly as possible - say 50Hz update rate, but I also want to get as close to the actual signal frequency as possible. The signal I am looking for varies between 100Hz and 300Hz and I am sampling the signal at 1Khz.

Currently I am using a STFT of length N applied every N/2 samples (so hop size N/2) and using a Hanning window on the input. I am then averaging the 3 overlapping frames via Welch's method to achieve a final estimate of the signal I am looking for.

I am aware that the maximum frequency resolution I can get is 1000/N so I am using Jain's method to interpolate between bins. I am also aware that the maximum time resolution I can get is 1000/N Hz. I know that these two are in tension.

Sometimes the signal is missing. In order to detect this I am calculating the standard deviation as well as the average of the three overlapping frames and using the ratio to determine the "noisiness"

My questions are:

Is this a valid procedure?
Is using Welch in this context valid given that the input signal varies?
Am I averaging the right data?
Is the time resolution really 1000/N or is the fact that I can produce output at 1000/2N Hz better?
Would more overlap help? If so how?
Is my noise calculation valid?

I have heavily relied on https://holometer.fnal.gov/GH_FFT.pdf

I think that you are splitting your signal $x[n]$ into $N$ STFT intervals which you are overlapping by 50%, then on each interval finding the Power Spectral Density (PSD) with Welch's Method.

If this is correct I think this method is valid,

but I think you might be confused by the 'overlapping' because there are two at play here. When you use Welch's method you are splitting the STFT interval into smaller intervals, from which the periodograms is calculated, and then these are averaged to give you your PSD for that STFT interval.

So there are some variables which can be changed here

The length of your STFT intervals, $N$
The overlap of your STFT intervals
The overlap of your welch's periodograms

The smaller the length of the STFT the less smooth your 'time resolution' will be because you are averaging over a smaller length of time which will allow you to pick out sudden changes in frequencies, but less data is obviously going to have some kind of negative effect (maybe lower frequency resolution, but im not sure)

I think how much you overlap the STFT is quite related to the 'time resolution' here because you can think of it like, for the first interval you've windowed you are finding the PSD over $[x[0] - x[N]]$ but then you could quite easily take the next interval to be $[x[1] - x[N+1]]$ these two clearly overlap depending on $N$, but its better to just think of you moving your STFT window along $x[n]$ and calculating the PSD at each one, so in this sense your 'time resolution' can just be the same as your sampling rate, but its smoothed by the overlapping.

The overlapping of the periodograms smooths the frequency resolution but also reduces noise

I'm not sure what Jains method is but if you wanted to look at higher frequency resolution methods I would look at instantaneous frequencies using the Hilbert-Huang Transform, see this paper for some comparisons


Difference between single channel and single microphone separation What is the difference between single channel separation and single microphone separation? Are they equal?

It really depends on what you mean by a channel. Normally a channel in a source separation scenario would the number of sensors detecting your unknown/partially known signal sources. In the audio domain sensors are referred to as microphones.

If this viewpoint is taken single channel separation and single microphone separation are equivalent. However depending on the domain the processing behind these algorithms is often different.

For the one channel (sensor) case in any signal domain, the assumptions of independent component analysis (ICA) do not hold (there are more signal sources than sensors and the number of sensors is not greater than one). Computation auditory scene analysis (CASA) can be used to address this case. Though note that the other answer is incorrect in saying that this is for only one signal source. CASA can be expanded to separate more than one source.


How many bits of precision are often used in fixed point DSP? I understand FIR & IIR digital filters are often implemented in fixed point DSP hardware. I have also seen them implemented in software. Do companies make special purpose digital filtering chips that use floating point precision (float) and (double)? When fixed point hardware is used, what are the choices for the bits of precision?

Almost all such programming is done in C, which means that 16 bits (shorts) and 32 bits (ints) make sense.  I doubt it will surprise you that most chips support those bit widths, with 16 bits being a particular favorite because it is wide enough to be useful and narrow enough to be power efficient, small, and fast.

I believe that there are a few chips that operate on 12 bits, but that is much rarer than the 16 bit chips.


What is the maximum frequency that can be captured by a given sampling rate? Let's say that we want to digitally sample an analog signal and have a sampling frequency of 8000 samples/second. Usually the maximum frequency that can be captured is the half of the sampling frequency; that is the maximum frequency that can be captured is approximately 4KHz for the example above. But what is the mathematical relationship between the two?
Thank you.

The maximum frequency is also limited by how long you sample at a given sample rate.  If you only sample for a finite length of time at 8000 sps, the bandwidth will be less than 4kHz.  A 4000 Hz signal (or any integer multiple thereof) will be ambiguous, or simply disappear.

A signal somewhat below the Nyquist frequency (depending on how long you sample at a given S/N or accuracy) can be "captured", as well as any signals between two successive integer multiples of Fs/2 given a-priori knowledge of no aliasing from any other spectra inside any other folds.


Compression Ratio in ECG compression using Wavelet Decomposition I have compressed and reconstructed an ECG signal using wavelet coefficients.
In the compression part, many of wavelet coefficients are set to zero. For example, if in original ECG $5000$ wavelet coefficients have been obtained, now $4900$ coefficients have been set to zero. I have reconstructed the ECG signal using only $100$ coefficient, but how do I measure the compression ratio?

You can get a crude wavelet "condensation" ratio of $100/5000 = 1/50$, but nothing more.

Indeed, you cannot get a compression ratio in that case, for several reasons. A stricto sensu compression ratio is given by the size of the compressed file divided by the file of the original file. But:

we do not have the size of the original file. Imagine it is a  raw raster file, whose size is: (number of samples $\times$ precision) + header size. The number of samples could be of $5000$ samples, but there is no certainty. It could be due to redundant wavelets, edge extension. And we have no idea about the original precision. Suppose $16$-bit.
we do not have a compressed file. First, the nonzero wavelet coefficients could be coded on $32$-bit floats. Then, if the data is $16$-bit, the condensation ratio is doubled. Second, you do not know the location (time-scale index) of the kept coefficients. They should be coded too. Without them, if you send your $100$ coefficients to somebody, he would not be able to reconstruct the data. Even with the location, we don't have a compressed file, which generally included quantization, entropy coding, etc.

You can have better estimates using first- or second-order entropy of the coefficients, that can give you an estimate of the theoretical limit.


What is a periodic signal in image processing? In the context of image processing (and computer vision), the concept of convolution comes up a lot. Convolution is quite related to the concept of Fourier transform and DFT. In the context of image processing, we are mainly interested in DFTs, given that signals (images) are discrete and finite. In this context, I have heard the expression "periodic signal", which reminds me of a sine or a cosine waves, which have a period of $2\pi$, that is, every $2\pi$ the output of these functions repeats. Now, in the case of images, it is not clear to me the meaning of a periodic signal. What would be a periodic image? How is this concept related to DFT and convolution?

You can think of the periodic signals like sine and cosine are basis for building other signal. See this figure for a 1-D example.

For images, sine and cosine extend to 2-D patterns as something like in here. As in 1-D signal, any images can be represented as a weighted sum of these basis images.  And roughly speaking, DFT are just weights in the weighted sum above.

Convolution connects with DFT by the nice convolution theorem. Namely, if you have two images (one of them is typically a filter, but you can consider any filter as an image patch as well) A & B, F(A * B)=F(A)F(B). That is, the Fourier transform of the convolved image of A and B is equal to the product of the Fourier transform of image A and the Fourier transform of image B.


Practical implementation of Expected Value? If i compute the average power of my input signal random variable $X(t)$ as
$$R = E[X^2(t)]$$
i.e. as the expected value of a random process, is this really just an estimate of average power? More generally, is it true that we don't have (cannot compute) the Expected Values of actual input signals, because we have only one instance of a random signal to work with? And so we are stuck with just estimates of Expected Values, instead of the real, actual Expected Values?
thanks,
js.

The answer is yes. Even if we could generate several realizations of the random process by conducting the random experiment a number of times. We could only have a limited number of realizations and for a limited period of time practically. Now to take expectation of the random process we will have to take freeze time and look at the distribution of the values across the several realizations we generated. But that would be an estimated expectation because we can only estimate the distribution, because although we have conducted the experiment several times, we cannot practically be sure to exhaust the sample space.


Looking to educate myself in dsp theory I'm  currently studying sound design, however over the past couple of years i started creating audio effects/instruments. I have also been programming in C++ for quite some years.
I noticed my shortcomings though; i can design compressors, envelopes and simple stuff like additive synthesis quite well along with basic signal knowledge (like samplerates, nyquist theorem), however i don't understand the theory behind some of the more dense mathematical subjects like filter designs, transfer functions, fourier transforms, complex numbers, time/frequency domain etc., or even the mathematical notation of such.

Without sounding too oblivious, what are the relevant fields for educating one in such topics? Also, if anyone has some book recommendations for easy introduction to required math/basic signal theory i would be most grateful :)

Regards

The following blog has download-able IPython notebooks that you can run and tinker with.

http://python-for-signal-processing.blogspot.com/

Lots of concepts covered that you can compute yourself using the scientific Python toolchain.

Good luck!


Simulate a pulse of specific power Suppose you produce a pulse x in discrete samples of duration dt each (e.g. in Matlab) according to a function.
I was given some notes in which it shows that if I want it to let it have a specific power I have to first normalize it by doing

energy = sum((x.^2).*dt) %first calculate the pulse energy
x = x./sqrt(energy) %normalize its energy to 1

then you calculate the target energy

power = 0 %target power in dBm
power = (10^(power/10))/1000 %target power in W
energy = power*len(x)*dt %target energy of the pulse

and finally you scale the pulse so that it has the target power

x = x .* sqrt(energy) %the pulse at the target power

My question is, how does that really work? Why dividing x by sqrt(E) normalize it to unit energy? And why multiplying the resulting normalized x by the new sqrt(E) brings it to the target power? Why not just divide by E for example? and then multiply by E? Does it have any relation to RMS values? And finally is it something that can be used in any arbitrary pulses, e.g. in sinusoids as well?

Well, for a signal $x[k]$, the energy is:
\begin{equation}
E=dt*\sum_{k=-\infty}^{\infty}x^2[k]
\end{equation}

You want the energy to be one by dividing $x[k]$ through a value $a$, hence:
\begin{equation}
1=dt*\sum_{k=-\infty}^{\infty}(\frac{x}{a})^2[k]
\end{equation}
which by multiplying with $a^2$ becomes
\begin{equation}
a^2=dt*\sum_{k=-\infty}^{\infty}x^2[k]
\end{equation}
Comparing this to the first equation, it can be easily seen that $a^2=E$ and since $E$ is always positive, we can conclude $a=\sqrt{E}$. And that is why you divide by the square root of $E$.

From there it should be easy to see why multiplying this with $\sqrt{E'}$ makes a signal of energy $E'$.