How to process/generate a phase shifted frequency varying sine wave? I want to phase shift an incoming sine wave with varying frequency but I am unsure how to go about doing so in practical terms.

A little more info regarding the requirements:
I have an encoder producing a sine/cosine pair with a fixed Peak to Peak but as the speed changes so obviously does the frequency. I want to effectively advance/retard the signal by altering the phase of the incoming signal.

Where I am unsure is how to go about doing this, if I use a look up table (LUT) for fast frequencies this would be ok but at slower frequencies I am going to get a very rough (digitised) signal. I was thinking I could calculate the signal 'on the fly' but again that would require having a LUT ie output = sin (input(amplitude) + offset[amplitude])

I think I maybe going about this the wrong way and there is possibly a simple solution to this?

Indeed, as Jason R sais, that's how all sine oscillators are designed for synthesizers, they must change frequency without changing phase, so you have to control them only using a counter, and you vary the counter increment speed to change the sine/cosine frequency, in that way the sine is always one count onwards from it's previous value regardles of speed changes.


Can velocity be calculated in first stage FFT? I understand that, For an FMCW radar, the velocity of a target is proportional to the phase changes at same time instants across the chirps.
I have a MxN data cube with N chirps(columns) & M samples across each chirp(rows).
So, why not first perform FFT across the chirps-to get the velocity and then perform a second FFT inside a chirp- to get the range?
In short, does it make sense to estimate velocity from first stage FFT across the chirps?

When forming a range-Doppler map using FMCW, you need to do both FFTs if you want to reliably detect a target.
Let's say we collect an MxN data matrix, where we have M fast-time samples (range dimension) and N collected chirps (Doppler dimension). Here we've simulated a target at 150 m in range moving towards the radar at 75 m/s. There is no noise and we're using I/Q mixing.
If we do only the Doppler FFT and then take the squared-magnitude (scale in dB) it looks like this:

Notice the constant vertical ridge at the 75 m/s line. Taking the FFT across every chirp's range samples, we get energy at the frequency corresponding to the target's velocity. However, we have no information on the target's range. Trying to estimate a target's velocity at this point is possible, but you don't know if this ridge is a result of one target or many.
If we do only the range FFT and then take the squared-magnitude we now get this:

A similar situation like the above, except now that we have a constant horizontal ridge at the 150 m line. Here we're showing that for each collected chirp, we've correctly determined the target's range (beat frequency) but we have no information on it's velocity.
Finally, doing both FFT's yields the target return at both the expected range and Doppler:

Now you can reliably determine velocity and be confident that it is due to a single target.


Question about the Haar Wavelet transform I'm trying to detect the presence of a sinusoid in 1/f noise conditions using the STFT and the DWT with Haar wavelets. I find this interesting phenomenon that I'm not able to explain (see plot at 1)

My sampling rate is 1200 Hz and the graphs show the one sided response (0-600 Hz). I am using 2 scales for the Wavelet Transform.

My question is, why in the FFT method can I clearly see the blips in about the correct frequency as my simulated data (180 Hz), while in the WT method, the increase in energy is spread visible across all three bands. Shouldn't I see a bright line only in the highpass segment of the second scale?

The FFT is much better at detecting a sinusoid than a DWT.  The FFT is approximating a periodic signal with a series of periodic signals.  The coefficients of the FFT will be maximum in the frequency bins (could be a single bin with sufficient resolution) where the component of the FFT series best matches the periodic signal to be detected.  The noise will distribute roughly evenly in all bins, so the signal to be detected will show up readily.

The DWT is representing the signal to be detected (in this case a periodic signal) with a set of non-periodic functions, so this requires lots of components (bands as you refer to them – but DWT really works with scale not frequency) to follow the periodic nature of the signal to be detected.  The DWT is better suited to detect / analyze non-periodic signals.


Why should the last point be excluded when performing a least-squares fit of a periodic discrete time signal? I fitted the function: f(t)=A_o+A_1 cos(wt)+B_1 sin(wt) to the following periodic discrete signal:
t=0:0.15:1.5;

y=[2.200 1.595 1.031 0.722 0.786 1.200 1.805 2.369 2.678 2.614 2.200];

Where w=2*pi / T, and T=1.5 seconds.
It happens that the fitted curve presents a phase shift with respect to the data. I consulted books on this topic and I found that I should exclude the last point of the series in order to get the right answer,...

I suspect that this matter is related with the exclusion of the last point in the calculation of the DFT, but do not find a mathematical argument that can prove this statement.
I would appreciate any mathematical explanation on this matter.

First (wrong) answer (for integrity) The $y$-value of the last point
is the same as the first one. As you apparently know the frequency,
this point comes in excess of the "fundamental period". It sounds like
this additional point comes like an implicit double-weight to the
first  point of the period.

Second take: I have tried to fit the data, with or without the last point. It seems to fit well, both cases.

%https://dsp.stackexchange.com/questions/71164/why-should-the-last-point-be-excluded-when-performing-a-least-squares-fit-of-a-d
clear;close all
%% Settings
T = 1.5; w=2*pi / T;
t = (0:0.15:1.5)';
y = [2.200 1.595 1.031 0.722 0.786 1.200 1.805 2.369 2.678 2.614 2.200]';

%% Fitting all points
ft = fittype(@(a1,a2,a3,x) a1+a2*cos(w*x)+a3*sin(w*x),'coefficients',{'a1','a2','a3'},'independent', {'x'});
f = fit( t, y, ft );
% Plot fit
subplot(2,1,1)
plot( f, t, y )
axis([t(1) t(end) 0.5 3])
grid on
title('Whole point set')
%% Fitting all points but the last
ft = fittype(@(a1,a2,a3,x) a1+a2*cos(w*x)+a3*sin(w*x),'coefficients',{'a1','a2','a3'},'independent', {'x'});
f = fit( t(1:end-1), y(1:end-1), ft );
% Plot fit
subplot(2,1,2)
plot( f,  t(1:end-1), y(1:end-1) )
axis([t(1) t(end) 0.5 3])
grid on
title('Whole point set minus 1')


Ending points of the root locus Let $$D(s) + KN(s) = 0 \tag{1}$$where $D(s)$ and $N(s)$ are polynomials of $s \in \mathbb{C}$ such that $\text{Deg}(D) = n, \ \text{Deg}(N) = m$ and $n\ge m$. The root locus method tells us how the solutions of $(1)$ changes as we change parameter $K$ from $K=0$ to $K = \infty$.
I'm trying to understand these extreme cases. Let $K \to 0$ and we have $$D(s) = 0$$ so in this case, the set of the solutions is $A = \{s \in \mathbb{R} | D(s) = 0\}$.
Now let $K \to \infty$, if we choose $s$ such that $N(s) \not = 0$ then the answer will be infinity. So we should choose $s$ such that $N(s) = 0$. In that case, if $(1)$ holds, we also have $D(s) =0$ which means $N(s)$ and $D(s)$ have the same factor but this isn't the result that should be obtained. Curiously, if we rewrite $(1)$ $$\frac{N(s)}{D(s)} = -\frac{1}{K} \tag{2}$$ and let $K\to \infty $, one possible case that $(2)$ holds is $N(s) = 0$ which gives us $m$ solutions and this doesn't require $D(s)$ have the same factor as $N(s)$! Why this happens? And why the first solution is wrong?
Example: Let $D(s) = s^2 - 4$ and $N(s) = s + 3$. So $(1)$ becomes $$s^2 - 4 + K(s+3) = 0$$If $K \to \infty $ and $s = -3$ then $9 - 4 = 0$ which is clearly wrong. On the other hand, rewriting the equation $$\frac{s+3}{s^2 - 4} = -\frac{1}{K}$$ If $K \to \infty $ and $s = -3$ then $\frac{0}{9-4} = 0$ which is true, of course.

The problem with your example is that $\infty\cdot 0$ isn't necessarily equal to zero. The only way to judge what is happening in the limit $K\to\infty$  is to divide the original equation by $K$:
$$\frac{D(s)}{K}+N(s)=0\tag{1}$$
Now it is obvious that for $K\to\infty$ the actual value of $D(s)$ is irrelevant, as long as it is finite. Consequently, the only necessary condition for $(1)$ to be true as $K$ becomes large is that $N(s)=0$.
This problem is a bit similar to the problem of determining the limit
$$\lim_{x\to\infty}\frac{x}{x+c}\tag{2}$$
The value of the limit $(2)$ is independent of the choice of the constant $c$, which becomes negligible compared to $x$. The same is true in the first equation of the question: the actual value of $D(s)$ becomes irrelevant compared to $KN(s)$.


How do I measure the amount of information within an image? I'm trying to find some information threshold that is required for performing detection of objects in images. However, I'm not sure how exactly to quantify the amount of information contained within an image.

Thoughts:

  I imagine that looking at the gradients/edges of the image would be a
  good idea, but simply counting the number of edge pixels doesn't seem
  enough. Maybe looking for patterns within the edge image?

  Perhaps I should look at the image in the frequency domain?

  I don't think I can just use Shannon's entropy because it doesn't
  really "look" at the big picture.

Any help/guidance would be appreciated.

My problem with Shannon's entropy is that it is calculated without considering
spatial structures. For example, the following two image have the same entropy score:

My problem with Shannon's entropy is that it is calculated without considering spatial structures.

That's not true!

Entropy is the expectation of Information, which is the negative logarithm of probability.

So, if you have correlated pixels (for example, in pictures of the first kind, say, a pixel in anything but the top row always has the same value as the pixel above), then that pixel's info drops to zero.

So, your statement is only true if you assume all pixels to be independent, which they aren't.


What is the minimum sampling duration to capture a given frequency? Due to the Nyquist rate, you need to sample twice as fast as the highest frequency you want to measure. Then when you do a FFT the number of samples determines the number of bins and thus the frequency resolution.

But what about the lowest frequency you want to measure?
I would think your sampling duration needs to be at least as long as one cycle of your lowest frequency.

Is this true, or can you register $20$ $Hz$ from a dozen samples at $44$  $kHz$?

If this is true, how long do you need to sample, and how do you limit the number of bins?

Using the DFT, the lowest bin is indeed corresponding to one period (i.e., the frequency which have exactly one period within the length of the DFT).

Conceptually, the DFT is correlating the input signal with different sinusoids. Now, to capture the phase information, these are "complex sinusoids", but if you knew the phase beforehand, you could just correlate with a real-valued sinusoid with the same phase to obtain the amplitude of the sinusoid.

Hence, it would conceptually be possible to register the existence of 20 Hz using a dozen sample at 44 kHz if you knew the phase of the signal (or by making multiple computations with different phase shifts). The drawback is that there are (probably/most likely) other signals which will be falsely detected as 20 Hz, because they correlate as well as the 20 Hz signal (and/or cancel the 20 Hz signal).

This is in no way a comprehensive answer, but more an insight that what the DFT and most other transforms do is to correlate the input signal to a set of basis functions, to see how they compare. With enough knowledge of your input signal you can find simpler basis functions for your particular case, but there is a trade-off between flexibility (ability to handle different cases) and simplicity (how little information/computation is needed).


Output signal length after fft I am trying to understand the functionality of fft.m in matlab .

For example: if I input a signal of length $7$, and use the FFT function, the output is of the same length.

How does that occur if the radix-$2$ takes the nearest power of $2$, which is $8$?

So how does it return an output of length $7$, assuming it operates the way I described?

There are lots of FFT algorithms, and they don't all require power of 2 sizes.  For instance, you can do any size (including prime sizes like 7) with Bluestein's (Chirp-Z) algorithm, and the algorithm is still fast ($N \log(N)$).

I've heard that MATLAB uses FFTW under the hood, and FFTW includes many (all?) of these other algorithms.  So MATLAB can do a length 7 FFT without rounding it to a power of 2.


Resampling from variable to fixed-rate I have a gyroscope which reports x,y,z values to my embedded device.
The rate that one sample appears is variable and has a jitter of around 10ms.

However I later want to have a log file with a fixed sample rate (higher than the input sample rate). I found this post What is an algorithm to re-sample from a variable rate to a fixed rate?

Unfortunately it seems way too complicated compared to what I want to do. For me it would be fine simply interpolating linearly between the samples, but I'm missing the correct algorithm to do so. Is there any common way on how to interpolate a variable sample rate to a fixed sample rate?

If you want to do it yourself, programming a linear interpolation is really simple. You interpolate between sample 1 and sample 2 by calculating the slope to generate a line between the samples (y are your measurements and x the times)

\begin{align}
m     &= \frac{y_1 - y_2}{x_2 - x_1} \\
y_1(x) &= m (x - x_1) + y_1
\end{align}

this is the linear interpolation between sample 1 and 2, it allows you to calculate the interpolated values for any point in time between sample 1 and 2. After that you just have to repeat this for sample 2 and 3 and so on. It's always the same, that means it comes down to writing one function and a loop. Just google linear interpolation and you can probably find some nice animations and youtube videos that visualize the process.

If you don't want to write the code yourself Scipy offers you multiple ways of interpolating data you can find a nice documentation here:

https://docs.scipy.org/doc/scipy/reference/tutorial/interpolate.html

the first way is probably the simplest (i just stole the examples from the page there):

from scipy.interpolate import interp1d

x = % array with the timestamps of your measurements
y = % your measurements

f1 = interp1d(x, y) % default is linear
f2 = interp1d(x, y, kind='cubic')

% xnew contains all points at which you want to sample the interpolated function
% linspace should go from your start time to your end time
xnew = np.linspace(0, 10, num=1001, endpoint=True)
ynew1 = f1(xnew) % generate the y values for all x values in xnew
ynew2 = f2(xnew)

now you can just plot(xnew, ynew) and it should look interpolated.

spline interpolation is another popular choice you can find it a bit further down with another example you could try that too if you'd like to. I personally like spline the most but none of them will magically generate data that means there really isn't a magical solution that is best, it depends on what you want to do with the result.


Filtering a signal to have a constant exponential moving average Problem

I wish to filter a discrete-time input signal, $x[n]$, so as the output signal's, $y[n]$, exponential moving average, $a[n+1] = a[n] + \alpha(y[n] - a[n])$, tracks around a constant reference $\mu_x$ such that $\mathbb{E}\{y[n]\} \approx \mu_x$, and also an approximately constant variance, $\text{var}\{y[n]\} \approx \sigma_x^2$, while maintaining as much information about the incoming signal as possible.

Naive Motivating Solution

1) Maintain an average of the input, $b[n+1] = b[n] + \alpha(x[n] - b[n])$ and an approximation of the signal variance, $c[n]$.

2) Define the output to be $y[n] = \frac{x[n] - b[n]}{\sqrt{c[n]}}$ which tracks a reference of 0 with unity variance.

Question

There are a number of problems with this approach which make it not very numerically robust, such as when low values of variance occur.

Is there a more intelligent way of doing this? Perhaps using a dynamic linear filter? Or high pass filter?

Other Points

The inputs can be considered independent, and are drawn from one of several dynamic distributions.

The application, if you are interested, is comparing the amplitudes of recent inputs with each other by assigning them a "reward" output. This is high reward for relatively higher input values and a low for lower lower values. However, I desire the mean reward to be an approximately constant output mean. Furthermore, there may be some medium-time dynamics in the inputs, such as decaying toward zero over time.

I'll draw on my (limited) experience in audio mixing and mastering. Speech recorded through a microphone requires the same statistical properties (constant moving average and constant variance), and two straightforward processing stages ensure this.

An "exponential moving average" is another word for a 1-pole low-pass filter. This moving average can be brought close to zero by passing the signal through any high-pass filter, such as a "DC killer" or "remove subsonic rumble" filter that an audio processing package such as Audacity or Audition may provide. The stronger the high-pass characteristic, the closer the moving average will stick to the center line.

Second, you want to make the variance constant. That's a job for a level compressor. Find the RMS (square root of variance) of the signal over the next few milliseconds, and adjust the gain (amplification) of the signal over that period to bring the variance up or down to the variance you want. Audio processing packages provide this as well.

You mention numerical stability. The standard DSP literature describes ways to implement a robust high-pass filter. For instance, FIR filters tend to be very stable.


SNR for the magnitude of a complex signal The SNR for a complex valued signal is presented by Dilip Sarwate in: https://dsp.stackexchange.com/a/42359
Following the same steps I would like to estimate the SNR for the magnitude of a complex signal. The received sample signal $y[n]$ is consisted of the transmitted signal $x[n]$ plus noise $w[n]$, such as,
$$y[n]=x[n]+w[n].$$
Also, it is given that, $w$ is AWGN $\sim CN(0,2\sigma^2_w)$.
When the received signal is treated as magnitude,
\begin{align}
E\{|y[n]|\}
&= E\{ |x[n] + w[n]| \} \\
&= E\{ \sqrt{ x^2_I[n] + 2x_I[n]w_I[n] + w^2_I[n] + x^2_Q[n] + 2x_Q[n]w_Q[n] + w^2_Q[n]} \}
\end{align}
And this is were I hit a road block (assuming that my math is correct). How can you extract the SNR out of this equation?
Which I don't know if it is solvable since the square rot expands all over the parameters.
Another solution I had in mind is to transfer the signals to FFT and use integrals to compare the two different magnitudes. Does that make more sense?

The magnitude $|y[n]|$ follows a Rice (aka Rician) PDF. You can find analytical expressions for the average, variance and SNR in Voltage Signal-to-Noise Ratio (SNR) Nonlinearity Resulting From Incoherent Summations and in The Squaring-Loss Paradox. These two references further investigate the SNR of the addition of independent magnitudes.


Basic Wiener application I'm trying to filter this $g(x,y)$ WGN added
image
My Wiener filter is:
$$H(u,v) = \frac{P_{f}}{P_{f}+P_{n}}$$
Where the image power spectral density $P_{x}$ is estimated as $\frac{|X(u,v|^{2}}{N^{2}}$, $X$ being the FFT and $N^{2}$ its size.
Since image and noise are uncorrelated $P_{g} = P_{f}+P_{n}$, thus
$H(u,v) = \frac{P_{g} - P_{n}}{P_{g}}$
The noise variance is known, and the noisy image PSD is estimated as shown above. The application of the filter though yields disappointing results:

Any ideas on what I may be doing wrong?

To have the desired effect of image denoising, the Wiener filter for can only be implemented with a point-spread function (blurring filter) that filters out the noise by averaging in some neighborhood of image pixels. As you do not disclose your program text, one cannot say for sure, but most probably your disappointing result comes from failure to apply the averaging in the vicinity of each pixel.
The typical derivation of the Wiener filter in Fourier domain arrives at the expression
$$
H(u,v) = {\frac {B^*(u,v)S_{ff}(u,v)} {|B(u,v)|^2 S_{ff}(u,v) + S_{nn}(u,v)}}
$$
where $H(u,v)$ is the Wiener filter in agreement with your designation, $B$ is the point-spread function (blurring filter), and I rename your $P$'s into $S$'s in order to reserve the character $P$ for a common PSD. In keeping with your use of $f$ subscripts for the original image and $g$ for this image with WGN added, $S_{ff}(u,v)$ is the (supposedly unknown) power spectrum of the original image before the noise was added.
Using the periodogram technique to estimate the power spectrum from the DFT of the observation $X(u,v)$, one can write $S^{periodogram} = |X(u,v)|^2/N^2$ and then arrive at the expression similar to your $H(u,v) = (P_g - P_n)/P_g$:
$$
H(u,v) = {1 \over B}{\frac {S^{periodogram}(u,v) - S_{nn}}{S^{periodogram}(u,v)}}
$$
with the remarkable appearance of the ${1 \over B}$ factor.
Now, there can arise problems with the inverse of the point-spread function $B(u,v)$. The solution is to use matrix pseudoinverse.
Back to your question: you failed to reconstruct a denoised image because you omitted the averaging in the pixel neighborhoods. Experimentation with MATLAB function wiener2(I,[m n]) for the pixel-wise adaptive low-pass Wiener filtering with different sizes of the pixel neighborhood would help you understand the working of the Wiener filter on noisy images. You would see how the choice of optimal m,n values depends on the local image variance.


How to read off power at different frequencies In MATLAB, consider a signal y and its Fourier transform Y:

y = A .* sin(2*pi*f_c*t + Phi);
Y = 10*log(abs(fftshift(fft(y))) / length(y));

Now, I plot the signal and its spectrum:

figure
subplot(2,1,1)
plot(t,y)
subplot(2,1,2)
plot(t,Y)

my question: what exactly did I plot in the second subplot? I know my $x$-axis gives frequency in $\textrm{Hz}$. But what does the $y$-axis give? I want to be able to read off the power of my signal at different frequencies - if possible in $\textrm{dBm}$. Since my signal-amplitude A is used to define y and not its power, I am unsure of what my $y$-axis yields.

I know my x-axis gives frequency in Hz.

That is not the case.

A DFT is a mathematical transform between discrete signals in time and frequency domain. Just as as the number of your sample doesn't give you time in seconds, the number of the DFT bin doesn't give you frequency in Hertz.

Instead, the number of the sample gives you time in "sampling intervals $T$", and thus, the bin number in the $N$-point DFT gives you frequency in "Nyquist rate fractions" $\frac {T^-1}N = \frac{f_{sample}}N$, whereby you need to pay attention – the DFT should be understood as something that shows you both negative and positive frequencies, so therange is $-\frac{f_{sample}}2$, to  $+\frac{f_{sample}}2$ (for odd and even lengths, take care of the boundaries when calculating what the frequencies of the DFT bases are). The zeroth DFT base is always the average/sum/DC of the signal.

  if possible in dBm.

Not possible, unless you know the proportionality between your digital input signal to physical power, ie. you've calibrated your measurement setup and know which power a digital $1+0j$ has. Then: Parseval's theorem is your friend.

Note that the DFT is nothing but a linear operator – a matrix that maps $\mathbb C^N\mapsto \mathbb C^N$; the question how to interpret complex frequency domain values is absolutely up to your understanding of your system. Sometimes, the values' magnitude squares can be interpreted as estimate for the Power Spectral Density at the frequency of the bin – at other times, that is not the case (especially if you don't know whether your signal is weak-sense stationary).

Sometimes, the bin of a DFT of a system is an impulse (ask a physicist about the meaning of the Fourier transform for solid state electronics), sometimes it's something like a pseudo-angle (ask an antenna designer), sometimes, it's proportional to speed of something (ask a radar engineer), sometimes, it's really just that — a complex, $N$-dimensional vector, that has no physical significance.

I'd like to quoute:

  One can Fourier transform anything—often meaningfully.
                                                          ~John Tukey


Finding $E_b/N_0$ when coding rate is $R_c$ I would like to calculate $E_b/N_0$  when the code rate is $R_c$ for BPSK signal. What I got this

$$
E_b/N_0= \frac{1}{2R_mR_c \sigma^2}
$$
and that is the case also according to this page.

$R_m=1$ for BPSK and $\sigma$ is standard deviation of the noise. Therefore I can say new $E_b/N_0$ is equal to old $E_b/N_0$ multiplied by $1/R_c$, therefore in $\log$ scale
$$
{E_bN_0}_{[\rm db]}={E_bN_0}_{[\rm db]}(\textrm{old})-10\log10(R_c)
$$
however this is different from what  MATLAB help says since according to them

CodedEbNo = UncodedEbNo + 10*log10(codeRate);

Can anyone explain where might be the problem?
And why my equation is not the same as Matlab help?

The energy per bit, $E_b$, is independent of the coding rate. Note that $E_b$ measures the energy per transmitted information bit, not per transmitted symbol.

Let's say you're willing to spend one joule per information bit, so $E_b=1$. You use uncoded BPSK, so that each transmitted symbol carries one bit of information and so it also has energy one. Let us say that $E_s=1$, where the subscript $s$ denotes "symbol". Let us further assume that you're operating at an SNR per bit of 10 dB; that is, $E_b/N_0=1/0.1=10$.

Now you decide to try a repetition code, with three repetitions per bit. With this code, a bit 1 is transmitted as three consecutive positive symbols, and a bit 0 as three negative symbols. If you're unwilling to increase the energy per bit, then each symbol must have energy equal to 1/3 of a joule: $E_s=1/3$. In this way, you still have $E_b=3E_s=1$.

Now, in the channel, the SNR of each symbol will be reduced: since each symbol has energy 1/3, then each symbol has SNR $E_s/N_0 = (1/3)/0.1 = 5.2 \text{ dB}$. However, the SNR per information bit hasn't changed; it remains  $E_b/N_0=1$.


what is the output of MFCCs? I had calculated 13 MFCCs coefficients of a speech signal and i got output in the form numerical vectors please see below figure.
My confusion whether it is in time domain or frequency domain.
suppose if i need to plot a graph of MFCCs, what should be in x axis and y axis?

TL;DR the x-axis is time and the y-axis is coefficients which means they have no units.

Let's go through the process of manufacturing the MFCC features while emphasizing measure units tracking. We start from a sampled audio signal $s(n)$ with amplitude units on the y-axis and time units on the x-axis.

First step is to frame the audio signal. An audio signal is statistically independent, meaning it is constantly changing and unmanageable. To make it approximately constant (statistically speaking) we frame it into $20_{ms} - 40_{ms}$ frames. This does not change the units so we are staying with seconds on the x-axis for each $s_i(n)$, where $i$ denotes the frame.
Next step would be to calculate the Periodogram estimate of the power spectrum for each frame $s_i(k)$. This represents the power distribution of the signal per frequency bin. the x-axis of the output is in units of frequency (Hz or $\frac{\pi \times rad}{sample}$ depending on normalization) y-axis will then respectively be $\frac{power}{frequency}$ ($\frac{db}{Hz}$ or $\frac{db}{\frac{\pi \times rad}{sample}}$). The intuition here is to imitate the human hearing system which is Identifies which frequency regions have more energy in them.
The human ear is more sensitive to frequency changes in the low-frequency region over the high one. In other words, it is easier for us to distinguish the $100_{Hz}$ frequency difference $400_{Hz}-300_{Hz}$ over $1400_{Hz}-1300_{Hz}$. To imitate this behavior, we use a filter bank to cumulate the energy over a number of bins (typically 26). We use triangular filters, to sum energies in a certain region. We use 26 filters and the higher the frequency the wider the filter. This leaves us with power units.

An example of the triangular filters:

On the next step, we approximate the results using Discrete Cosine Transform (DCT). The output of that process leaves us with 26 coefficients. We usually discard the lower and are left with 12-13. This means we have 12-13 coefficients per frame, without any units.
We produce these coefficients for each time frame and stack them together. as presented in the image below:

Note here, the x-axis represents different time frames and is measured in seconds. The y-axis is the discrete values for 12 coefficients. the level of each coefficient for each time frame, is represented by the color scheme.

An additional big advantage of this representation is also being a very compact representation of sound resulting in 12-13 numbers per time frame.


Kalman filtering of signals In the framework of the Kalman filter

How can we check errors are Gaussian and independent?
If the processes are not stationary, could that be a problem?
If I have the unobservable variables that are latent. How could I get the transition matrix or the covariance matrix for the equations?

1) How can we check errors are Gaussian and independent?

Check out this answer for a way to check for independence (whiteness).

This answer gives one way (not necessarily a good way) to check for Gaussianity: calculate the kurtosis and check that it's close to 3.

  2) If the processes are not stationary, could that be a problem?

Only if your signal model (state transition matrix, input matrix, and output matrix, process noise covariance, measurement noise covariance) do not reflect the sort of nonstationarity that your signal has.

  3) If I have the unobservable variables that are latent. How could I get the transition matrix or the covariance matrix for the equations?

How do you get the transition matrix ever? Make it up!

How do you know that there are unobservable states if you don't already have a signal model?


Up-sampling using ZOH The basic idea of up-sampling is to add $M-1$ zeros after each sample.
I need to make upsampling(factor of 3) of the function $rcosdesign(\beta, M, T)$ using the ZOH filter in Matlab.
$\beta=0.25$,
$M=6$,
$T=4$.

I didn't get a full understanding of the upsampling using this filter, so i need help on how this should be done.

Thanks in advance!

Upsampling is basically padding zeros in between signal samples and then interpolation. The interpolation itself is a low pass filter. The rcosdesign is a low pass filter so after zero padding by factor of 3, you need to convolve the resulting zero padded sequence with this filter.

Note that: it is not an ideal filter, so you would get varied attenuation for different frequencies for the discrete time frequency repsonse after convolving with this filter.

Further, remember you would be "drawing in" frequency components from other periodic repetitions of the discrete time frequency response of the original sequence immediately after zero padding, these would appear in the high frequencies of your $-\pi$ to $\pi$ spectrum hence roll off and spectrum of original sequence would decide the amount of alisaing and distortion of the original signal's frequency response.

As an example if your original sequence is bandlimited in $-\pi/4$ to $\pi/4$. Then this would shrink to $-\pi/12$ to $\pi/12$ after upsampling by 3 but you would also get frequency components to the left of $-7\pi/12$ and right of $7\pi/12$. So in this example make sure your roll of begins at $-\pi/12$ and $\pi/12$ and frequency resonse of filter decays very much to 0 before $7\pi/12$ and $-7\pi/12$

Sample MATLAB code below:

Sampling period  implied and ignore the time scale:

close all;
clear all;
clc;

t=-5:0.1:5;
T=5;
RECT=rectpuls(t,T);
FFT_size=1024;

signal_fft=fft(RECT,FFT_size);

signal_up=upsample(RECT,3);
signalup_fft=fft(signalup,FFT_size);
h=rcosdesign(0.25,64,2);
H=fft(h,FFT_size);
Product=H.*signalup_fft;
filtered_sequence=ifft(Product,FFT_size);

figure;
subplot(6,1,1);
plot(t,RECT);
ylim([-2 2]);
xlim([-5 5]);
title('original signal');
subplot(6,1,2);
plot(fftshift(abs(signal_fft)));
title('frequency response of original signal,periodicity implied');
subplot(6,1,2);
plot(fftshift(abs(signal_fft)));
title('frequency response of original signal,periodicity implied');
subplot(6,1,3);
plot(fftshift(abs(signalup_fft)));
title('frequency response of upsampled signal,periodicity implied');

subplot(6,1,4);
plot(fftshift(abs(H)));
title('frequency response of rcosdesign filter,periodicity implied');

subplot(6,1,5);
plot(fftshift(abs(Product)));
title('frequency response of filtered sequence,periodicity implied');

subplot(6,1,5);
plot(filtered sequence);
title('recovered signal, ignore time scale');


Finding length of period in time domain data I have a series of measurements of a signal source, which emits a periodic signal at an unknown interval time of p seconds. Detecting the signal is not easy so I am missing quite a few signals in the data. Data is also noisy so there are errors in the data of up to +/- 10% per signal.

Example measurements x might look like this:

72.3, 185.1, 364.2, 570.2, 679.2, 1060.7

The result I am looking for, is the period hidden in the measurments, p = ~100 in this example case and possibly the best offset c = ~70.

Typically I have 4-7 measurements in a series I would like to analyse, so it seems an analytic answer is more appropriate than a sound analysis (FFT).

Ideas I have considered:

Non-linear optimization such that we optimize for x = n * p + c. Levenberg-Marquart should do the trick.
Building a histogram of difference between consecutive measurements and picking the average of the fullest bin.
Floating point GCD https://stackoverflow.com/questions/445113/approximate-greatest-common-divisor

But I am very unsure whether I might have missed an obvious solution from the sound analysis camp (autocorrelation, Harmonic Spectrum Product, DFT, e.g. https://stackoverflow.com/questions/4716620/algorithm-to-determine-fundamental-frequency-from-potential-harmonics).

So I am turning to the wisdom of SO. How would you go about solving this problem most elegantly? Library suggestions are fine (C++).

This really reminds me of the problem of beat induction in music signals, in which you have to infer a latent periodic grid from note onset measurements. The grid is "latent" - it's not because a track is at 120 BPM that it strictly has a quarter note every 500ms - the composer can choose to phrase the melody with half notes, then eighth notes... And then there's accentuation (which would change the intensity of pulses) and swing (which would shuffle them a bit temporally)...

If you really want to work from what you describe in your question - the list of timestamps - you need to compute a histogram of time deltas, convolve it by a gaussian window to smooth it (to consolidate together deltas which are close to each other), and then evaluate as possible period candidates integer divisors of the histogram peak. This is described in section 4.3 of this paper.

While it might work well in your case, I strongly encourage you not to detect events, but rather try to estimate the period from your underlying raw signals. The reason is that the detection process is discarding some valuable information which might help the period estimation process. For example, you record "spikes" at 72.3, 185.1, 364.2, 570.2, 679.2, 1060.7... But what is there at 270? It will be useful for the periodicity detection system to be informed that there might be a little spike at 270 just below your event detection threshold.

My own experience with beat tracking / induction is that it works much better "pre-detection" (with a continuous onset detection function keeping "weak" source of informations such as ghost notes) than "post-detection" (where all "weak" events have been thresholded).


Help me understand the stages involved in filtering a signal using Discrete Fourier Transform I have a series of discrete values measured from a sensor. I want to filter the frequencies coming from this sequence of values. Then, if I understood the process correctly this is what I do:

I create a discrete fourier transform of the values
I identify the bins that correspond to the frequencies I want to remove from the original signal using the formula freq = (k * FPS)/N, where k is the bin number (starting at zero), FPS is the frames per second the signal is being captured and N is the number of samples.
supposing I want to remove from the signal every frequency below 10 Hz and the 9th bin is equal to 10 Hz, then I zero, all the real and imaginary bins from 0 to 9 of the DFT result.
then I reconstruct the signal using this inverse DFT results with certain bins zeroed (filtered).

If this process is correct, I do not understand one thing:

In my original signal I get only real values. I input these real values into the DFT algorithm using zeros for all imaginary parts. I get real and imaginary from the DFT. I filter the whole thing and do an inverse DFT. The final results is real and imaginary.

How do I get a real only signal after the inverse DFT? How do I get rid of the imaginary part and get a real only result?

A purely real time domain signal has a frequency transform such that positive and negative frequencies are complex conjugated of one another.  If you are scaling the positive bin, but not the negative, you are implementing a complex filter kernel, resulting in a complex output.

Side note, you can’t just change the magnitude of a frequency bin in a frame and do an inverse transform in real time processing.  The phases of adjacent bins won’t line up.  This is a more advanced topic beyond the scope of your question.


How to find peaks and periods? I am looking at data like this:

Source data: http://pastebin.com/raw.php?i=L6cd8d5K
The data is quasi-periodic, there is a distinct cycle but it is not an exact repeat, and the period changes considerably.  I would like to find accurate peak positions, and period durations, for each cycle.  Ideally the peak positions would match the intuitive notion of where the exact peak is.
What makes this hard:

noise: although it may not be entirely noise!  the variation is much greater when the values are high, and so the highest spikes tend to occur near where the "real" peak is.

the waveform is not symmetrical: neither left-to-right nor up-down.  doing something like a highpass filter will move peaks up and down by an amount that depends on their duration and the shape in between.

What is unlikely to work well:

bandpass filter combined with conventional peak picking: the position of the peaks will shift since values from the (asymmetrical) up and down slopes is added to where the peak is.  The really high values will not have a high weight in the result since the lowpass part of the bandpass filter will spread their energy to adjacent values.

autocorrelation: the shape of successive cycles is too different.

period durations measured from the leading (or trailing) edge of one cycle to the next: again shapes are too different.

Personally, I would just use a lowpass filter  try and keep it as small as possible to reduce the effect of asymmetry.

The fact that there is significant noise already implies some uncertainty in the peak position, I expect this will be of similar or greater magnitude to any bias the filter introduces.

If you are really concerned about the bias you could look at using morphological opening and closing filters. This is basically running a ball over/under the surface so either one will over/underestimate the peak height. You can then get 2 lines to give an envelope of where the peak may lie.  The average of the two curves is probably a decent estimate of the "true" curve prodived the noise frequency is suitably high.


Fourier transform of unit step I was reading pdf by caltech and in one of its section, Fourier transform of Unit step signal is calculated but I am confused, how this can be possible if region of convergence for Laplace transform ($1/s$) of unit step signal does not contain imaginary axis?

And if above case is possible then if it given that impulse response of a system is Unit step then frequency response of it should also exist and equal to          $H(ω)= πδ(ω) + 1/jω$ and then can we
calculate Fourier transform of output by computing $H(ω)X(ω)$ where $X(ω)$ is Fourier transform of input?

The Fourier transform can be generalized for functions that are not absolutely integrable. We can define a Fourier transform for functions with a constant envelope (e.g., sine, cosine, complex exponential), and even for functions with polynomial growth (but not with exponential growth). In these cases we must be prepared to deal with generalized functions in the expression of the Fourier transform, such as the Dirac delta impulse or its derivatives. This is also true for the Fourier transform of the step function.

The multiplication property of the Fourier transform remains true, albeit with certain restrictions, so we can generally compute the Fourier transform of the convolution of two functions by multiplying their Fourier transforms, provided that the convolution exists.


How was this result on discrete Fourier series achieved? I was trying to do the question 10, part b of the following document (https://ocw.mit.edu/resources/res-6-007-signals-and-systems-spring-2011/assignments/MITRES_6_007S11_hw10.pdf)

I was going through the solution(https://ocw.mit.edu/resources/res-6-007-signals-and-systems-spring-2011/assignments/MITRES_6_007S11_hw10_sol.pdf). Can someone explain to me how they reached from step 3 to answer to part b) of the 10th question. I understood till step-3 of the solution.

Note that $x[n]$ has period $N$ and the sum you refer to has $NM$ terms, so the sum has the form

$$\sum_{n=0}^{NM-1}x[n]f[n]=x[0]f[0]+x[1]f[1]+\ldots+x[N-1]f[N-1]+x[0]f[N]+\ldots +x[N-1]f[2N-1]+x[0]f[2N]+\ldots=\\=x[0]\big(f[0]+f[N]+f[2N]+\ldots+f[MN-N]\big)+x[1]\big(f[1]+f[N+1]+\ldots+f[MN-N+1]\big)+\ldots+x[N-1]\big(f[N-1]+f[2N-1]+\ldots+f[MN-1]\big)=\\=\sum_{n=0}^{N-1}x[n]\sum_{l=0}^{M-1}f[n+lN]$$

With $f[n]=e^{-jk2\pi n/NM}$ we obtain

$$\frac{1}{NM}\sum_{n=0}^{NM-1}x[n]e^{-jk2\pi n/NM}=\frac{1}{NM}\sum_{n=0}^{N-1}x[n]\sum_{l=0}^{M-1}e^{-jk2\pi (n+lN)/NM}\tag{1}$$

For the other sum over $y[n]$ you can do exactly the same. The final result is obtained as follows. The right-hand side of Eq. $(1)$ can be written as

$$\frac{1}{NM}\sum_{n=0}^{N-1}x[n]\sum_{l=0}^{M-1}e^{-jk2\pi (n+lN)/NM}=\frac{1}{NM}\sum_{n=0}^{N-1}x[n]e^{-j(k/M)2\pi n/N}\sum_{l=0}^{M-1}e^{-jk2\pi l/M}\tag{2}$$

The last sum over $l$ in Eq. $(2)$ is given by

$$\sum_{l=0}^{M-1}e^{-jk2\pi l/M}=\begin{cases}M,&k=mM,\;m\in\mathbb{Z}\\0,&\text{otherwise}\end{cases}\tag{3}$$

Consequently, for $k=mM$, Eq. $(2)$ evaluates to

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

where $a_k$ are the discrete Fourier series (DFS) coefficients of $x[n]$.
For other values of $k$ the expression in Eq. $(2)$ equals zero. Exactly the same manipulation can be done with the sum over $y[n]$.


What are the types of noise (FM,AM, MSK modulated) that I can use for jamming GSM signals? I need to jam GSM signals in-order to prevent mobiles in near-by vicinity from connecting to the BTS. I am employing a noise based jamming strategy. What type of noise will be most appropriate for this task?

Almost any kind on in-band interference (like Gaussian noise) will mess up the GSM signal and prevent its correct reception. Note, however, that:

Transmitting on the GSM band without a license is illegal in many locations.
Tampering with cell phone operation in general is also illegal.

If you insist on doing this, one approach would be to set up a GSM base station using OpenBTS in an environment where your signals won't interfere with other base stations and won't be otherwise detected (for example in a laboratory with RF shielding). You can use unlocked phones to connect to your OpenBTS base station, and then measure the effect of your interefernce on their performance.


Quantization error in sampling Consider the analog signal,
$$x(t) = 2\cos(3000\pi t)+3\sin(4000\pi t)+7\cos(6000\pi t)$$
If the sampling rate is 8000 samples per second and quantized at 8 bits, find:

Discrete values at any two points
Quantization error at those points

I have been trying to solve this with no success so far. How can I solve this?
So far what I have done is convert $x(t)$ which is continuous time signal to $x[n]$ which is discrete time signal.

$$x[n]= 2\cos\left(\frac{3\pi}{8}n\right)+3\sin\left(\frac{3\pi} {2}n\right)+7\cos\left(\frac{3\pi}{4} n\right)$$ where $n$ is integer.
By any two points what I mean is for any two values of $n$. For example, suppose I calculate $x[n]$ for $n=1$ and $n=2$. How do I find quantization error at those points?

From your question it is unclear what is the dynamic range exactly, or is it needed to normalize the signal or not? Anyway, I assumed it to be normalized (since somehow this happens in ADCs). To quantize the normalized signal, you can use the following equation :
$$x_q[n]=\frac{\bigg\lfloor\left(\frac{x[n]}{\textrm{Dynamic Range}}\right)\cdot2^\textrm{Number of Bits}\bigg\rfloor}{2^{\textrm{Number of Bits}-1}}$$
where, $\textrm{Dynamic Range}$ here is $2$, $\textrm{Number of Bits}$ is $8$.

Quantization error in each sample would simply be the difference of the original signal and the quantized version of it.

$$Q_e=x[n]-x_q[n]$$

I can't provide the complete solution, however the following code gives you the idea of application of above discussions.

fs=8000;
Ts=1/fs;
f1=1500;
f2=2000;
f3=3000;
t=(0:512)*Ts;
x=2*cos(2*pi*f1*t+13124)+3*sin(2*pi*f2*t+pi/8)+7*cos(2*pi*f3*t);
x=x./max(abs(x)); %Normalization

Range=2; % Dynamic Range ; from -1 to 1 I assume
Q=8; %Number of Bits

Quantized_Signal=floor(x./Range*2^Q)./2^(Q-1);


Matched Filter of Gaussian Signal How can we evaluate the matched filter's impulse response of a Gaussian function $x(t) = \exp(-\frac{t^2}{2})$. As far as we know, for a signal of finite duration $T$, the impulse response of it's matched filter is $h(t) = x(T-t)$. But as Gaussian function is of infinite duration, how can we find the matched filter ?

If you look at:

https://en.wikipedia.org/wiki/Matched_filter

the derivation does not depend on the filter being either causal or have a finite duration. The derivation is fully satisfactory for convolution for time index $k$ over $-\infty,\infty$.  It's all about what filter will maximize the Schwarz inequality.

In this case
$$
x(t)=\exp(-\frac{(t-0)^2}{2}) = x(0-t)
$$

Functions that are symmetric are their own matched filter.

The $T$ parameter enters the discussion in terms of when does one sample the filter for the purpose of detection.  If $x(t)$ was delayed by some unknown time(nonsynchronous detection),  I would use the matched filter to estimate it's delay.


Inverting a sampled system I'm doing some self-study for an upcoming exam and came across the following question:
My first idea to solve it was using the bilinear transform to get some approximation of $H(Z)$ (or just using the plain sampled version of the $H(s)$) then $V(Z)$ would be $H^{-1}(Z)$ to invert the effect of $H$.
Is this the correct approach? If yes, is there any other approach to this type of problem, like not using inverses?
Thanks for the help!

Can't be done.  Next question?
First, the author does not define $x[n]$.  Let's assume, then, that $x[n] = x(nT_s)$, where $T_s$ is the sampling interval; i.e. $f_s = 1/T_s$.
So assume that $x_1(t) = x_1[n]\ \delta \left(t - \left(n + \epsilon\right) T_s \right)$, where $\epsilon$ is any number where $0 < |\epsilon| < 1$, and $x_2(t) = x_2[n] \left( u \left(t - T_s \left(n - \frac{1}{2}\right) \right) - u \left(t - T_s \left(n + \frac{1}{2}\right) \right)\right)$.
Basically, $x_1(t)$ is a string of Dirac delta functionals that do not land on $t = n T_s$, where $x_2(t)$ is a string of rectangular pulses (or a "stair-step wave") where the changes don't land on $t = n T_s$.  I contend without proof* that for any given signal $b[n]$, both an $x_1[n]$ and an $x_2[n]$ could be found that would cause the system to generate $b[n]$.  These $x_1[n]$ and $x_2[n]$ would have the following properties:

In general, $x_1[n] \ne x_2[n]$.
In the case where $x(t) = x_1(t)$, $x[n] = x(nT_s)$ would be zero for all $n$
In the case where $x(t) = x_2(t)$ (if I've dotted my 'i's and crossed my 't's right), $x[n] = x(nT_s) = x_2[n]$ would be true for all $n$.

Since you can have more than one version of $x(t)$ that leads to the same $b[n]$, there is no inverse.  You can define some constraints on the form of $x(t)$ that depend on some $x[n]$, as I did above, and if you do it properly you can then find an inverse system from $b[n]$ to $x[n]$ -- but because the sampling process loses information, you can never reconstruct some generic $x(t)$.
It could be that the author was just implicitly thinking (or stated explicitly, and you missed it), that $x(t)$ was constrained to be bandlimited to $f_s/2$ (i.e., the Nyquist criteria).  If that's the case, and if the relationship between $x(t)$ and $x[n]$ were defined then you could invert the system.
* Proof is left an an exercise for the reader, heh heh heh.


Audio descriptor of data Is there a way to generate a waveform from data (any kind of data: integer, hexadecimal, string, char, binary), in a common audio player, without converting the data to an audio format?

I'm looking for a way to "code" an arbitrary file in to audio (obtain an audio-descriptor of the data).

If you want a simple way to transform a numeric array in sound, I suggest you to use Python, with the wave module and Numpy numeric package.

Now if you want to analyze the spectrogram of the numbers without listening to anything, then Python's Matplotlib does what you want.

Below is a sample working code that generates a simple tone, saves it to a .wav file, AND displays a spectrogram:

import numpy
import wave
from matplotlib import pyplot as plt

samplerate = 44100

# generate simple sound
sound_data = numpy.sin(numpy.linspace(0,20000,100000))*1000

# converts to a string representation (I suspect there might be a more natural way to do this)
raw = "".join((wave.struct.pack('h', item) for item in sound_data))

# saving to .wav file
filename = 'aaa.wav'
sound = wave.open(filename, 'wb')
sound.setparams((1, 2, samplerate, 1, 'NONE', 'noncompressed'))
sound.writeframes(raw)
sound.close()

# plotting spectrogram
plt.specgram(sound_data)
plt.show()


Zero padding effect on a FFT of gaussian noise I have a gaussian noise $\nu(t)$ with variance $\sigma^2$. After a FFT I get $X(\omega)$. If now I do the IFFT on the $X(\omega)$ can I say that the result is still a gaussian noise of variance $\sigma^2$? What is the effect of zero padding on FFT? How the zero padding affects the statistics of the noise after FFT and IFFT?
Thanks.

Think about both questions separately. First of all, the (I)FFT is just an implementation of the (I)DFT, so I'm going to generalize all this to the DFT.

Does the zero-padded IDFT retain variance?

Parseval's theorem says power out = constant factor · power in, and the power of the zero-padded sequence is the energy of that sequence divided by it's length – and that length is larger than the original length, whereas the energy stayed the same.

Does the zero-padded IDFT retain gaussianness?

Long story short: yes. This is a result from the fact that the DFT is effectively a large sum of sufficiently identically distributed random variables. What you need is independence (which would imply whiteness), but as mentioned below, most people mean "WGN" when they say "GN".

Other effects?

When people say "Gaussian noise" they often mean "white Gaussian noise", but since that would have a constant PSD, and you explicitly made it so that your noise realization's Fourier transform is anything but constant, but comes in a boxcar shape, you obviously lose the whiteness.


Solving Scale and Permutation indeterminacies in Blind Source Separation Many papers using blind source separation algorithms solve the inherent scale and permutation indeterminacies of the estimated mixing matrix using an "optimal method", however this method is never explained. How does this work?

Thanks in advance

What i know in implementations of BSS (e.g jade) is that they first normalise the input data (to variance = 1) and then order the data according to scale. So this is unique and will produce consistent output ordering for the same input data (will add references shortly).

Have in mind however that this is just a consistent heuristic (any consistent heuristic will do) as BSS methods by definition are clue-less about scaling and permutation of the inputs.

References of permutation problem for jade:

The following heuristic is described in the comments section of the matlab version of jade:

% Usage:
%   * If X is an nxT data matrix (n sensors, T samples) then
%     B=jadeR(X) is a nxn separating matrix such that S=B*X is an nxT
%     matrix of estimated source signals.
%   * If B=jadeR(X,m), then B has size mxn so that only m sources are
%     extracted.  This is done by restricting the operation of jadeR
%     to the m first principal components.
%   * Also, the rows of B are ordered such that the columns of pinv(B)
%     are in order of decreasing norm; this has the effect that the
%     `most energetically significant' components appear first in the
%     rows of S=B*X.


Filtering out a specific frequency in an analog signal I'm wondering how I would go about filtering out a regular sine-wave signal with a constant frequency and voltage.

I'm really new to the whole DSP-area but imagine I would have a regular sinewave at a specific frequency (let's say 30Hz), and on top of that signal I have a bunch of other small peaks (which are the ones I would want to sample and read). What method would I choose?

Here is a visualization of what I meant:

Excuse my terrible Paint-skills but how would I filter out that sinewave, and just keep the small spikes in voltage? What might be worth sharing is that the time between the spikes does have an impact on my measurement of the spikes, or to rephrase: The amount of spikes within x-units of time does have an impact.

Any guidance to the right direction would really be appreciated.

Thank you.

A notch filter is a good approach in general for removing an unwanted frequency. But if you are certain of the exact frequency, phase and amplitude of the unwanted sine wave, then you may also consider adding the noisy signal to a sine wave of the same frequency and amplitude but 180 degrees out of phase with the unwanted sine wave.

To illustrate the notch filter example, the plot below shows an attempt to generate a sine wave with some odd impulses (black), and then filter it with a notch filter knowing the exact frequency (red).

R code for implementing it below.

# 24782
# install.packages('signal')

T <- 1024
t <- 0:(T-1)

fs <- 1000
omega <- 2*pi*30/1000
phi <- 0.9279835

x <- sin(omega*t + phi)

changed <- sample(1:T, 10)

y <- x
y[changed] <- y[changed] + 3

alpha <- 0.9
num <- c(1, -2*cos(omega), 1)
den <- c(1, -2*alpha*cos(omega), alpha*alpha)

yf <- signal::filter(num, den, y)

plot(t,y, type="l", lwd=10)
lines(t,yf, col="red", lwd=4)

legend(850, 4, c("Original", "Filtered"), lwd=c(2.5,2.5),col=c("black","red"))


LoRa sensing using USRP GNURADIO I made a simple RF sensing using USRP N210. I am trying to detect LoRa transmission.
A LoRa Node is operating at 915 MHz with Spreading Factor 8 and bandwidth 500 kHz.
I have used a USRP UHD Source with sampling rate 500 kHz, I plotted a Waterfall diagram to show the real time visualization.
Only when I tune the USRP center frequency to 930 MHz I see some LoRa spreading spectrum. However it does not match the bandwidth at which it is operating.
The variation in center frequency maybe due to clock mismatch, however I am not sure about the bandwidth issue.
Any reason why the result is showing like this?

500 kHz (watch the capitalization: Khz is wrong!) is a very low sample rate for the N210. In fact, the N210 has to decimate by a factor of 200 from its ADC's sampling rate, and that might introduce unexpected roll-off and aliasing, because that over-stresses the CIC decimator (if this doesn't make too much sense to you: don't worry, it's an implementation detail. Just: rather go for 1 or 2 MS/s sampling rate and do the rest of the decimation in software).
Also make sure you're actually seeing the signal you think you're seeing: your USRP AND your LoRa transmitter might both have frequency offsets. But: 15 MHz difference at 915 MHz? That would be an error of 1/61 ; the USRP has a frequency accuracy of ca 10 ppm (that's 1/100000), your LoRa transmitter maybe of 50 ppm. No way you're observing a 915 MHz signal at 930 MHz. Either this is something else, or you're misinterpreting something and your transmitter isn't trying to transmit at 915 MHz!
It's hard to tell from your slightly unusual choice of color scheme whether your spectral broadening happens from an overdrive of the receive amplifier. Reduce gain – if it gets better from that, that's happening.
Also make sure you're not actually seeing what you want to see – after all, LoRa is chirp-spread spectrum, and you're sampling these chirps at more or less irregular intervals. How does your spectrum change when you change the FFT size?


Important differences between BCJR and SOVA algorithms? BCJR algorithm needs noise variance for Gamma (branch metrices) calculation. So I am wondering does SOVA (soft output Viterbi algorithm) needs noise variance?

For an AWGN channel, the exact computation of the log-likelihood ratio requires knowledge of the SNR. However, in Hagenauer's original paper on the SOVA there is an approximation which allows the implementation of the algorithm without estimation of the SNR.


LTI system output given input and frequency response The question I'm trying to understand is as follows: A linear time-invariant continuous-time system has the frequency response function $$H(\omega)=\frac{1}{j\omega+1} $$

Compute the output response $y(t)$ for $-\infty<t<\infty$ when the input $$x(t)=\cos(t),\quad -\infty<t<\infty$$

I know the answer, but I don't understand how to get there.

I know $H(\omega)=\frac{Y(\omega)}{X(\omega)}$, and from that $y(t)=h(t)*x(t)$. Taking the inverse Fourier transform of $H(\omega)$ I get $e^{-t}$. If I try to compute the convolution I just get a divergent integral. Is this wrong? The answer is supposedly $$y(t)=|H(1)|\cos(t+\angle H(1))$$ but this seems to come from nowhere. Can someone please explain?

Apply the eigenfunction property of the LTI system after decomposing the sinusoidal input by Euler identity.

The eigenfunction property of the LTI system states that

$$ x(t) = e^{j \omega_0 t} \implies y(t) = H(\omega_0) e^{j \omega_0 t} $$

where $H(w)$ is the frequency response of the LTI system.

Euler identity states that :

$$ x(t) = \cos(t) = 0.5 \{ e^{jt} + e^{-jt} \} $$

Then due to system being LTI, you can define the output as:

$$ y(t) = 0.5 \{ H(1) e^{jt} + H(-1) e^{-jt} \} .$$

Since the system impulse response was a real function, its frequency response $H(w)$ will be conjugate symmetric; i.e., $H(w) = H^*(-w)$.

Let's call $H(1) = A e^{j\phi} $ where $A = |H(1)|$ is the magnitude and $\phi = \angle{H(1)}$ is the phase of the frequency response at the frequency $\omega = 1$.

Hence the output will be :
$$ y(t) = 0.5 \{ H(1) e^{jt} + H^*(1) e^{-jt} \} $$

consider the term $z = H(1)e^{jt}$  as a complex variable, then the second term will be $z^*$, hence the output is:

$$y(t) = 0.5 \{ z + z^* \} = \mathcal{Re} \{ z \} $$

$$y(t) = \mathcal{Re} \{ A e^{j\phi} e^{jt} \} = A \cos(t + \phi) = |H(1)| \cos(t + \angle{H(1)})$$


Why Wavelet developed when we already had Short-time Fourier transform If we already had Short-time Fourier transform for better analysis of a signal than Discrete Fourier Transform, then what was the need that leads to development of Wavelet Transform ?

The short-time Fourier transform doesn’t offer better analysis of data than the discrete Fourier transform, it offers a different kind of analysis.  The DFT offers an exact decomposition of data to a frequency representation. The STFT offers an approximate decomposition to a time/frequency representation.  Which is better depends on what you are after.  The same holds true of the Wavelet transform.  Wavelet transforms can be thought of as decomposition to a time/frequency representation, but wavelet transforms generalize the concept of decomposition.  Different wavelet functions have been created so you can choose a decomposition that suits your needs.


Implementing overlap add method for convolution with a slice of the filter In overlap add method for convolution a filter $x$ of length $M$ will get convolved with a signal $y$ of length $N$ where $M << N$, i.e., $z = x * y $.

Here the signal $y$ is sliced into subparts each of length $M$.

Can I slice $y$ into subparts of size $P < M$ and then do the convolution by overlap add method. What should be the algorithm for this.

The method I have tried for overlap add method is given here(in MATLAB).
In the above code, in the for loop, I wish to have filterLength < 16.
How should I modify this code.

What you are describe is typically called "Partitioned Convolution".  The benefit is reduced latency at some cost to efficiency (unless the larger size FFT is considerably less efficient e.g. due to cache size).  This reduction in efficiency may be outweighed by the increased ability to use multiple cores for the parallel tasks.


Direct sequence spread spectrum modulation When the receiver gets the signal $m_t$ sent by the transmitter, the signal will have an additive interference $i_t$. Thereby signal received will be : $$r_t = m_t + i_t$$

I know that to get back the data from the signal, we have to multiply by the pseudo-noise sequence (PN).

That will take out the pseudo-noise from the signal $m_t$ allowing us to get the data.

What do we do with the noise $i_t$?

It has been multiplied by the PN but I supposed we cant take out the data from the $r_t$ until we eliminate the noise. I don't really see what kind of signal do we get multiplying the additive interference by the PN. I've read that we use an low pass filter but i dont really undestand why.

Multiplying $r_t$ by the pseudo random spreading code does different things to $m_t$ and $i_t$: Because $m_t$ is correlated to the spreading code, it is despread while keeping its power, so the power spectral density will be higher in a small bandwidth (due to remaining modulation). On the other hand, $i_t$ is not correlated and is translated to wideband noise again.

This gives you the opportunity to filter with a bandpass and just use the band, where the signal to noise ratio is above unity.

The noise is not eliminated, but its overall power is reduced, because your bandpass rejects some of the noise (that is still wideband) and passes all of the signal (that is now narrowband).


STFT overlap add method - What am i doing wrong? i am trying to implement a short time fourier transform of an audio signal to do some filtering (subtract another stft signal using spectral subtraction). I am using a hamming window function with 50% overlap. Additional the size of my FFT is twice the window size, zeros at the end of the signal in time domain.

mic is the input and output signal for in-place filtering

here is my code for so far:

    stftWindowSize = 256;
    stftWindowSizeOver2 = stftWindowSize/2;
    stftN = 512;
    stftLog2n = log2(stftN);
    stftNOver2 = stftN/2;
    stftSetup = vDSP_create_fftsetup (stftLog2n, FFT_RADIX2);

//length is the complete signal to be processed
for(i = 0;i<(length/stftWindowSizeOver2)-1;i++){

    //apply hamming window and copy to complex arrays
    vDSP_vmul(&mic[i*stftWindowSizeOver2],1,stftHamming,1,stftTmpStorage,1,stftWindowSize);
    vDSP_ctoz((DSPComplex *)stftTmpStorage, 2, &stftMic, 1, stftWindowSize);

    //fft
    vDSP_fft_zrip(stftSetup, &stftMic, 1, stftLog2n, FFT_FORWARD);

    //here shold the filtering happen right?

    //inverse fft
    vDSP_fft_zrip(stftSetup, &stftMic, 1, stftLog2n, FFT_INVERSE);

    //scale back
    float scale = 2*stftN;
    vDSP_vsdiv(stftMic.realp, 1, &scale, stftMic.realp, 1, stftNOver2);
    vDSP_vsdiv(stftMic.imagp, 1, &scale, stftMic.imagp, 1, stftNOver2);

    //mic data back to real, clear original signal and add overlapping data back in
    vDSP_ztoc(&stftMic, 1, (COMPLEX *)stftTmpStorage, 2, stftWindowSizeOver2);

    //set old signal to zero in non overlapping area and add
    if(i==0)
        vDSP_vclr(mic, 1, stftWindowSize);
    else
        vDSP_vclr(&mic[(i+1)*stftWindowSizeOver2], 1, stftWindowSizeOver2);

    vDSP_vadd(stftTmpStorage,1,&mic[i*stftWindowSizeOver2],1,&mic[i*stftWindowSizeOver2],1,stftWindowSize);
}

basically what im trying to do here is:

get chunk of float data of length: stftWindowSize (only half the data is "new" due to 50% overlap with previous chunk) from input/output signal mic
apply hamming window
compute FFT of the chunk + zeros at the end of signal with overall length: stftN
compute inverse FFT
scale back as a result of the FFT
clear non overlapping area in innput/output signal mic.
add my processed iFFT data back into mic with a 50% overlap with previous chunk data.

My Problem:
as i'm not doing anything with my FFT-data the sound should be the same as before, but unfortunately it isn't. It sounds somewhat modified and especially low frequencies look significantly changed when i look at the spectrum of the signal before and after my STFT.

What am i doing wrong? have i misunderstood the zero padding? is this method suited for filtering (especially spectral subtraction)?
i hope you can help me!

I think, the problem lies here:

  scale back as a result of the FFT
  clear non overlapping area in innput/output signal mic.

Why are you doing this? You are using a Hamming window with 50% overlap, summed up, this yields a constant one. If you delete the "non overlapping" part, you effectively change the window you are using to a weird "half hamming" one, this is bound to yield incorrect results. Skip these two steps, and you should be fine.


How to apply Butterworth High Pass Filter I have a problem with applying Butterworth High Pass Filter to my data. I would like to print filter for Bx and By matrix. As you can see I have both positive and negative values,

  how to apply math.fabs() to Bx and By to get only positive
  values?

For my high pass filter I have those requirements:

Fc = 2 Hz

I would like to cut off values below 100 pT.

A part of my current code is:

plt.ylabel('Pico Tesle [pT]')
plt.xlabel('Time [ms]')
plt.grid()
plt.plot(time[51:-14], Bx[51:-14], label='Canal 1', color='r', linewidth=0.1, linestyle="-")
plt.plot(time2[1:-14], By[1:-14], label='Canal 3', color='b', linewidth=0.1, linestyle="-")
plt.legend(loc='upper right', frameon=False, )

And a chart:

UPDATE:
I have used this function to generate absolute values for Bx and By matrix.

plt.subplot(413)
np.absolute(fft1)
plt.plot(time[51:-14], np.absolute(fft1), color='r', linewidth=0.1, linestyle='-')
plt.grid()

plt.subplot(414)
np.absolute(fft2)
plt.plot(time2[1:-14], np.absolute(fft2), color='b', linewidth=0.1, linestyle='-')
plt.grid()

What I received thanks to that is that all my measurements (here are atmospheric discharges) are only in positive values (pT - picoTesla). First plot shows Canal 1, second plot shows Canal 3, and the third plot shows both Canals (Channels) combined.

No, I need (I guess) use High Pass Filter to cut off all measurements below 100 pT. Any ideas?

If you just want to do thresholding, then you don't have to use any filter. Just throw away all those samples whose values are below your threshold (in this case, $100 \ \mathrm{pT}$).

In MATLAB, you could do something like this:

data(data < 100) = [];

If you want to keep the array size, then you might consider replacing those values with NaN or 0 instead of [].


czt of a short sinusoid I am using a Chirp-Z transform (aka. czt) to analyze a short length sample of a sinusoid. In scilab the code that I use to generate the signal is:

sampleRate = 4000000
points = 2^8
t = (0:(points-1))/sampleRate
y1 = sin(2*%pi*t*510000)

When I take the transform I get:

I expected to see only one non-zero magnitude sample. Obviously this is not the case. Why the sinc(x) waveform?

Here is my scilab code:

clear()
clc()
clf()

function a = angle(c)
// return the argument a in rad of complex number c
   a = atan(imag(c), real(c))
endfunction

function [Y, f] = chirpZ(y, sampleRate, f1, f2, m)
    w0 = 1
    phi = 2*%pi*(f2-f1)/(m*sampleRate)
    w = w0*exp(-%i*phi)

    a0 = 1
    theta = 2*%pi*f1/sampleRate
    a = a0*exp(%i*theta)

    Y = czt(y,m,w0,phi,a0,theta)

    f = (((0:(m-1))*(f2-f1)/m) + f1)
endfunction

sampleRate = 4000000
points = 2^8
t = (0:(points-1))/sampleRate
y1 = sin(2*%pi*t*510000)

f1 = 300000
f2 = 700000

m1 = points

[y2, f2] = chirpZ(y1, sampleRate, f1, f2, m1)
y3 = y2/points

subplot(211)
plot(f2,abs(y3),'g.')
subplot(212)
plot(f2,angle(y3),'b.')

[mag1,ndx1] = max(abs(y3))

subplot(211)
plot(f2(ndx1),abs(y3(ndx1)),'r.')
subplot(212)
plot(f2(ndx1),angle(y3(ndx1)),'r.')

You are not analyzing a sine wave, you are analyzing a windowed sine wave.

$$
y_1(t) = \sin(2\pi f_0t /f_s) w_{\rm rect}(n)
$$
where $f_0 = 510,000$, $f_s = 4,000,000$, and
$$
w_{\rm rect}(n) = \left\{  \begin{array}{cl}
1 & \mbox{where } 0 \le t \le (2^8 - 1)/fs\\
0 & \mbox{otherwise}
 \end{array}\right.
$$

One way to get closer to what you expect is to choose the duration so that your sine wave is "continuous" at each end.  That would mean setting

points = 40000; //2^8

for example.

If I do that, then I get

See this question and its answers for addressing a similar problem.

The reason I chose 40,000 points is to ensure that there is no discontinuity at the ends of the data. This is what happens when you plot [y2 y2] for 40,000 points vs 256 points as in your example.

Notice the discontinuity in the 256 case (second plot).

I agree that choosing a larger number of points is "interpolating the problem away".  There is no real way around that, however.  There will only be 255 zero points in the FFT of a signal of length 256, so if you're expecting to see zeroes everywhere except $f_0$, you're going to be disappointed.

Let's look at the analytic expression for the DFT of a windowed sine wave:
$$
X(\omega) = \sum_{n=0}^{N-1}\sin(2\pi \hat{f}_0 n)e^{-j2\pi nk/N} \\
= \sum_{n=0}^{N-1} \frac{1}{2}\left[e^{+j2\pi \hat{f}_0 n} - e^{-j2\pi \hat{f}_0 n} \right]e^{-j2\pi n\omega/N}\\
= \sum_{n=0}^{N-1} \left [ e^{+j2\pi n ( \hat{f}_0 - \omega )}  - e^{-j2\pi n ( \hat{f}_0 + \omega )}  \right]\\
= \frac{1 - e^{+j2\pi ( \hat{f}_0 - \omega ) N} }{1 - e^{+j2\pi ( \hat{f}_0 - \omega ) } } + \frac{1 - e^{-j2\pi ( \hat{f}_0 + \omega ) N} }{1 - e^{-j2\pi ( \hat{f}_0 + \omega ) } }\\
= \frac{e^{+j\pi ( \hat{f}_0 - \omega ) N}}{e^{+j\pi ( \hat{f}_0 - \omega ) }}
\frac{(e^{-j\pi ( \hat{f}_0 - \omega ) N} - e^{+j\pi ( \hat{f}_0 - \omega ) N})}
{e^{-j\pi ( \hat{f}_0 - \omega )} - e^{+j\pi ( \hat{f}_0 - \omega )} }
+
\frac{e^{-j\pi ( \hat{f}_0 + \omega ) N}}{e^{-j\pi ( \hat{f}_0 + \omega ) }}
\frac{(e^{+j\pi ( \hat{f}_0 + \omega ) N} - e^{-j\pi ( \hat{f}_0 + \omega ) N})}
{e^{+j\pi ( \hat{f}_0 + \omega )} - e^{-j\pi ( \hat{f}_0 + \omega )} }\\
= e^{+j\pi ( \hat{f}_0 - \omega ) (N-1)} \frac{\sin(\pi ( \hat{f}_0 - \omega ) N/2)}{\sin(\pi ( \hat{f}_0 - \omega )/2)}
+
e^{-j\pi ( \hat{f}_0 + \omega ) (N-1)} \frac{\sin(\pi ( \hat{f}_0 + \omega ) N/2)}{\sin(\pi ( \hat{f}_0 + \omega )/2)}
$$
which is where the sinc-like shape comes from.


Transfer function of a frequency shifting system There is a system which shifts frequencies of input by $-F_c$ such that:

  $$Y(S) = X(S).H(S)$$

But $X(S)$ has value zero from $0$ to $F_c$.

I am confused on how the product of $X(S)$ and $H(S)$ becomes a positive value in $Y(S)$ in that frequency range, for any $H(S)$?

How transfer function of the system $H(S)$ will look like in the frequency domain?

The system you're looking for cannot be described by a transfer function because it is a time-varying system. Only linear time-invariant (LTI) system can be fully characterized by a transfer function. However, there is no LTI system that can shift frequencies. The output of a (stable) LTI system can only have frequency components that are already present in the input signal.

Frequency shifts are usually achieved by modulation, i.e., by multiplying the input signal with a sinusoid or with a complex exponential. Such a system is linear but not time-invariant.


What input to use for an AR model of a vowel sound? I've recorded a 2-sec pronunciation of a vowel sound. The first 0.12 or so seconds of the signal are shown below.

Now, I've constructed an auto-regressive (AR) 8th-order model to compress this signal. (Actually, I'm just modelling 160 samples or 0.02 sec at a time.) The ar function in Matlab's System Identification Toolbox can estimate the parameters for an "optimal" spectrum fit.

My problem is choosing the stochastic input to the model filter. I suppose there's something better than white noise. The periodicity (14 periods per 0.02 seconds) leads me to think that an impulse train with the same period would be suitable.

If so, how would I choose the amplitude, and how would I find the periodicity? ACF and PSD estimations are quite noisy. Am I even on the right track?

A pitch estimator is commonly used to find the vocal periodicity.  Common pitch estimators include cepstrum/cepstral analysis, harmonic product spectrum, and composite algorithms, such as YAAPT.


How can I determine the complex envelope of a signal when knowing the analytic signal? Consider a signal $x(t) = \sum A_i \sin(2\pi(f_c+f_i) t + \theta_i),\, i=1\ldots N,$
and its analytic signal $Z_x(t) = x(t) + j\mathcal H[x(t)]$.

I want to describe the complex envelope (equivalent lowpass signal) as $x_{eq} (t) = r(t)\exp(j\Phi(t))$.

How do I know the amplitude $r(t)$ and the phase $\Phi(t)$ of the complex envelope?

As an example, let's consider: $x(t) = 1 sin(2\pi*11.5e9*t) + 2 sin(2\pi*11.6e9*t)$

In Matlab I would start the code with:

T = 1e-9;
f_s = 1e12;
t = [0 : 1/f_s : T]';
f_c = 11.5e9;

x_passband =  1*sin(2*pi*(f_c+0)*t) + 2*sin(2*pi*(f_c+0.1e9)*t);

x_analytic = hilbert(x_passband);

x_envelope = x_analytic .* exp(-1i*2*pif_ct);

The complex envelope should also be equal to:

amplitude = abs(x_analytic);

phase = angle(x_analytic);

x_envelope = amplitude * exp(j*phase);

With constant amplitudes and phases, this becomes a simple exercise:

First, go to the spectral domain
$$
X(f) = \mathcal{F}\{x(t)\}=\sum A_i \frac{1}{2j} (\delta(f-(f_c+f_i)-\delta(f+(f_c+f_i))\exp(j2\pi\tau_if)
$$
with
$$
\tau_i=\frac{\theta_i}{2\pi (f_c+f_i)}
$$
being phase offset translated to the equivalent time shift.

Now, you can shift this such that the carrier becomes DC:

$$X_d(f)=X(f-f_c) = \sum A_i \frac{1}{2j} (\delta(f-2f_c-f_i)-\delta(f+f_i))\exp(j2\pi\tau_i(f-f_c)$$

Applying a low-pass filter kicks out the high-frequency Dirac at $2f_c$. Then, transforming back to time, you get

$$
x_d(t) = -\sum A_i\frac{1}{2j}\exp(j2\pi f_i(t-\tau_i))\exp(-j2\pi\tau_if_c).
$$

This should be your complex envelope.

In case yo want to do this in Matlab, I would go for the following approach (which is completely different from the math above. It is the straight-forward implementation of shifting in frequency by multiplication in time):

T=1e-6;
fs=1e11;
t=[0:1/fs:T]';
fc=11.5e9;

B = 0.5e9;  % Signal bandwidth

carrier = exp(2j*pi*fc*t);

xpassband=1*sin(2*pi*(fc+0)*t)+2*sin(2*pi*(fc+0.1e9)*t);

xbaseband = carrier .* xpassband;

[b,a] = butter(6,B/(fs/2));
xbaseband_lowpass = filter(b, a, xbaseband);

subplot(2,1,1);
hold off;
plot(t, real(xpassband));
hold on;
plot(t, real(xbaseband_lowpass), 'r', 'LineWidth', 2);
plot(t, imag(xbaseband_lowpass), 'k', 'LineWidth', 2);
xlim([1e-7, 2e-7]);


$9\rm V$ Battery with $45^\circ$ phase How can you implement a $9\rm V$ battery with a phase of $45^\circ$?  (As a black box with a DC Voltage of $9\rm V$ and a phase of $45^\circ$)

Please preface your answer with spoiler notation by typing the following two characters first ">!"

OK, this is a slightly constructed situation, but as far as I can see, the following is the only thing that makes sense, kind of ...

  If $A$ is your DC voltage, then a phase of $45$ degrees means multiplying it with $e^{j\pi/4}$, i.e., you get $Ae^{j\pi/4}=\frac{A}{\sqrt{2}}(1+j)$. So scale the voltage with $1/\sqrt{2}$, and - apart from the ground wire - use two wires connected to '+' coming out of that box. Stick a tag with $j$ on it on one of the two.


Simulating noise in OFDM I'd like to simulate a certain $E_s/N_0$ in a baseband OFDM MATLAB modem to generate SER curves.

I randomize 50 QAM symbols with $N$ subcarriers each in frequency domain, add null frequencies (guard bands) to each symbol, perform IFFT and then add CP. Then I convolve the "time domain" symbols with a certain channel with taps of my choice and then I want to add noise, but I'm not sure about how to normalize it so that I obtain the correct SNR as generally reported in theoretical curves.

I get funky results with just MATLAB's agwn function with the measured flag.

In theory, all you'd need to do is add Gaussian noise of the appropriate variance. Without seeing your code and without knowing what "funky" means, I would check for:

when you do the convolution, scale the result by $\Delta t$, the time differential
make sure the FFT scaling that Matlab performs is not throwing you off

One thing you can do, which is simple but not very elegant, is to add the noise to each symbol after the FFT in the receiver, but before running any decision rule on them. You can simply add one sample of complex noise of variance $N_0$ to each symbol; since you control their average energy, you can get any SNR you want.


Correcting for Doppler shifts at the speed of sound for a 800 MHz signal I'm working with a system which is planned to reach velocities of about 340m/s. It's likely that we'll be using a frequency around 800 MHz for rx and tx. At these speeds, supposing a fixed source, the Doppler effect may lead to a frequency shift of about 890 kHz (supposing my math is correct).

This question did provide some information, but does not specify the higher frequency or higher vehicle velocity.

Supposing we don't have to worry about interference, how should I go about correcting this frequency shift? I was thinking about using an FLL or a PLL, but I'm not sure where whether these could support such important shifts.

frequency shift of about 890 kHz (supposing my math is correct)

$$\begin{align}
\Delta f_\text{Doppler} &= f_0 \frac vc\\
&=800\cdot 10^6\,\text{Hz}\frac{3.4\cdot10^2\frac{\text m}{\text s}}{3\cdot 10^8\frac{\text m}{\text s}}\\
&\approx 800\cdot 10^6\,\text{Hz}\cdot 1.13\cdot 10^{-6}\\
&\approx 906 \,\text{Hz}
\end{align}$$

So, your shift expectation is off by a factor of 1000.

  Supposing we don't have to worry about interference, how should I go about correcting this frequency shift? I was thinking about using an FLL or a PLL, but I'm not sure where whether these could support such important shifts.

You can design the bandwidth of your PLL or FLL to accomodate any shift. It just gets a bit tricky keeping other signals out of the observation.

Even your ~ MHz shift isn't all that much to deal with. Normal observers will simply have a  larger bandwidth than just necessary to observe your communication signal, anyways, and need

to have frequency correction anyways, since no two clocks / frequency generators are perfectly identical, physically (and if someone says "I've built this communication system. For now, I neglected synchronization needs", you know that they omitted the hard part)
When knowing what your communication signals look like, it's often not that hard to find them in the spectrum.
When you observe a communication partner for a while, you can easily write down a Doppler rate, i.e. the rate at which the Doppler shift changes. From that, you can predict future shifts and make your detection even easier (typical Kalman-Filter application)


Extract the count of particular macro block types(I, P and S) from the frames in the video How can I extract the count of particular macro block types(I, P and S) from the frames in the video (H.264 and MP4)?
Also, I would like to extract the prediction residual of each of the P and B frames (computed from its corresponding reference frame) in the video from the encoded version itself while the package is decoding it.

There are several tools to analyze video streams which which can provide slice and macroblock level information:

Intel Video Pro Analyzer
Tektronix Video Analyzer
H264bitstream
ffmpeg (page Analyzing Macroblock Types)
Avidemux

From How extract Intra-Prediction modes out of h264 stream, it seems that ffmpeg or ffprobe can do a lot of jobs.


How does a bandpass filter work? My question is about the action of a bandpass filter. First, I have a temporal signal. I convert it in a spectral signal via a FFT. Now, I want to apply a bandpass filter on that spectral signal. But my question is : how to do it ?

My temporal signal is a real signal so, I understand how to multiplicate by a coefficient. But my spectral signal is complex so, what do I have to do ? Multiplicate the real part and the imaginary part by the coefficient ?

Thank you for your answers,

Dr_Click

You're wanting the magnitude of your spectrum to decrease in certain areas (since you're trying to filter out those frequencies). So, yes, you could just multiply your real and imaginary components by your coefficient. For example, if you have a number $5 -5i$, the magnitude is $\sqrt{5^2 + (-5)^2} = 5\sqrt{2}$. If you want to decrease this magnitude by 50% you would just multiply your real and imaginary numbers by .5. This can be seen by $\sqrt{(5*.5)^2 + (-5*.5)^2} = 2.5\sqrt{2}$.
Note that if you're wanting to filter your signal in the frequency domain, and then transform it back to the time domain, you probably don't want to do a rectangular window, where you multiply the frequencies you want to filter out by 0, and everything else by 1. This would create some wide band noise in the time domain.


Magnitude and phase response and cut-off frequency of a moving average filter The frequency response of a typical moving average filter of length $N$ is given by $H(\omega)=\frac{1}{N}\frac{\sin(\omega N/2) e^{-j \omega ((N-1)/2)}}{\sin(\omega/2)}$. Firstly, isn't the cut off frequency of an average filter equal to 0 Hz (an average filter passes the dc value and filters the rest out)? By that notion i will get the value of $H(0)=1/N$ (using L hospitals rule).
But i don't think that's right. What am i missing here? Care to explain.

The cut-off frequency is whatever you define it to be. One standard definition would be the frequency $\omega_c$ at which the attenuation is $3$dB, i.e., $\big|H(\omega_c)\big|=\frac{1}{\sqrt{2}}$.

As far as we know there is no analytical formula for the exact computation of the $3$-dB cut-off frequency of a moving average filter. More than you ever might want to know about approximately computing the $3$-dB cut-off frequency can be found in the answers to this question.

Note that you made a mistake evaluating $H(0)$. The moving average filter satisfies $H(0)=1$ (and not $H(0)=1/N$). The real-valued amplitude function of a length $N$ moving average filter is

$$A(\omega)=\frac{\sin\left(\frac{N\omega}{2}\right)}{N\sin\left(\frac{\omega}{2}\right)}\tag{1}$$

Using $\sin(x)\approx x$ for $x\approx 0$ gives

$$A(\omega)\Big|_{\omega\approx 0}\approx\frac{\frac{N\omega}{2}}{N\frac{\omega}{2}}=1\tag{2}$$

A more straightforward way of seeing this is to note that

$$H(0)=\sum_nh[n]\tag{3}$$

where $h[n]$ are the filter coefficients. Since for a moving average filter we have $N$ filter coefficients with values $1/N$, the sum in $(3)$ equals $1$.


Emulating a Variable Delay I am trying to emulate a satellite link where the delay time between the satellite and ground station changes over time due to the satellite's motion. I plan do this by converting to the frequency domain via FFT, multiplying it by $e^{-2\pi \omega_{delay} }$, and converting back to the time domain via IFFT.

However, when I implement this in GNURadio, it fails to work as I expect.

I generate a signal, pass it through the delay, and then match filter the original signal with the delayed signal. I then change the delay with a GUI slider. When I use the frequency-domain-phase-shifting method that I implemented, the correlation peak does not move when I change the delay. However, if I substitute in gr-baz's variable delay block, the peak does move when I change the delay.

Why does my implementation of the shifting method not work? Is it because my delay block only works with a finite window of the signal?

Please see the answer to this question How to implement Polyphase filter?
on how to implement a polyphase interpolator, and then note that the same structure without actually commutating the output can be used as a discrete variable delay; each output is the input within the passband of the filter at a different delay in uniform discrete steps. By stepping through the output continuously as we do in the interpolator, we achieve a higher sampled version of the input. By staying at one output we get the input at a specific delay based on that filter, by changing to another output we get another delay. The design of the filter as given in the link will properly create the filters such that each is a different delay of the other (with the same ampltitude response) in uniform steps.

Continuing with the example shown of "Resampling" further shows how a variable delay device results:

Here is shown a IS95 CDMA waveform with 2 samples per symbol (as indicated by the red dots on this eye diagram). The analog waveform that is represented by these samples is shown in blue.

If we used the polyphase structure as an interpolator, we would get all the samples in between the 2 samples per symbol we have, evenly spaced, given by the number of filter banks in the polyphase structure.

Therefore it should be clear that we can select any filter output that is closest to the desired delay we want to achieve (in this case within one symbol duration).

As mentioned in the comments a Farrow Filter is another approach to implementing a variable fractional delay in an FIR filter.


Fourier transform of t*rect(t) In my previous post I asked for help for a Fourier transform of $$ t \text{rect} ( t- \frac{1}{2} ) $$ and I think I’ve understand the process. Now I’ve 2 another similar Fourier transform to do , I already solved both , but I don’t have the correct result. Can someone tell my if the 2 results I obtained are wrong ? Thank you

If $$ x(t) = t \text{rect}(t) $$ I obtained $$ X(f) = \frac {\text{sinc}(f) - \frac{1}{2} [e^{i \pi f } + e^{-i \pi f } ]}{i 2 \pi f } + \frac{1}{4} \delta(f) $$

And for $$ x(t) = \Big(t + \frac{1}{2}\Big) \text{rect}(t) $$ I obtained $$ X(f) = \frac{ \text{sinc}(f) - e^{- i \pi f } }{ i 2 \pi f } + \frac{1}{4} \delta (f) $$

Hint:

Derivative in frequency property of FT:

$$\frac{d}{d\omega}jF(\omega) \rightarrow tf(t)%$$

Time shift property of FT:

$$f(t-\tau) \rightarrow F(\omega)e^{-j\omega t_o}$$

FT of rect:  $F(\omega) = sinc(\omega/2)$

For practice if you are learning FT I recommend deriving the properties above to confirm you can get them from the FT directly using (practice makes perfect!):

$$F(\omega) = \int_{-\infty}^{\infty}f(t)e^{-j\omega t}dt$$


Confusion regarding plot of linear convolution vs fast convolution via FFT I am reading Proakis Digital signal processing using Matlab, 3rd edition
I am performing example 5.23 but i am not getting same plot/results as book although i am using almost same code
First i am attaching my code and output and then i will attach snapshot of book containing code and output
The issue with book code is that , i get error for "nu" variable and matlab says "NI" is undefined(i have underlined NI in book snapshot with red paint pencil)  , so i used "NL" in my code and "NL"=N*L ,but i am not getting output as book
My code is as follows
clc;clear all;close all;
conv_time = zeros(1,150); fft_time = zeros(1,150);
%
for L = 1:150
tc = 0; tf=0;
N = 2*L-1; nu = ceil(log10(N*L)/log10(2)); N = 2^nu;
for I=1:100
h = randn(1,L); x = rand(1,L);
t0 = clock; y1 = conv(h,x); t1=etime(clock,t0); tc = tc+t1;
t0 = clock; y2 = ifft(fft(h,N).*fft(x,N)); t2=etime(clock,t0);
tf = tf+t2;
end
%
conv_time(L)=tc/100; fft_time(L)=tf/100;
end
%
n = 1:150; subplot(1,1,1);
plot(n(25:150),conv_time(25:150),n(25:150),fft_time(25:150))

My output/plot snapshot is follow

[![enter image description here][2]][2]
Following is book snapshot
[2]: https://i.stack.imgur.com/1wPVM.jpg

The variable NI in the book is just a typo, it should be N instead (not N*L as in your code). Apart from that, remember that book was written about $25$ years ago, and the code was run on a $33$ MHz $486$ PC. So in order to see some effects on today's computers, you should crank up the value of L.
I've modified the code a bit (see below). Now the figure illustrates the expected result.

L = 200:500:10000;
K = length(L);
conv_time = zeros(1,K); fft_time = zeros(1,K);
Nav = 10;

for k = 1:K
    tc = 0; tf=0;
    Lk = L(k);
    N = 2*Lk-1; nu = ceil(log10(N)/log10(2)); N = 2^nu;
    for I=1:Nav
        h = randn(1,Lk); x = rand(1,Lk);
        t0 = clock; y1 = conv(h,x); t1=etime(clock,t0); tc = tc+t1;
        t0 = clock; y2 = ifft(fft(h,N).*fft(x,N)); t2=etime(clock,t0);
        tf = tf+t2;
    end
    conv_time(k)=tc/Nav; fft_time(k)=tf/Nav;
end

plot(L,conv_time,L,fft_time)


Estimating Damping Factor (Q) from Noisy Measurements I have a damped, tuned circuit and want to measure its Q factor. The hardware sends an impulse 'ping' and samples the output as it rings down.

Is there an efficient way to fit an equation of the form $A e^{t(i\omega  + b)} $ to a series of samples? I know $\omega$ up front, and $A$ isn't relevant. I am primarily interested in the value of $b$.

For my system the Q factors I want to measure are in the range 2-5, and the amplitude of the noise is about 10% of the signal.

My current solution uses a Goertzel filter to measure the energy at the frequency of interest for each of the successive cycles, then looks for the point where the magnitude drops by a quarter. With some interpolation between the cycles before and after the quarter point, this gives reasonable but not great results. It is especially bad for the low end of the Q range.

Our system has the impulse response (why we changed to imaginary instead of real? :) ):
$$h(t)=\mathbb{Im}(Ae(iwt)e(-bt))=Ae^{-bt}sin(wt)$$

With the following structure (ref.) (discretized):
$$H(s)=\frac{Kw}{(s+b)^2+w^2}$$

Under matlab, you only need to use procest:

procest(data,'P2U');

Which models the same structure under the form:
$$H(s)=\frac{K}{(1+2\zeta T_ps+T_p^2s^2)}$$


Inverse FFT of complex exponential with time delay - window leakage? I am working in MATLAB with a frequency domain signal that is a sum of complex exponentials, e.g.

$X(\omega) = \sum_i e^{-j\omega t_i} $ where $t_i > 0$ $ \forall $ $ i $

I expect the output to be a sum of delta functions, e.g.
$x(t) = F^{-1}[X(\omega)] \approx \sum_i \delta(t_i)$

I have noticed that when the time bin lines up exactly with the $t_i$, I get an excellent $\delta(t_i)$ function. However when the time bin is off slightly, there is significant energy spread throughout the rest of the time axis, as below:

Does anyone have a clearer idea why this happens? Or even better, a solution? (Perhaps a window that better treats these delta functions?, or a better way to do this IFFT?) Code used to generate image is below

Nfft = 2^12;
Tend = 20e-9;%20 ns of runtime
t = linspace(0,Tend,Nfft); %time
Ts = t(2) - t(1); %sampling rate, s
Fs = 1/Ts; %sampling rate, Hz
f = (-Fs/2):Fs/Nfft:(Fs/2 - Fs/Nfft); %frequency

%create complex exps. with delays 4ns, 6ns, 10ns
%4 ns exists as a time bin
ts = [4e-9 6e-9 10e-9];

%create exponentials
exps = zeros(length(ts),length(f));
for k = 1:length(ts)
    exps(k,:) = exp(-1j .* 2*pi.*f .*ts(k));
end
X_f = sum(exps); %create f-domain signal

x_t =real(ifft((X_f))); %ifft

plot(t.*1e9,20.*log10(x_t)); %plot
xlabel('Time, ns'); ylabel('dB');

Yes, this is leakage.

Your problem is in the discrete domain and your nomenclature is for the continuous domain.  Often times the two are conflated, as you have done, and this causes confusion.

The inverse DFT and the forward DFT are basically the same animal.  The only difference is the convention used in the sign of the exponent and the normalization factors you use.  So a window function can "shape" your sidelobes (leakage) somewhat, but can't eliminate them.

What is your larger purpose?


Compute ADSR: How? Does it make sense? I'm writing a Master degree thesis on Instrument Recognition using pyAudioAnalysis. I was wondering if results could be improved by also analyzing the ADSR shape of the samples.
All I find on the internet about ADSR is related to sound synthesis, not analysis.

Do you think it may be possible to compute ADSR, maybe manually by monitoring the energy over time and detecting when energy is kind of stable? Would it make sense?

... Instrument Recognition using pyAudioAnalysis. I was wondering if results could be improved by also analyzing the ADSR shape of the samples.

In a pattern recognition setting, additional features could lead to better classification results, provided that they can better "differentiate" the phenomenon being "measured".

Because of psychoacoustics, human beings are capable of distinguishing timbre from the first few milliseconds. As an extreme example, if you listen to a guitar and a similar piano tone, half way through they have started playing, it is very difficult to tell them appart. For example, here is Yngwie Malmsteen sounding like a flute, a violin, a guitar (or two) just by modulating the Attack via the guitar's volume (notice his right hand after the "without echo" part).

Now, for some piano tones and some guitar tones it might be difficult to tell them appart just by their spectral content because they are in similar ranges. But if you hit some bass tones on the piano, really hard, they will Sustain differently because the strings are now less of a string and more of a rod, which vibrates differently (than a taut string). In that case of course, the timbre will start having differences too.

So, ADSR based features would help in instrument recognition when the spectral characteristics of two or more instruments are near enough to be mixed up by a classifier using just spectral features. (For more information, please see this link).

How you are going to do this is another issue because it's not as easy as just running some metrics on some recordings. The instruments would have to be recorded "clean" in a room with minimal reflections and under controlled conditions. Certain instruments will sound differently for the same tone if it was to be played softly or loudly and the range of each instrument would have to be taken into account to decide on the Sampling Frequency. All of theses parameters will have an impact to your classification performance.

So, assuming a set of clean recordings, you could obtain the envelope of the waveform  to then try and infer the ADSR phases from that curve and use these figures in your classifier. For a brief overview of automatic ways of extracting the envelope please see this link while an example (using the analytic method of deriving the envelope) in Python is available from this link.

Hope this helps.


Detect when signal starts oscillating I have a signal which consists of oscillations riding on square waves as shown in the image below. I wish to detect the time at which oscillatory events occur (as shown by the red star).

Each square wave is a different height and the oscillatory sequences riding on each sqaure wave  are not identical. Here I have zoomed in on a small chunk of one oscillation.

A second order low-pass butterworth filter with cutoff at 500 Hz can be used to clean the data if this is a useful preprocessing step.

How can non-stereotyped oscillations be detected in a signal like this? Thanks!

As far as I can tell from the graph, the variance of the signal goes up substantially under "oscillation" conditions. So, monitor the variance over a rolling window. High variance indicates oscillation. To choose the window width, consider:

if the window is too short, the computed variance will be too noisy
if the window is too long, the monitor will be slow to respond

Also note that, given the information you've provided, there is no way to tell the difference between the step function signal itself and the beginning of an oscillation condition. You should therefore impose some auxiliary condition on the signal to avoid false positive detection. For a square wave signal, I'd suggest detecting extremely large time derivative. If the derivative exceeds some critical value, you probably have a square wave step, not an oscillation. So the criterion is:

oscillation if variance exceeds $V_{min}$ and max derivative is smaller than $D_{max}$.


Doppler Frequency Resolution for Radar Doppler frequency is calculated by the formula $f_{d} = 2*v_{r}/\lambda$, where $v_{r}$ is the radial velocity and $\lambda$ is the radar wavelength.

Normally FFT is used to detect the doppler frequency(A bank of FFT filters is used for Doppler frequency estimation).

Assume that sampling frequency is fixed as $4*f_{max}$, where $f_{max}$ is the maximum $f_{d}$ corresponding to the maximum radial velocity.

Now I have to determine the number of FFT points to be employed in the FFT calculation. How can I determine this.

In other words, which target parameter determines the Doppler frequency resolution.

Doppler resolution usually depends on the dwell time $T_{dwell}$, i.e., the time for which the radar is going to stare (look) at the target. If the pulse/sweep repetition interval is given by $T_{PRI}$ then $N = T_{dwell}/T_{PRI}$ is the number of pulses/sweeps in one dwell, which is also the limit of your FFT size. Your maximum Doppler corresponds to $N/2$th bin, therefore the Doppler resolution is $f_d/(N/2)$. Note that more pulses you are able to send to the target, better is the Doppler resolution.