SNR estimation: signal with unknown amplitude and Gaussian noise I would like to know can anyone suggest me techniques to estimate the SNR for a given noisy signal. I do not know the amplitude of my signal but I do know that noise is Gaussian.

I have tried to do estimation in frequency domain where I simply take the FFT of the signal and find the dominant peak and find the power associated with that as the signal power and noise power as the power in all other bins.

However, this technique seems to work well when I have sinusoidal signal and fails easily when I just have a noisy signal and do not know the nature of the signal since I do not know what peaks in FFT corresponds to my signal and how many of them.

A common technique relies on filter banks and robust statistics. The idea is to isolate some frequency subbands where the transformed signal is sparse, and the rest is (filtered) noise. From here, you can use a median estimator.

In the context of orthogonal wavelets (a form of dyadic filter banks), if the signal is sufficiently sampled, the highest frequency subband often satisfies the above requirements. Or one can choose a suitable union of subbands where the signal is considered sparse.

Then, if $c_i$ denote the subband coefficients, an estimator of the Gaussian noise variance $\sigma$ is:

$$ \hat{\sigma} = \frac{\textrm{median}|c_i|}{0.6745}\,.$$

The following picture displays the average noise estimation, with standard deviation, from several realizations of a Gaussian noise added to a line of an image.

This is described in Penalized threshold for wavelet 1-D or 2-D de-noising. More details are given in S. Mallat, A wavelet tour of signal processing, section 11.3 Thresholding Sparse Representations,
Noise Variance Estimation. The denominator factor is discussed in Rousseeuw and Croux, Alternatives to median absolute deviation, 1993.

A very basic avatar consists in differentiating the data with the basic $[1,-1]$ finite difference, and using the above estimator, further divided by $\sqrt{2}$ to account for its energy normalization.


Is the decimation process time invariant? I'm stuck at exercise 1.18 in Understanding Digital Signal Processing by Richard G. Lyons:

  There is an often-used process in DSP called decimation, and in that process we retain some samples of an $x(n)$ input sequence and discard other $x(n)$ samples. Decimation by a factor of two can be described algebraically by
  $$y(m) = x(2n) \tag{P1-3}$$
  where index $m = 0,1,2,3,. . .$ The decimation defined by Eq.
  (P1-3) means that $y(m)$ is equal to alternate samples (every other
  sample) of $x(n)$. For example: $y(0) = x(0)$, $y(1) = x(2)$, $y(2) = x(4)$,
  $y(3) = x(6)$, $. . .$ and so on. Here is the question: Is that decimation
  process time invariant?

  Illustrate your answer by decimating a simple sinusoidal $x(n)$ time-domain sequence by a factor of two to obtain
  $y(m)$. Next,create a shifted-by-one-sample version of $x(n)$ and call it
  $x_{shift}(n)$. That new sequence is defined by (P1-4):

  $$x_{shift}(n) = x(n+1) \tag{P1-4}$$
  Finally, decimate $x_{shift}(n)$ according to Eq. (P1-3) to obtain $y_{shift}(m)$. The decimation process is
  time invariant if $y_{shift}(m)$ is equal to a time-shifted version of $y(m)$. That is, decimation is time invariant if
  $$y_{shift}(m) = y(m+1)$$

I'm somewhat confused by the notation $y(m) = x(2n)$--according to the definition of the decimation, I would hve written $y(n) = x(2n)$. Anyway, if I try to resolve the problem as it is stated:

"Illustrate your answer by decimating a simple sinusoidal () time-domain sequence by a factor of two to obtain ()"

Let say $x(n)=sin(\omega n) \text{ with } \omega=2\pi f_0 t_s$
$\implies y(m) = x(2n)=sin(2\omega n)$
"create a shifted-by-one-sample version of $x(n)$ and call it $x_{shift}(n)$"

$x_{shift}(n) = x(n+1) = sin(\omega n+\omega)$
"decimate $x_{shift}(n)$ according to Eq. (P1-3) to obtain $y_{shift}(m)$"

$y_{shift}(m) = x_{shift}(2n) = x(2n+1) = sin(2\omega n +\omega)$
"That is, decimation is time invariant if $y_{shift}(m) = y(m+1)$"

$y(m+1) = x(\text{???})$

Eq (P1-3) doesn't really give a relation between $m$ and $n$. So what is $m+1$? According to the decimation definition given in the exercice, I would be tempted to say that $m+1=2n+2$. But that leads to $y(m+1)=sin(2\omega n +
2\omega) \neq y_{shift}(m)$

Could you pinpoint my mistake or clarify my misunderstanding in this exercise?

Sample rate conversion systems (expansion or decimation) are time-varying operators.

For your example system of decimation by $M$ :

$$ y[n] = T\{x[n]\} = x[Mn] $$

you can easily see that results of shift in the input and output are not the same; i.e.,

$$ y_1[n] = T\{ x[n-d] \} = x[Mn - d] $$

and

$$ y_d[n] = y[n-d] = x[M(n - d)] \neq y_1[n]$$


Python MNE - reading EEG data from array I have EEG data that comes in the form of a 3D numpy array (epoch * channel * timepoint). timepoint is a 256 element array containing each sampled timepoint (1s total, at 256Hz). epoch is an experimental trial.

I'm trying to import the numpy array into a form Python-MNE (http://martinos.org/mne/stable/mne-python.html) understands, but I'm having some trouble

First, I'm not sure if I should be importing this raw data as a RawArray or an EpochsArray. I tried the latter with this:

ch_names = list containing my 64 eeg channel names
allData = 3d numpy array as described above

info = mne.create_info(ch_names, 256, ch_types='eeg')

event_id = 1

#I got this from a tutorial but really unsure what it does and I think this may be the problem
events = np.array([200, event_id])  #I got this from a tutorial but really unsure what it does and I think this may be the problem

raw = mne.EpochsArray(allData, info, events=events)

picks = mne.pick_types(info, meg=False, eeg=True, misc=False)

raw.plot(picks=picks, show=True, block=True)

When I run this I get an index error: "too many indices for array"

Ultimately I want to do some STFT and CSP analysis on the data, but right now I'm in need of some help with the initial restructuring and importing into MNE.

Whats the correct way to import this numpy data that would make it easiest to complete my intended analyses?

mne.EpochsArray is for 3-D data (epochs * channels * times). mne.RawArray is for 2-D data. Use EpochsArray.

events is an n * 3 integer array. The 3 columns are: time (in sampling points), length (you can put a dummy here - it is almost never checked - but you still need 3 columns), value (e.g. condition). You fed it a 1 * 2 array.

Try it:

import numpy as np
from mne import create_info, EpochsArray
n_epochs = 100
channels = ["a", "b", "c", "d"]
n_channels = len(channels)
events = np.array([np.arange(n_epochs), np.ones(n_epochs), np.ones(n_epochs)]).T.astype(int)
d = np.random.random((n_epochs, n_channels, 256))

info = create_info(channels, 256, "eeg")
epochs = EpochsArray(d, info, events)

epochs.plot(show=True, block=True)


FIR filter design by the Fourier transform method I am having some problems understanding how the Fourier transform method is used to determine the FIR filter.

As far as i have understood, you start by using the ideal impulse reponse for the specific filter design you want to have.  Those are given in this table

As the the ideal impulse reponse contains infinite impulses, we truncate the infinite to finite a number by $m$, (which is confusing how this is done).

But my question is how do predetermine the order of the filter, without having to try randomly, using this method?

Some papers provide rule-of-thumb formula for the filter length (and coefficient quantization), for instance:

Accurate estimation of minimum filter length for optimum FIR digital filters, 2000, Ichige et al.
On computational complexity in digital filters, 1981, Bellanger
FIR Digital Filter Design Techniques Using Weighted Chebyshev Approximation, 1975, Rabiner et al.

The one I generally  use in lectures is borrowed from the second paper, for a low-pass filter with the following design:

Then, the estimated filter length $N$ (apparently $2M+1$ for you) is, with sampling frequency $f_s$:

$$ N = \frac{2 f_s}{3(f_2-f_1)}\log_{10} \left(\frac{1}{10\delta_1\delta_2}\right)$$


Getting a weired behavior when windowing an audio signal I am processing audio signals using Essentia, an open-source library and tools for audio and music analysis, description and synthesis.
Here is a waveform corresponding to a frame extracted from an audio signal. The plot is obtained with the following instruction:

plt.plot(essentia.array(audio_pool['frames.actual-frames'][22]))

... and here is the previous waveform, windowed. The plot is obtained with the conjunction of these two instructions:

w = std.Windowing(type='hamming',
                  size=len(audio_pool['frames.actual-frames'][22]),
                  normalized=False)
plt.plot(w(essentia.array(audio_pool['frames.actual-frames'][22])))

We can see that the resulting waveform is somehow inverted, i.e. the beginning of the original signal occurring in the middle.

I tried other window types (rect, hann, etc.) as well as w/ and w/o zero-padding the original signal, I get the same behavior.

It turns out that the windowing function that is used performs by default "Zero-Phase Zero-Padding" on the input signal.
You can disable this feature by setting the zeroPhase parameter to False as follows:

w = std.windowing(type='hamming',
                  size=len(audio_pool['frames.actual-frames'][22]),
                  normalized=False,
                  zeroPhase=False)


Noise removal in medical segmented image Can anyone suggest methods for removing the noise (indicated inside the red square) from every where in the following image, while keeping the white lines?

One solution I have found is as follows:

Thresholded on grayscale value.
Remove objects on the basis of size.
Some more morphological operations.


downsample signal to a specific number of samples in Python The Python Scipy library provides several functions to downsample signals, but they all have limitations:

The resample function is based on Fourier method, which means it assumes periodic signals.
The resample_poly function assumes "values beyond the boundary of the signal to be zero."
The decimate function can only downsample signals by some integer factors, not to  specific number of samples.
The upfirdn function requires FIR filter coefficients as inputs, which I am not sure how to get.

I have a signal that's NOT periodic and the values beyound the boundary of the signal are NOT zeros, and I want to downsample the singal from 611 samples to 100 samples.

Is there a simple way to do that in Python?

Thanks.

EDIT: The reason why I need to downsample signals to a fixed number of samples is that I have signals with different length, and I want to use those signals to train a 1D CNN for classification/regression tasks, which requires inputs of fixed size; I prefer not to use RNN/LSTM, although those models can accept inputs with various sizes.

I am doing something similar to your application with 1D CNN. I think scipy.resample_poly is the most versatile function since it allows both upsampling, downsampling, or a combination of both.

Also your comment about "values beyond the boundary of the signal to be zero" can be solved by using the option "line" in padtype, as shown in the function documentation, in order to eliminate spurious effects at the head/tail of the signal.

In general if the integer number of resampling factor is still an issue to you (but it shouldn't as you can play with the length of the signal) you can play with a combination of moving average followed by interpolation on the grid.


get convex hull of isolated shapes with openCV/C++ I have segmented CT images, something like the following one:

I would like to get a convex hull that contains all the white regions.
Something like this circle, but more fitting, and of course not necessarily a cicle.

I've already tried this example https://stackoverflow.com/questions/11813352/finding-the-convex-hull-of-an-object-in-opencv but I've only got this output:

Try dilation, a technique used to enlargen light spots over darker spots in an image.

Start by dilating the image such that many of the enlargened white blobs will overlap and form a complete connected contour, as such:

Now, when OpenCV finds contours from this large figure to make convex hulls, it will find a complete circle rather than loose seperate blobs.

From here, you can find the largest contour, which is the (outer) largest light-purple contour.

Here is the code for this (which currently computes the contours, but not the largest area - this is just to demonstrate the dilation):

#include "opencv2/highgui/highgui.hpp"
#include "opencv2/imgproc/imgproc.hpp"
#include <iostream>
#include <stdio.h>
#include <stdlib.h>

using namespace cv;
using namespace std;

int main( int argc, char** argv )
{
 Mat src; Mat src_gray;
 src = imread( "convex.jpg", 1 );
 //resize(src, src, Size(640,480), 0, 0, INTER_CUBIC);
 cvtColor( src, src_gray, CV_BGR2GRAY );

 namedWindow( "Source", CV_WINDOW_AUTOSIZE );

 //Dilate
 int erosion_size =10;
  Mat element = getStructuringElement(cv::MORPH_ELLIPSE,
        cv::Size(2 * erosion_size + 1, 2 * erosion_size + 1),
        cv::Point(erosion_size, erosion_size) );

 dilate(src_gray,src_gray,element);

 imshow( "Source", src_gray );

 // Convex Hull implementation
 Mat src_copy = src_gray.clone();
 Mat threshold_output;
 vector<vector<Point> > contours;
 vector<Vec4i> hierarchy;

 // Find contours
 threshold( src_gray, threshold_output, 200, 255, THRESH_BINARY );
 findContours( threshold_output, contours, hierarchy, CV_RETR_TREE, CV_CHAIN_APPROX_SIMPLE, Point(0, 0) );

 // Find the convex hull object for each contour
 vector<vector<Point> >hull( contours.size() );
 for( int i = 0; i < contours.size(); i++ )
 {  convexHull( Mat(contours[i]), hull[i], false ); }

 // Draw contours + hull results
 RNG rng;
 Mat drawing = Mat::zeros( threshold_output.size(), CV_8UC3 );
 for( int i = 0; i< contours.size(); i++ )
 {
  Scalar color = Scalar( rng.uniform(0, 255), rng.uniform(0,255), rng.uniform(0,255) );
  //drawContours( drawing, contours, i, color, 1, 8, vector<Vec4i>(), 0, Point() );
  drawContours( drawing, hull, i, color, 1, 8, vector<Vec4i>(), 0, Point() );
 }

 // Show in a window
 namedWindow( "Hull demo", CV_WINDOW_AUTOSIZE );
 imshow( "Hull demo", drawing );

 waitKey(0);
 return(0);
}


Upsample / Downsample a Signal from 1920.93Hz to 1920Hz I've just learned that the A/D converter I am using is not sampling at my desired sampling rate.  I request 20 seconds of data at 1920Hz and I receive 38400 points of data that were actually sampled at 1920.93Hz.

So for example, if I sample a 30Hz sine wave and chop the signal up into 300 sections, the phase of the sine wave in section 1 differs from the sine wave in section 300.  If my sampling rate was 1920Hz, as expected, the phase of each section would line up perfectly.

Now, to correct for this, is there any way of processing the signal such that it appears as if it was sampled at the correct sampling rate through interpolation or other means?

Assuming the signal was reasonably band-limited before sampling, you can just use a high quality interpolation, such as a Sinc interpolation kernel, to calculate new equally spaced points at any lower density.

For small frequency deltas, even a long interpolation kernel or high degree polynomial interpolation can be a lot more efficient than naively upsampling and then downsampling.


Discrete signal vs digital signal? Is there any difference between the two?

When we say discrete signal, does it only imply that the x-axis or time is discrete? Or it can also mean discrete amplitudes on the y-axis?

Or is the signal is called digital only when both of its axes are discrete?

For most in signal processing, a digital signal (source: wikipedia)

  is a representation of a physical signal that is a sampled and
  quantized

So it represents  data as a finite sequence of  discrete values. Often, digital data is embedded into a specific format to record/store values (for instance 16-bit words), possibly with a specific coding of symbols (think about sign conventions, or floating point standards), to end up in a concrete bit-wise or symbol-wise representation.

A discrete signal is a more abstract thing. It is an abbreviation for a discrete-time signal. It may possess a countable number of samples, and their values are not restricted to be discrete. Most often, the time interval between samples is constant (or uniform). To me, lacunary signals that are uniformly sampled, but for which some samples are missing,  are still considered as discrete.


HMM's - Possibility? I have a basic understanding of HMM's and as a final year project at school wanted to expand my knowledge further.

My basic idea was to create an algorithm that determines if someone is saying either Yes or No. As interesting as this sounds, it's relatively easy and I want something that challenges me even more and wanted to know if my idea was possible.

My idea:

The algorithm will be able to answer certain questions (given samples).. E.g.

Q: Is this person saying either Yes no No?

A: Yes

Q: Is this person male or female?

A: Male

Q: Is this animal a bat?

A: No

Obviously, my data-set will be massive, but is this possible using HMM's? So, in theory, instead of handling whether someone is saying yes or no, the algorithm becomes universal and I can determine many many things.

I hope someone can help me.

HMM are useful for sequence modeling and classification - problems for which your observations unfold on a 1-D axis in time or space. Hence their usefulness for speech recognition, because a word is a sequence of heterogeneous states corresponding to its various phones. But the problem of recognizing whether a speaker is male or female doesn't really have this sequential aspect, so HMM would be indeed an odd approach for it. One good way of finding out whether HMM are a good fit for a problem is the following: would you still be able to perform the classification if the data is shuffled? If you take a speech recording, break it into segments and shuffle them, you won't be able to recognize the words spoken, but you will still be able to recognize whether the speaker is male or female, and indeed identify this speaker.

I am not sure about the bat question - are you talking about recognizing animals from the sound they emit? Or else, from what kind of data? Does this data have a sequential structure?

It seems to me that what you want to do is, from a broad point of view, machine learning or pattern classification. HMM is a very specific solution to machine learning problems (those related to classifying processes unfolding along a 1-D axis), but this is by far not the only tool! Maybe you should learn more about the other approaches (bayesian methods, support vector machines, neural networks, decision trees...), and the specifics of each data types (features for image, audio, text...).


What is the moment when all oscillators aligned to make a jump called? Say we have a square wave and its Fourier series. When the wave doesn't jump, its oscillators aren't aligned:

But if they are aligned, the wave will jump:

What is this moment called? They might not be necessary in the same angle, and I don't think "resonating together" is correct, since they have different frequencies. What is it called when all oscillators aligned to make a jump?

Meta discussion: Why are questions asking for terms off-topic?

I think the word/phrase you are looking for is "discontinuity" / "discontinuity in the source signal".  Although a jump is always a discontinuity, a discontinuity is not always a jump.  When it isn't a jump, it is known as a pluggable discontinuity, i.e. you can define a value to plug a hole.  The sinc function at zero is an example of this.  So to be more specific you would say "an unpluggable discontinuity".

Whether a series of continuous functions could converge to a discontinuous function was a question that led to some very serious, and non-intuitive, math.

As far as I know, there is no word to describe the state of the components at this moment.  The closest word I can think of is a secondary meaning of "concordance" = "agreement".

Hope this helps.

Ced


zero padding shifting signal peak? I am trying to find the fundamental frequency of a low-frequency signal. I need an estimate that is precise to .01 Hz based on only a few cycles, so I'm trying to code up a fft in Python. The signal itself will be almost a pure tone, so I'm testing it out on a simple simulated wave.

I'm not seeing the results I expect, which may just be a dumb coding mistake. I'm new enough to DSP, though, that it's also possible I'm fundamentally misunderstanding what the transform is doing. Any help would be very much appreciated!

from math import *
import numpy as np
import matplotlib.pyplot as plt

I start by setting some system parameters

fs = 10000 # sampling rate (Hz)
T = 0.1 # length of collection (s)
assert fs*T % 1 == 0 # make sure we have an integer number of samples
windowlength = int(fs*T) # total number of samples

Then I create a sampled cosine wave for testing. Its period isn't an integer multiple of my window, so I'll need to zero-pad

f0 = 20.1 # fundamental frequency in Hz
checkwave = [0]*windowlength
for x in range(windowlength):
    checkwave[x] = cos(2*pi*f0*x/fs)

Plotting my wave in the time domain looks fine

t = np.linspace(0,T,int(fs*T)) # time domain x axis
plt.figure()
plt.plot(t,checkwave)

I then carry out the zero-padding and take the dft with numpy's rfft module

zp = 10 # factor to use in zero padding
dft = np.fft.rfft(checkwave,zp*len(checkwave))
fmaxindex = np.abs(dft).argmax() # find the peak frequency's index
fmax = fs*fmaxindex/(zp*len(checkwave)) # calculate the peak frequency
print(fmaxindex) # Higher than I expected
print(fmax)

But my plot in the frequency domain is off. I thought I should be looking at a sampled sinc wave centered on 20.1Hz, and since I'm sampling well above the Nyquist limit the centering should be perfect, right? But, if that were the case, 20Hz should be the largest bin (translates to bin 20). Instead, the largest bin is 21Hz. Plotting, we can confirm that's the case.

f = np.linspace(0,fs/2,int(fs*T*zp/2)+1)# frequency domain x axis
freqspan = 4 # number of samples to plot around the peak, so we get some detail
plt.figure()
plt.plot(f[fmaxindex-freqspan:fmaxindex+freqspan+1],abs(dft[fmaxindex-freqspan:fmaxindex+freqspan+1]))

My first thought was I made an off-by-one indexing error, but when I reproduce my last two blocks of code with an extra factor of ten on the zero padding, I see that things get even worse.

zp = 100 # factor to use in zero padding
dft = np.fft.rfft(checkwave,zp*len(checkwave))
fmaxindex = np.abs(dft).argmax() # find the peak frequency's index
fmax = fs*fmaxindex/(zp*len(checkwave)) # calculate the peak frequency
print(fmaxindex) # Higher than I expected
print(fmax)

Plotting in the frequency domain:

f = np.linspace(0,fs/2,int(fs*T*zp/2)+1)# frequency domain x axis
freqspan = 6 # number of samples to plot around the peak, so we get some detail
plt.figure()
plt.plot(f[fmaxindex-freqspan:fmaxindex+freqspan+1],abs(dft[fmaxindex-freqspan:fmaxindex+freqspan+1]))

Now the peak frequency is 20.8 (bin 208), about 6 bins away from where the sinc should be putting it.

Thanks in advance for any thoughts.

The FFT of a strictly real input signal is conjugate symmetric.  What that means is that, for a single sinusoidal frequency, the FFT result includes not just one positive frequency peak, but also a another peak at the negative frequency.  For signals that are not strictly integer periodic in the FFT length, the peak can be a very broad Sinc after zero-padding.  For extremely low numbers of signal periods (and for frequencies that are very near Fs/2), the negative frequency peak is close to and thus can interfere with the positive frequency peak and distort its exact location.

So, for very small numbers of periods that are not known to be strictly integer periodic in aperture, it may be better to use a parametric estimator of some sort, rather than an FFT, even windowed or zero-padded.


Why isn't this AM Demod GNU Radio flow graph working with different sampling rates? I am trying to demodulate an AM signal in GNU Radio, and I am having problems with the sampling rates.  I will start off with an example where the flow graph was working. The sampling rate is 256kHz, and it works just fine:

Next, I tried a sampling rate of 192kHz, and then it didn't work.  Is there a problem with this flow graph?

Next, I modified the 2nd flow graph (changed the sampling rate, but didn't need to), and it worked almost perfectly, other than that the voice sounded very slightly too low pitch/slow.  This sampling rate suffered from the same problem as that of figure 2, however I tweaked the rational resampler block until the voices sounded normal.

As some background, I am running an RTL SDR (DVB-T) dongle, on an Odroid U3 (ARM architecture).  I have a "Ham it up" upconverter, so that the AM (810Kkz) is shifted up by 125MHz, into the VHF band that the dongle can actually receive.

UPDATE:  There were two issues with the flow graph:

The filter cutoff frequency needed to be calculated using normalized frequency (see MBaz's detailed answer).
The R820T chip does indeed, as MBaz ascertained in the comment section, have only a limited set of sampling rates that it supports.

In order for the flow graph to output natural sounding audio, it must use one of the acceptable sampling rates, and also have the proper filter setting.  The 2nd and 3rd flow graphs didn't work because they had the wrong sample rates.

In order to make the filter calculations easer, you should specify its sampling rate equal to the source's. Assume that:

$D$ is the filter's decimation factor
$f_s$ is the source sampling rate
$f_c$ is the filter's specified cutoff frequency
$f_f$ is the filter's specified sampling rate

The actual cutoff frequency from the point of view of the filter's output, $f_a$, is calculated as: $$f_a=\frac{2Df_c}{f_f}\frac{f_s}{2D}=\frac{f_cf_s}{f_f}.$$ This is because the filter doesn't care what the actual sampling rate is; what it needs to know is the cutoff frequency as a fraction of the normalized Nyquist frequency, which is $2Df_c/f_f$ after accounting for the decimation.

In your first flowgraph, the actual cutoff frequency is $5\cdot256/64=20$ kHz, not 5 kHz (the same problem occurs in the second flowgraph; in the third, the cutoff is 25 kHz). If you specify the filter sampling rate equal to the source rate ($f_s/f_f=1$), then your signal will actually be filtered at 5 kHz.

You also need to make sure that the product of the source sampling rate, times all your interpolation, divided by all your decimation, is equal to the sink sampling rate. In the third flowgraph, the sampling rate into the sink is 16 kHz, not 48 kHz.

Having said that, I think there's another problem in your flowgraph. I assume that you're demodulating AM DSB-LC. If the RTL-SDR device is centered exactly at the carrier frequency, then you only need the low-pass filter; no complex to mag conversion is necessary. You only need to keep the real part of the signal (send the imaginary part to a null sink -- you can do this right after the RTL-SDR and save some computations).

You only need to do envelope detection when the AM signal is not in baseband. In that case, you calculate its absolute value and then low-pass filter it. You also need a final high-pass filter to eliminate the DC bias present at the LPF output.


How to write e^(-1000|t|) in matlab title? I wanted to write e^(-1000|t|) in matlab title of a plot. When I entered this and ran it actually, it showed ( as a superscript and every thing else in the line in regular font.  Now I wanted the whole thing -1000|t| to appear as a superscript.  I tried bracketing etc and finally I achieved what I wanted by inserting ^ before every letter. Is there a better way?

You can write e^(-1000*|t|) in Matlab in colored form if you want as

title('\ite^{-1000 |t|}','Color','b')
Here,"b" is for blue color ,"r" for red,"g" for greeen.

For subscript use _

For superscript use ^

For exponential symbol use \ite ,

For lambda symbol use \lambda

For alpha use \alpha.


How to interpret the effect of different windows in short time fourier transform? There are many categories of windows, e.g., rectangular, Gaussian, and triangular. What are their effects on STFT?

A widnow $w[n]$ truncates and weights (tapers) an input signal $x[n]$, to produce $v[n] = x[n]. w[n]$., for subsequent spectral analysis of $x[n]$. A windows's effect on the input signal's true spectrum $X(e^{j\omega})$ is described by a convolution of $X(e^{j\omega})$ with $W(e^{j\omega})$ (window's Fourier transform);
$$V(e^{j\omega}) = \frac{1}{2\pi} \int_{-\pi}^{\pi} W(e^{j\theta})X(e^{j (\omega -\theta)}) d\theta $$
due to this convolution, two main effects are observed on $X(e^{j\omega})$.
1 - smoothing (smearing) of $X(e^{j\omega})$ due to the main lobe width of $W(e^{j\omega})$,  which results in a loss of spectral resolution.
2- spectral leakage due to peak side lobe of $W(e^{j\omega})$, which results in a loss of weak components shadowed by nearby stronger ones.
Main lobe width of any window type is primarily determined by its length. Increasing the length of any window will therefore decrease its main lobe width (hence increase its spectral resolution capability)
A rectangular window has the narrowest main lobe width and the highest peak side lobe, compared to all other windows. And the remainig window types perform a tradeoff between mainlobe width and the peak side lobes.
Peak side lobe is primarily determined by windows's shape. So by changing the shape (type) of the window you adjust it's leakage amount.


How does filter works in noise reduction? On page 22 of the slide, the 1D array is [1, 1, 1, 1, 1], but the dots are uneven. Then why their values are the same?
What's the value of the point that the arrow is pointing at? Though it gives Moving average in 1D: [1, 1, 1, 1, 1]/5, what exactly is that value?

The illustration is showing your data (labelled as "original") and the result of convolving your moving average filter with your data (labelled as "smoothed"). The filter/kernel is [1 1 1 1 1]/5 and gives each of the five samples equal weight. Effectively, it takes the average of the five samples it overlaps with during each step of the convolution process. It can also be written as [0.2 0.2 0.2 0.2 0.2] or [1/5 1/5 1/5 1/5 1/5]. The values in this filter are the same, but you can also choose different values depending on what you wish to accomplish. In the next slide, they give the filter [1 4 6 4 1]/16. This filter gives more weight or importance to the samples that are near the center of the filter.
I hope this helps.


Generating a time domain signal using the fundamental frequencies I currently have the power spectrum (in dB) of a signal through which I obtained the fundamental frequencies and amplitudes. My question is: how to generate a time domain signal using these amplitudes and frequencies (using C++ for example)?

No way. You cannot recover a time domain signal from its fundamental frequency. Although, it is not clear what does a plural means here, "the fundamental frequencies and amplitudes", so maybe I am not quite understanding you.
Also, the power spectral density (PSD) alone, whether in dB or else, is not sufficient to recover a time domain signal: you need an energy spectral density for this deed. Citing a Wikipedia article:

The above definition of energy spectral density is suitable for
transients (pulse-like signals) whose energy is concentrated around
one time window; then the Fourier transforms of the signals generally
exist.

The Fourier transform of the signal is what you need in order to recover a time domain signal, because the PSD lacks the phase information. Moreover, a periodogram, if it is what you hint by mentioning "fundamental frequency", does embrace averaging, either in time or frequency, which makes the periodogram still less suitable for a time domain signal recovery.
Although it seems that your question is about PSD, you may be interested in retrieving phase information from energy spectral density. There is a number of techniques, the classical one is use of the Kramers Kronig relations for analytical signal. See also a review of three techniques (including KK) used in optical measurements, in Deep learning as phase retrieval tool for CARS spectra. Section 2.2 of this article presents the formula for deducing the phase from the squared modulus (formula 15), with stipulations made about a noise and discrete nature of measurements.
A word about implementation: in order to write a C++ implementation, you need a thorough understanding of signal processing algorithms. Searching github.com for open source implementations of phase retrieval techniques, I've found a Python implementation of phase retrieval for Raman spectroscopy, with Jupyter notebook, explaining the implementation details, see CRIKit2: Hyperspectral imaging toolkit. You may browse this presentation, at least out of curiosity.


Direct Correlation (DC) Time Delay Estimation: Variance Keeps Decreasing for Increasing SNR? I'm trying to reproduce the results from this paper "Discrete time techniques for Time Delay Estimation" doi:10.1109/78.193195, for both the correlation and Least Squares.

I've generated a random (statistically gaussian distributed) signal and filtered it with a gaussian low-pass, to make it band-limited. Then I added increasing noise to the signals.

I've tried for 0, subsample and 1< sample delays. In all cases, the variance of my DC keeps decreasing as I increase SNR (actually, the behavior closely matches that of the LS or average square difference function, as they call it).

Is there any reason as to why I don't see that asymptotic behavior predicted by the paper?
Does it have to do with the signal properties or the way I'm processing in some way?

Also, if possible, can someone please make sense to me of why that limit exists for the correlation, conceptually? I can't seem to understand from the paper.

Here is a simplified version of the code that shows the problem.

Any help would be greatly appreciated.

This is a plot from the article, showing the effect I can't reproduce.

EDIT/UPDATE: Also tried for different reference correlation/least squares window sizes, as represented on the code. Essentially, I considered both the case where I also compare the signal with the surrounding data of the time window under analysis, and for the case where I ignore that data.

My reasoning is that a time delay in a continuous signal will bring uncorrelated samples into the signal. If I consider the surrounding window the number of uncorrelated samples will be constant, and not dependant on the delay (as long as my delayed window is still contained within the larger reference window.). This makes least squares more intuitive, also, since I always divide by the same number of points.
Either way, I get the same result.

I can reproduce a plateau value of about $0.0013$, compatible with the figure, as the error due to parabolic interpolation (2nd degree polynomial interpolation) around the peak of noise-free cross-correlation.

Without loss of generality, we assume $T=1$ and a uniform distribution of the delay $d$ between $-0.5$ and $0.5.$ Also we assume a perfectly white critically sampled signal and its delayed copy. Their cross-correlation is the critical samples of a time-shifted sinc function. We do a parabolic fit on three successive samples around the sample with the largest value. Calculating the peak of the parabola will estimate the delay as (calculations omitted):

$$\hat d = \frac{d}{2-4d^2}\tag{1}$$

Figure 1. Estimate $\hat d$ and true delay $d.$

The estimate is unbiased, as the expected value of its error is zero due to the symmetry seen in Fig. 1:

$$\operatorname{E}\left[\hat d - d\right] = 0$$

The variance of the error of the estimate will be:

$$\operatorname{Var}[\hat d - d] = \int_{-0.5}^{0.5}\left(\frac{d}{2-4d^2} -
 d - \operatorname{E}[\hat d - d]\right)^2dd \approx 0.0072$$

If we treat the left and right half-ranges of $d$ separately, the estimate is actually biased. The error can be reduced by subtracting from the estimate the expected value of the error conditional on the sign of $d$:

$$E\left[\hat d - d\,|\,-0.5 \le d < 0\right] = 2\int_{-0.5}^0\left(\frac{d}{2 -
 4 d^2} - d\right)dd = -\frac{\ln(2)}{4} + \frac{1}{4}\approx 0.07671320485\\
E\left[\hat d - d\,|\,0 \le d < 0.5\right] = \frac{\ln(2)}{4} - \frac{1}{4}$$

$$\hat{d_2} = \begin{cases}\hat d + \frac{\ln(2)}{4} - \frac{1}{4}&\text{if }-0.5 \le d < 0,\\
\hat d - \frac{\ln(2)}{4} + \frac{1}{4}&\text{if }0 \le d < 0.5.\end{cases}$$

Figure 2. Estimate $\hat{d_2}$ and true delay $d.$

$$\operatorname{E}\left[\hat {d_2} - d\right] = 0$$

$$\operatorname{Var}[\hat {d_2} - d] = \int_{-0.5}^{0.5}\left(\hat{d_2} -
 d - \operatorname{E}[\hat {d_2} - d]\right)^2dd \approx 0.0013$$

In practice you just add the constant $\pm 0.07671320485$ to the estimate $\hat d$ of the same sign, giving you $\hat{d_2}$.

Your signal is low-pass and not critically sampled so a parabolic estimate of the peak of cross-correlation will be more accurate even without any such correction.

For a sinc peak shape, Eq. 1 can be inverted and used to remove the bias in the estimate.


Rate vs throughput vs capacity What is the difference between Rate , throughput and  capacity ? Their unit are all bps/s,but it seems that there are something different between them,can anyone tell me?

My research is optimizing something in energy harvest and SWIPT,and in the internet,it said the throughput and capacity have some difference,but i did't understand it.

WIKI:

Capacity: capacity is the tight upper bound on the rate at which information can be reliably transmitted over a communication channel.

throughput:
throughput is the rate of production or the rate at which something is processed.

rate formula : $time \times log_2(1+SINR)$

Because the unit of them are different,so i think maybe they should be the same,just apply to different place,but i am not sure about this!
When i maximize the rate,do i also maximize the throughput and  capacity ?

Rate: the actual rate you transmit data at.

Throughput: the percentage of signal data that is actually delivered correctly to the receiver in the overall data (signal data + coding + overhead ... etc).

Capacity: maximum rate that a channel allows for a given SNR such that the received signals can be reliably decoded at the receiver.


Visualize visual words in bag of words model I am implementing visual bag of words through these steps:

Find interest points using SIFT
Calculate SIFT descriptor
Build codebook through kmeans clustering of SIFT descriptors.

How can I visualize these visual words? For example in the following lecture, how do they get those patches?

Each SIFT descriptor corresponds to a region of the image. You take these from a bunch of images and group them into some number of clusters. I think what he's showing in the slide there is just a few samples from each cluster where he chose human-meaningful names for the clusters after the fact.


Source and channel encoding In information theory, there are two concepts that are widely used: source encoding and channel encoding.

I'm afraid I may not be capable of expressing my doubts clearly, but the thing is that I don't understand they interact between each other.

To try to make it clearer, I'll try to give an example of my doubts:

Every codeword for the channel encoding must have the same length. Why is this? If the source from which the channel is reciving its input used a Huffman code, for example, lengths might be different for each symbol. Isn't that a contradiction?

Every codeword for the channel encoding must have the same length.

That is not inherently true, but it's usually handy to consider the information that goes into the channel coder to be of constant length, especially when considering block coders.

  Why is this? If the source from which the channel is reciving its input used a Huffman code, for example, lengths might be different for each symbol. Isn't that a contradiction?

No. The point is that channel and source coding are two independent things. You take source data, and the source coder maximizes the average information per bit (by converting input words of constant length into source-coded words of varying length). Then the source-coded words are again simply considered as sequences of bits (or to be more general, sequences of fixed-length units of information) and channel-coded, optimizing robustness.


What is the Difference between Texture and Image? I am beginner in the image processing

I want to know what is difference between Texture and Image? or whether they are same? can anyone explain it.

Also,if i want to do texture segmentation,can i use same algorithms as applicable to the image segmentation?

An image (2D) is the projection of a natural scene or man made drawing on to a plane. The output of an optical sensor when displayed is an image.

Texture is a picture created by repeating a basic element called 'texel'. The texel can be a natural image or an artificial one.

The edge detection algorithms works on matrices. So if it is a normal image or texture image, you can do in the same way.

In texture images, recognizing the texel makes the rest of processing easy. Because you need to work on a smaller sized matrix.


Bandwidth expansion by computing missing harmonics in music? In old recordings or highly compressed audio files (lossy), the higher frequencies are lost, and we hear a distinctive muffling effect. (Just like if the song was played through a wall)

Is it possible to compute the missing harmonics using any method (fourier, bilinear, etc.) and add those back to the sound/song? What is currently feasible?

*It is true that those frequencies are lost, and I am not proposing to boost the higher frequencies (as it would raise the noise floor), but I thought there would be some way to approximate the missing frequencies, since we know what all the instruments should sound like (we know their corresponding harmonics)...*

I'm well aware of the fact "garbage in, garbage out", but even so, our brain itself have very advanced processing techniques. (3D visualization using only 2D sources in vision and audition) Wound't it be possible to fool the brain in "hearing" a better quality sound, and give old recordings back their original sharpness?

Edit: To be more clear, I'm proposing to make audio "sound" better, not to recover lost information.

Bandwidth expansion is actually a problem of missing audio data recovery, which is an active area of research.

Usually you assume that your signal follows some statistical model in the time-frequency domain, and the whole challenge is to infer the parameters of your model, based only on the available data.

In the case of bandwidth expansion, you would infer your model from the lower frequency part of your spectrogram, in order to recover the higher frequencies.

You can have a look at some of Paris Smaragdis' work on missing audio data and bandwidth expansion, using latent variable models. Some demos are available here.


Does delta-sigma ADC also reduce Gaussian noise on input signal to ADC or just quantization noise? The motive for posing this questions arises from a difference of analysis between a colleague an myself.

Our general environment is in the construction of an analog front-end which takes in signals to an instrumentation amp in the low microvolt range. followed by a gain stage which then feeds a 24 bit delta-sigma converter.

My understanding is that we can achieve a certain number of noise-free bits from the ADC, where the noise is the quantization noise. For any Gaussian noise on the input to ADC, this noise will not be reduced by the ADC itself. And in this regard, we need to do some additional averaging of the sampled data points for reducing the Gaussian noise as an additional processing step.

My colleague thinks that the delta-sigma ADC will also reduce the Gaussian noise on the input signal, and therefore additional averaging is not needed.

As a second question, is it the case that we should arrange the gain so that as much of the signal fills the dynamic range of the ADC so as to improve the SNR of the signal due to having a large signal quantized as a function of quantization noise. In other words, the quantization noise is a smaller fraction of the signal?

Could you help us out by commenting on this question?

No!  The ADC (delta sigma or not) can not reduce the uncertainty in the input.  It sounds to me that your friend has not made up a real signal flow diagram and then formed equations.  The answer to your second question is affirmative.  In fact it's generically  true that putting as much gain, without introducing extra, in the front end is a good idea.  But... think long and hard about input overloads distorting the signal and possibly overloading the amplifier circuit causing distortion and such; either in the front end amplifier or during subsequent signal processing.  In the "good"  old days overloading due to the input could cause numerous problems; including long term effects.  OTOH: detailed knowledge of the signal and noise spectral (or even time transient) characteristics can reveal ways to optimize your S/N .


Image Filtering by Averaging Similar Areas of the Same Image - Patch Based Processing I heard of filtering algorithms which, to filter an image specific area, find similar areas in the image and average them to reduce the noice of the original area. E.g. an edge of a house can be more or less the same in the whole image, so one can make many patches from that area and average them. It should be similar to using multiple shots of the same scene and averaging them to make a final image.

What's this algorithm called? I'd like to learn more about it as I found it pretty interesting. But couldn't find them on the Internet so far except some general averaging algorithms.

Those kind of algorithms are called Non Local algorithms.
The most known algorithms of this family is the - Non Local Means which is a decent Noise Reduction (Denoising) algorithm.

Until the Deep Learning boom, this approach has been extended and usually means working in the Patch Space of image - Patch Based Models and Algorithms for Image Denoising: A Comparative Review Between Patch Based Images Denoising Methods for Additive Noise Reduction.

All those are based on a model in patch based - Average by some kind of a metric / Similarity or some kind of Low Rank model


Power of a real exponential signal I am having difficulty with the case when a is real.

But after this, I am not going anywhere.

I think the correct expression for the power is
$$ P = \lim_{T\rightarrow\infty}\frac{e^{-aT}-e^{aT}}{2aT}. $$
When $T \rightarrow \infty$, $e^{aT}$ tends to infinity "faster" than $1/T$ and the limit does not exist.
However, an interesting thing happens when $a=0$. The signal $x(t) = e^{-at}$ becomes $e^0=1$, which is a power signal (with power equal to 1 watt).
I'm sure that (using L'Hopital or some other trick) it can be shown that $P=1$ when $a=0$, but I have not tried it. In any case, I would say that the premise of the question is wrong, and in fact $x(t)$ is a power signal for at least one (real) value of $a$.


Fourier Slice Theorem - Reconstruction Fourier Space I've stuck in one problem. I need to perform Fourier Slice Theorem on sinogram of medical image.

I read a lot about this theorem. I write a matlab code but results are always non-sense after inverse fourier transform. But Fourier Space seems alright.

I fill Fourier space with converting polar coordinates to cartesian coordinates. I apply 1 Dimensional Fast Fourier Transform with ftt command on every sinogram lines. After for each angle (I used 180 rotation angle) and radius value I convert it to cartesian coordinates and for each corresponding values G(theta, radius) ==> F(x,y) I filled up Fourier Space. It seems logical but results is not. How can I correct my code to run this algorithm?

Circular shape is Fourier Space, other nonsenseless image is ifft. Original image is binary liver but I can't add it.

Thank you!

Plus, my matlab code!

[![I=imread('binaryliver.png');

% I = im2bw(I, 0.1);
%Sinogram was calculated before for 180 angle!
\[w,h\] = size(I);

theta = 0:1:179;][1]][1]

xorg = floor(h/2);
yorg = floor(w/2);  %for find origin of matrix

F = zeros(w,h);     %fourier space assigment

for i = 1:length(theta)

    E(:,:,i) = fft(sinogram(i,:));  %calculate FFT for each line of sinogram

end

for i = 1:length(theta)

    for r = 1: length(E(:,:,1)) %Convert polar coordinates to cartesian coordinates

        x = xorg + (r-h/2+1)*cosd(-theta(i));
        y = yorg + (r-w/2+1)*sind(-theta(i));

        if x == 0 && y == 0   %

        else
            yy = round(y); xx = round(x);

            if yy <= 0

                yy = 1;

            elseif yy > h
                yy = h;

            end

            if xx <= 0

                xx = 1;

            elseif xx > w

                xx = w;

            end
        value = E(1,r,i);
        F(xx, yy) = value;
        end
    end

end

Im2=abs(ifft2(F));
figure; imshow(log(1+abs(F)),[]);

figure; imshow(Im2, []);

Fourier slice theorem states that Fourier Transform of your projections are equal to slices of 2D Fourier transform.

When you sample the projections over discrete angles the FT of projections become samples of 2D Fourier transform.

So you have to use your obtained samples to interpolate the remaining points.

in your code you put the samples in their position but the remaining points remain zero.


I need advice about how to make an audio frequency analyzer I just need some clarification of some very basic aspects of creating a frequency analyzer.  I am writing a c++ program, so far it is just a media player, with a custom explorer.  You can make playlists, shuffle, repeat etc.

But my goal is to continue to add functionality.  So far, I just have a signal level meter, and I am now trying to add a frequency analyzer.

I am using fftReal.  So I am at the stage where I have the output, and I am scaling to dbfs like so which I think is correct, but not exactly sure.

for (int i=0; i < fftsize/2; ++i)
    db[i]=20*log(sqrt(imaginary[i]*imaginary[i]+real[i]*real[i]));

So I know I should apply a window function somewhere, and I'm trying to understand how to do that.  But I am wondering this as well.  If I want to display the data as rectangles popping up and down rather than plotting points and connecting with lines, and I want to limit the amount of rectangles to about 50 or so.  Do I take the average db of a range.

I notice that in many applications, the frequency scale is skewed, where the higher frequency, the small distance between units.  What kind of scale is used, and how can I break up the spectrum so that I have a nice visual?

I'm sorry if this has all been discussed before, I have been googling for a while now, and overall I am overwhelmed with information, but there are a just a few key thing which I need to understand first.

Thanks

This really depends on how fancy you want to get. A "good" analyzer will typically do the following things.

Split the input into frames. The frames typically overlap and are windowed. Good choices are an overlap of 50% and a Hanning window.
Do an FFT
Select the center frequencies of bands. For audio good choices are octaves or third octaves (see http://www.engineeringtoolbox.com/octave-bands-frequency-limits-d_1602.html)
For each integrate the FFT energy around the center frequency. There are different choices for the integration windows. Good choices are trapezoids or butterworth windows. A really good window choice will make sure that all the energy is accounted for, i.e. that the sum of the band energies equals the sum of all FFT bins.

If the number of frequency bands is relatively small, it is more efficient and easier to simply use a bank of parallel bandpass filters and an RMS detector. This also allows to dial in the temporal behavior of the display nicely.


Other techniques like overlap save/overlap add I know overlap save and overlap add are used for long data sequence filtering. Are there any other similar or better techniques like these?

The traditional overlap-add technique cannot be used as is in musical applications (reverberation) where latency is critical. It is a common practice to break down the impulse response into non-uniform chunks, and to run in parallel different convolvers for each segment of the impulse response - some of them being naive FIR implementations (for the short blocks at the "head" of the response), and some of them being high-latency overlap-add FFT convolvers.

See this paper.


Filtering kinda "half-toned" image for OCR processing I have a scanned PDF material to which I want to add hidden text layer, so I could index the document. I used ghostscript black and white tiff output device (tiffg4) to extract pages as tiff images, and here is example of what they look like:

Processing this image with tesseract, does not give good results.
Changing ghostscript output DPI (600, 300, 150, 96) shows that image at 96 DPI gives best result from tesseract but it's still not satisfactory.

Now I thought to ask for advice which filter would enhance this image for OCR processing.

I could use imagemagick, or numpy/scipy/ndimage

What you really require is probably some morphological operation like dilation followed by erosion. This is called as Closing operation. May be in your case- just dilation itself might be good.

There was a similar question asked previously - which can help with other aspects.

monochrome (1-bit black-and-white) image conversion

How do I reconstruct text from an image using only morphological operations?


DFT and FFT application in signal processing Can anyone give me some hints to solve this problem

I prefer to give tips by comments but i don't have the right jet…
on the questions a, c and e you have look to the delta between 2 frequencies.
what is the space between each frequency and how much is the delta frequency of you're sampled time.
Probably questions B,D and f also but i'm not for shure because that a bit more complicated and had never done it myself...


Pre-calculation of Geortzel coefficient I have found an implementation of the Goertzel algorithm that contains code that I don't understand. And more specifically the pre-calculation of coefficients.

At the top of the code there is :

#define DTMF_697Hz  27980
#define DTMF_770Hz  26956
#define DTMF_852Hz  25701
#define DTMF_941Hz  24219

These defines are used inside the goertzel calculation like such :

for (i = 0; i < N; i++)
{
  v0 = ((cos_fact * v1) >> 14) - v2 + *x;
  x++;
  v2 = v1;
  v1 = v0;
}

x is my buffer.

v0, v1 and v2 are long int

cos_fact is the one of defines I showed above and passed as a function parameter.

My questions would be: What do these define correspond to? How to pre-calculate them for other frequencies? Is it normal that the values are so far apart?

ps: This code is part of ampm_open_lib and my goal would be to use it for detecting other frequencies than DTMFs tones.

EDIT:

Following the answer of Peter K. Here is a little C program that output the desired values.

#include        <stdio.h>
#include        <stdlib.h>
#include        <math.h>

#define         TAB     {697.0f, 770.0f, 852.0f, 941.0f, 1209.0f, 1336.0f, 1477.0f, 1633.0f}
#define         LEN     8

int             main(int ac, char **av)
{
  float         tab[LEN] = TAB;
  float         sr;

  if (ac >= 1)
    printf("Usage: ./a.out sample_rate\n");
  else
    sr = atof(av[1]);

  if (sr > 0)
    {
      printf("+-----------------------+-----------------------+\n");
      float coef = 0;
      for (unsigned int i = 0; i < LEN; i++)
        {
          coef = cos(2.0f * 3.14159f * (tab[i] / sr)) * 256 * 128;
          printf("|     %f      |       %f      |\n", tab[i], coef);
          printf("+-----------------------+-----------------------+\n");
        }
    }
  else
    {
      printf("You cannot use a sample rate that is smaller or equal to 0\n");
      exit (0);
    }

}

You can use it like this:

gcc yourFile.c -lm && ./a.out 8000

Thanks for the help.

Check out the comment in this code:

// ----- Parameters and Variables for Tone Detector -----
// cos = (cos (2*PI*(DTMF_Freq/8000.0))) * 256*128;
#define DTMF_697Hz   27980        // DTMF Row Frequency
#define DTMF_770Hz   26956
#define DTMF_852Hz   25701
#define DTMF_941Hz   24219
#define DTMF_1209Hz  19073        // DTMF Column Frequency
#define DTMF_1336Hz  16325
#define DTMF_1477Hz  13085
#define DTMF_1633Hz   9315
#define DTMF_1394Hz  15014        // DTMF Row Frequency    2nd harm
#define DTMF_1540Hz  11583
#define DTMF_1704Hz   7549
#define DTMF_1882Hz   3032
#define DTMF_2418Hz -10565        // DTMF Column Frequency 2nd harm
#define DTMF_2672Hz -16503
#define DTMF_2954Hz -22318
#define DTMF_3266Hz -27472

It explains how the cos_fact is being calculated pretty directly. For example, for the first item in the list:

$$
\cos( 2 \pi \times 697 / 8000) \times 256 \times 128 = 27979.5762684 \approx 27980
$$


Features for audio levels imbalance in films between speeh/music/effects What would be the key features of an audio track or set of tracks (like 5.1) for identifying film audio tracks that have a clear unbalance between speech levels and sound effects levels (apart from the obvious rms diference between tracks)? And how could they be extracted or calculated?

It's for a machine learning project about automatic identification of such files and ulterior processing/fixing. The imbalance could be on the same track, or in a combination of several ones.

Thanks in advance.

I would try a sum and difference calculation between the L and R channels to start. Usually when film audio is mixed, the speech is in both channels in phase while the special effects is out of phase.

So if you do:

X = L+R;
Y = L-R;

X should give you the signal with more speech and Y should give you the other effects. Then if you do an RMS calculation on X and Y, you can find how much are they off by.


What Is the Point of Doing the Zero Padding? What are the advantages and disadvantages of doing Zero-padding, in particular the case of speech signals?

You may think of it as efficient way to apply Dirichlet Window Based interpolation in the Fourier Domain.

The advantage of applying the interpolation using Zero Padding in the Time Domain is very simple - it is simpler and more computationally efficient.


Accurate measurement of object's position on a bayer-filtered and grayscaled images This is an experience-related question about cameras and image processing.

I have a laboratory setup. Under some simplifications it looks like a red ball rolling on a white table (the real setup is shown on the video). A computer vision system should measure X and Y positions of the ball.

For now I use a color camera (with Bayered output) and I find the ball by applying a simple color filter (InRange for R, G, and B components) to the frames. It works, but I would like to increase precision of the measurements.

Question:
Would accuracy of measurements be better in a case I used a grayscaled camera with the same parameters as bayer-camera?

I guess it should be so in a perfect case, i.e. when the ball is definitely highlighted on the frames. But in reality brightness filter usually works worse than color filter. On the other side I can change the shooting conditions, because this is laboratory setup. I can change the color of the objects, I can change lighting.

The key here is how well your acquisition system can differentiate the ball from everything else.  The other key is that the larger number of pixels on the object the better.

From a look at the video, and without seeing the grayscale version, I would think that a color camera (with the same pixel density as the greyscale camera) would do better.

One trick is:  greyscale cameras are often more sensitive to IR and, as a result, have an IR filter inserted just above the image sensor.  Your ball might be close enough to that so that, by removing the IR filter, the greyscale camera would do as well at seeing the ball as the color camera.

Another trick is:  Try illuminating the field of view of either camera with a different colored light.  You may find that the ball shows up better with one light source over the others (in one camera or the other).

Finally: try zooming in on the ball.  That will get you more pixels and either the color or greyscale camera will benefit from the extra pixels.


how to plot fft with imaginary and real part in 3d or calculate the degree of rotation I would like to make a plot like 1 and see the real and imaginary part in a 3d space. I dont want to make exactly the same plot. for me it is okay if i see the peak for both signals shifted. it is also ok for me if i can only calculate the rotation without plotting.
I tried it with a 3d plot python matplotlib
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
from mpl_toolkits.mplot3d import Axes3d

//create example Signal, sinus and cosinus-> they are shifted in time so i would like to see the rotation in complex plane

time = np.arange(1000)/1000)
signalSinus= np.sin(100*2.0*np.pi*time)
signalCosinus = np.sin(100*2.0*np.pi*time+np.pi)

fftSinus = np.fft.fft(signalSinus)
ftCosinus = np.fft.fft(signalCosinus)

fig=plt.figure()
ax=fig.gca(projection='3d')
ax.plot(np.fft.rfft(len(fftSinus),d=1/1000),np.imag(fftSinus),np.real(ffSinus)) // three axis plot but how to?
plt.show() //

i dont get right results it looks pretty strange to me :D if i plot both signals they look the same.

(1) nparange -> np.arange, and
(2) np.fft.rfft(len(fftSinus),d=1/1000)
You want the frequency values. You can start with just np.arange(1000) for normalized frequency, or frequency indices, and worry about actual values later.
Organize your ax.plot() as follows:

x = frequency index
y = imag
z = real


Fourier Transform: $\omega$ vs $f$ as frequency variable I try to understand how the Fourier transform changes when I try to compute $X(\omega)$ or $X(f)$. Can someone work me through the maths please for two examples, $x(t) = \exp(j\omega_0 t)$ and $x(t) = \cos(\omega_0t)$.
I know what their Fourier transform is in terms of $\omega$ but when I try to compute $X(f)$ from first principles, I am stuck.

You have the obvious relationship $\omega=2\pi f$. So going from a Fourier transform expressed in terms of $\omega$ to a Fourier transform expressed in terms of $f$ is just a matter of substitution.
The problem you've encountered is very likely in the argument of the Dirac delta impulse. Clearly, $\delta(\omega)=\delta(2\pi f)$, but you need to know the following relationship in order to obtain an expression in terms of $\delta(f)$:
$$\delta(ay)=\frac{1}{|a|}\delta(y),\qquad a\neq 0\tag{1}$$
So we have
$$\delta(\omega)=\delta(2\pi f)=\frac{1}{2\pi}\delta(f)\tag{2}$$
E.g., the Fourier transform of $x(t)=e^{j\omega_0t}=e^{j2\pi f_0t}$ is given by
$$\begin{align}\mathcal{F}\big\{e^{j\omega_0t}\big\}&=2\pi\delta(\omega-\omega_0)\\&=2\pi\delta\big(2\pi (f-f_0)\big)\\&=\delta(f-f_0)\tag{3}\end{align}$$
I'm sure that from this you can derive the corresponding expressions for $x(t)=\cos(\omega_0t)$.


Applying a Filter Kernel Defined by a Continuous Equation in the Frequency Domain I have a filter kernel $K(\omega)$ completely defined in the frequency domain as a continuous function of angular frequency $\omega$. I know $K(\omega)$ defined as a continuous mathematical equation.

Using the discrete FFT and IFFT, I would like to apply this kernel to a discrete time-domain signal using convolution via multiplication in the frequency domain.

What is a good way to deal with this convolution?  Should I generate the $U(\omega)$ as a continuous function in the frequency domain, take the IFFT, and then truncate, shift and window the signal as mentioned in the book Steven W. Smith - The Scientist and Engineer's Guide to Digital Signal Processing?

Alternately, should I take the continuous inverse Fourier transform of the function $K(\omega)$, and generate $k(t)$ completely in the time domain as a discrete signal by evaluating $k(t)$ for discrete values of $t$?  Then $k(t)$ can be treated as a FIR filter kernel and taken into the frequency domain.

Although this is a discrete 1D signal, are 2D image filtering kernels generated in the frequency or spatial domains?  What is the difference between the application of the kernels in the time and frequency domains?

Given an Image $ I \in \mathbb{R}^{m \times n} $ I would solve it as following (Using the Convolution Theorem):

I would sample the function in Fourier Domain into a grid of $ m \times n $.
I would apply DFT transform on the image to get the $ m \times n $ representation in the Fourier Domain.
I would apply Element Wise multiplication between the 2 $ m \times n $ arrays.
I would apply Inverse DFT on the result to get the filtered image in the spatial domain.

Pay attention that this will result in Cyclic Convolution and not Linear Convolution (Which requires more padding).


What is the maximum dynamic range of an N bit signal? In this presentation on slide 15 there is a corollary:

  Rule of thumb: the dark noise must be larger than 0.5

  Corollary: With a N bit digital signal you can deliver no more*) than
  N+1 bit dynamic range.

  *) You can if you use loss-less compression

They give an example:

  Example : A 102f camera with 11 bit dynamic range will deliver only 9
  bit in Mono8 mode. Use Mono16!

Why N+1?
How can an 8 bit signal deliver 9 bit dynamic range?

From taking a look at the presentation, it seems that they might be using a slightly different definition of dynamic range than is typical. Usually, as in the Wikipedia article on the topic, dynamic range is defined as:

$$
\text{dynamic range} = \frac{\text{largest possible representable value}}{\text{smallest possible representable value}}
$$

For an $N$-bit (unsigned) signal, this is equal to:

$$
\text{dynamic range} = \frac{2^N-1}{1} = 2^N-1
$$

However, their discussion of dynamic range is interspersed with discussion of quantization noise. Therefore, I posit that they instead define dynamic range as:

$$
\text{dynamic range} = \frac{\text{largest possible representable value}}{\text{largest possible quantization error}}
$$

For a uniformly-quantized quantity such as this, the maximum quantization error is equal to half of one bit. That leads to a dynamic range of:

$$
\text{dynamic range} = \frac{2^N-1}{0.5} = 2(2^N-1)
$$

The extra factor of 2 gives you an approximate increase of 1 bit in "dynamic range" when measured this way.  I assume that's what is meant by an $N$-bit digitized signal providing $N+1$ bits of dynamic range.


Effect of windowing on noise I understand that truncating a signal in time 'smears' the frequency response depending on the window chosen.  In general, the shorter the signal duration, the more 'flattened' the frequency response, as seen here (http://www.thefouriertransform.com/pairs/box.php):

However, how does the window length affect the frequency response of the (bandlimited additive white Gaussian) noise?  Assume a rectangular window of amplitude $A$, duration $T$, and a corresponding $\operatorname{sinc}(\cdot)$ main lobe in the frequency domain with amplitude $A\,T$ and width $\frac{2}{T}$:

$$\begin{align}
\mathscr{F}\bigg\{A \cdot \operatorname{rect}\left(\tfrac{t}{T}\right) \bigg\}
 &= \int_{-\infty}^{+\infty} A \cdot \operatorname{rect}\left(\tfrac{t}{T}\right) \, e^{-j2\pi ft} \ dt \\
\\
 &= \int_{\frac{-T}{2}}^{\frac{+T}{2}} A \, e^{-j2\pi ft} \ dt \\
\\
 &= A \, \frac{\sin(\pi fT)}{\pi f} \\
\\
 &= A \,T \, \operatorname{sinc}(fT) \\
\end{align}$$

If $A$ was fixed, and $T$ was halved, that would result in a $\operatorname{sinc}$ of halved amplitude but doubled main lobe width.  It would then seem that convolving this $\operatorname{sinc}$ would result in the 'same' amplitude of noise in the frequency domain because of the $\frac{1}{2}\cdot 2 = 1$ canceling.  That is, the effective noise bandwidth contributing to a given frequency is doubled, but the contribution per Hz of that bandwidth is halved.

Is that true? And in general, how does a window's duration and shape affect the frequency response of noise?
If (1) is true, does this imply that halving the window duration will also halve the SNR of a single sinusoid?  (Because the sinc magnitude of the signal is halved, but the noise floor remains constant)

Edit: One point I realized is that there may be destructive interference among noise components of different frequencies, and therefore this is not such a simple analysis as just convolving the fourier transform of the window function with the square root of the noise power spectral density.  Perhaps uniformly distributed noise phase at each frequency could be assumed?

I don't have access, but perhaps this paper is useful? http://ieeexplore.ieee.org/document/199437/

UPDATE: My previous response did not answer the OP's question. The following addresses the question directly:

Bottom line: Prior to windowing in time, each sample in frequency is an IID Gaussian random variable since the Fourier Transform of an AWGN waveform in time results in an identically distributed waveform in frequency (Gaussian distributed and to be white meaning each sample is independent of the next). After windowing in time, a dependence is created between the adjacent samples in frequency. But the overall frequency response will still be white (uniform and equal power overall) and Gaussian. The variance of a sine wave in relation to the variance/Hz of the white noise process (variance for an AWGN process must be given as a density in units of power/Hz as a truly white noise process has infinite power) will be unchanged in relation to each other; if the window caused the power of the sinewave to go down by one half, the power of the noise would also go down by on half. The actual values depend on how normalization is done in the computations, but for a straight power computation which is energy/time, reducing the window by one half (for example) would reduce the power by one half independent of what kind of waveform was involved (Sine, AWGN, etc).  This is in contrast to what would happen if we convolved with a rectangular window, which is covered in the second half of the post below (what was my original, but misguided, response).

Details:

For discrete time signals, consider the following from Parseval's Theorem which shows that the energy of the signal in time and frequency is the same:

When time goes from $-\infty$ to $+\infty$ which would be for the DTFT:

$$\sum_{n=-\infty}^{\infty}|x[n]|^2=\frac{1}{2\pi}\int_{-\pi}^{\pi}|X((e^{j\phi})|^2d\phi\tag{1}$$

Note when using normalized frequency (1) becomes the form below that is perhaps easier to follow:

$$\sum_{n=-\infty}^{\infty}|x[n]|^2=\int_{-0.5}^{0.5}|X(f)|^2df$$

When time is limited (windowed) would be for the DFT:

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

In the above DFT relationship using Parseval's Theorem we are comparing energy; if we further scale by M where M represents the total observation time in samples, we will then be comparing power under various rectangular window sizes of N samples which we can apply to both sinusoidal tones and white noise:

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

The DTFT case will not converge without any window applied (infinite energy) but we can get insight into the answer by considering an arbitrarily large window (the DFT) and then comparing that to what happens when we reduce it with a smaller window.

Sine Wave

Consider a sine wave with an arbitarily long window N with an observation time that also equals N:

If the window is indeed very large compared to a cycle of the sinewave, then the DFT of the sine wave will be well approximated by two impulses (as is the case exactly when the window is an integer number of cycles of the sinewave) each with a magnitude that is N/2 times the peak magnitude of the sine wave in time. Thus for a sine wave with an arbitrarily long window, Parseval's theorem results in the expected variance of a sine wave with peak $A_p$ (using M=N in Equation (3)):

$$\frac{1}{N^2}\sum_{k=0}^{N-1}|X[k]|^2 = \frac{1}{N^2}\left( \left(\frac{N}{2}A_p\right)^2+\left(\frac{N}{2}A_p\right)^2\right)=\frac{A_p^2}{2}=\sigma^2$$

As we reduce the window for the sinewave, the frequency response of the sine wave is indeed "smeared" to other bins; the impulses will become Sinc functions in frequency that will get wider as the window gets narrower, and the total power when considering the squared sum of all bins will go down as the ratio of N/M where M represents the original window size. Note that the total power of the original window size M will change in both domains if the residual fraction of a sine wave cycle becomes significant compared to the integrated area under one cycle squared, as is the case when the window duration is not significantly longer than one cycle of a sine wave. If we were considering a single complex exponential frequency tone, this variation as the window size became significantly reduced would not occur. However to be noted in either case, the power in time is equal to the power in frequency regardless of window duration and frequency of the tone (the power in both is equally effected).

AWGN

An additive Gaussian white noise process in time is an additive Gaussian white noise process in frequency, with the same distribution in both domains. (So therefore as far as a mathematical function it is just a change of variable from time to frequency when using a unitary Fourier transform). Let's also remind ourselves of what AWGN is conceptually: It is white, meaning it has equal power density over ALL frequencies (and therefore unlimited power and therefore not realizable), and Gaussian- meaning the distribution of its magnitude in time takes on a Gaussian shape. The Fourier transform of a Gaussian white process is also a Gausssian white process; what does that mean? In the frequency domain, the distribution in magnitude of the function versus frequency also takes on a Gaussian shape and in this case in terms of it being "white" it means explicitly that the transform of this function (the time domain function) has equal power over ALL time. Bottom line, as far as we are concerned, besides the variable defining the domain, the functions are identical.
With regard to Fourier transforms, multiplying by a window in one domain is convolution of the window kernel (Fourier Transform of the window) in the other domain. When we filter a signal, we convolve the signal with the impulse response of the filter, which is the inverse Fourier transform of the frequency response.  Further to be noted when working with the DFT as we have done above, the convolution itself is a circular convolution.

With that said, consider what would happen to the frequency response of an AWGN process when we window it in time: Prior to windowing, which is the case of an arbitrarily long window N with an observation time equal to N, the frequency response is indeed white, and as we noted above the "time response" is also similarly "white" in this case (meaning it extends over the full length with all the samples having a similar distribution). Also to note, relative to our sample time interval, each sample in time is uncorrelated from the next (therefore the resulting in a spectrum over our digital frequency interval that is indeed white). The variance of our time domain signal is equal to the variance of our DFT when we scale the DFT by N=M as shown in (3).

Just as in the case of the sine wave, if we reduce the rectangular window M to be less than M, the power (variance) will reduce by N/M, but what is interesting and pertinent to the question, is that the frequency response will remain white and Gaussian!  Why is this? By reducing the rectangular window to M, we are convolving the frequency response with a Sinc fucntion (or in our discrete system what well approximates a Sinc function for large M and is actually an "aliased" Sinc function), and as noted this is a circular convolution. Thus the frequency response would still be white, but to be noted we have created a dependence for each sample in frequency on adjacent samples due to the convolution operation. This means in frequency each sample is no longer independent from sample to sample, so in the time domain the transform will no longer be white- but in the frequency domain the amplitude distribution itself will still be Gaussian, and the power density will still be uniform over all frequencies within the digital frequency interval used and so therefore is is indeed still white in frequency.

Thus the impact of a rectangular window in time to the frequency domain is to remove the independence between the adjacent frequency samples, and reduce the overall power proportionally when compared over the same observation interval (equally as is done with a sine wave, so does not change SNR); but it does not change the statistical description of being white (in frequency) and Gaussian distributed.  The dependence between samples in frequency is similar to the effect of a dependence of samples in time: When we have a dependence between samples in time we have a band-limited (low pass filtered) process which we can therefore say is "frequency limited". When we have a dependence between samples in frequency we have a time-limited process; which is what the rectangular window is doing.

As a final point to help see what is going on; sometimes it easier to think in one domain instead of the other, so consider if we applied the rectangular window to any AWGN signal in frequency that is initially white (uniform density over all frequencies). Prior to windowing - the time domain signal would extend over our complete observation interval, and the DFT would extend over the complete frequency space defined by our sampling time interval.  When observing the signal in time, no matter how much we zoomed into the time domain waveform, it would appear as in the first plot below for AWGN, because every sample is independent of the next. And the historgram of the magnitude distribution is Gaussian. If we were to band-limit the frequency response (by multiplying the frequency response with a rectangular window), we would see in the time domain something similar to the second plot below; in that as we zoom in, we can see defined trajectories from one sample to the next! Note that the histogram of the magnitude (as long as we do it over enough samples) does not change and is still Gaussian. And important note that our time domain function still extends over our complete observation time with a uniform power- so it is "white" in time and Gaussian but it is no longer white in frequency. Thus we see directly what would happen to the frequency response in the case of the OP's question. Instead of the waveforms below being time, they would be frequency. The frequency response is still uniform in power (white) and Gaussian, but due to the windowing in time we would now be able to zoom in on the frequency response and observe the sample to sample correlation that would now exist that didn't exist prior to windowing. Prior to windowing each sample in frequency would be independent from adjacent samples so as we zoomed in on the frequency response it would continue to look like the first plot below. But if the time domain function was windowed, it would create dependence bewteen the adjacent samples in frequency and when we zoomed in to the frequency response in that case we would start to observe something like the second plot below: we would see a definite trajectory of the frequency response waveform as we move from one sample to the next- however it is still white (the power on average over all frequencies would be flat) and Gaussian distributed.

White Gaussian Noise (AWGN)

Band-Limited Gaussian Noise

A further way to prove that the frequency response remains white after multiplying the time domain function with a rectangular window is to observe the autocorrelation function in each case: The autocorrelation fucntion for an AWGN signal is an impulse, and the frequency response of an impulse is a uniform function.  Adding zeros to the AWGN fucntion (or equivalently windowing) does not change the result from being an impulse, and therefore the frequency response will still be uniform (white).  Adding zeros does interpolate between the existing samples in frequency, and thus the trajectories previously described are created... and to note from that, for a given window size of length T of an AWGN signal, the samples in frequency separated by 1/T will remain independent, but all samples in between will be dependent on the two adjacent samples separated by 1/T.

Previous post:
The following was initially given as a response but this is specific to convolving with a rectangular window which was not the question asked:

A windows duration and shape effects the spectral density of white noise based on the frequency response of the window directly. While noise will be reduced in power based on the relative length of the window; meaning as a sum of squares or $\int_0^T(x^2)dx$, while a sine wave within the correlation bandwidth of the window (meaning frequency < 1/T where T is the window length) would increase as a summation. I prefer to consider the window as a moving average such that the sine wave (if low enough in frequency) does not change and the noise is proportionally smaller. This just means we normalized the window to its length but is more intuitive that the window would not effect the sine-wave itself but would remove noise. The normalization if not used just results in an arbitrary scaling but the ratio of signal to noise is what is of interest in the end in either case.

Consider an example (digital) white noise process with total variance = 1

If we filtered this with a 10 tap unity gain filter (representing convolving the white noise process with a discrete rectangular window [1 1 1 1 1 1 1 1 1 1]), the noise from tap to tap in the filter would be uncorrelated, so would go up by the sqrt(10) in standard deviation (which represents its magnitude quantity), while a sine-wave that was within the filter bandwidth would be correlated and would increase by a factor of 10 in magnitude.

Observe the frequency response of such a filter, where the DC gain of 20dB represents the factor of 10 described above, as (20Log10(10)). This response shows exactly what would happen to the power level of a single tone at any frequency within the filters spectrum, while the power of multiple tones would be the sum of their individual powers (which is how we handle what happens to the noise, as in $\sum x^2 $ ) :

And the expected effect on the white noise

The noise is now shaped (colored) due to the lowpass nature of the window, and the overall noise after processing through this filter should only go up by 10log10(10) = 10 dB. Thus the SNR has increased 10 dB since the tone (signal) when up by 20 dB while the noise went up by 10dB, or if we normalize to the level of the tone, the noise has gone down by 10 dB or 1/10th in total power.

Testing this experimentally:

noise= randn(2^12,1);
var1 = std(noise);
noisefilt = filter(ones(10,1),1,noise);
var2 = std(noisefilt);
freqz(ones(10,1));  % frequency response

Results in var1 = 1.00355 and var2 = 10.64.

The increase is just a constant (and arbitrary) gain factor so what is important is how the noise is effected relative to a sine wave, in that the window reduces the noise power of white noise proportionally (in this case compare a wider window to one 1/10th in size and the smaller one removes 1/10th of the power) while reduces the sinewave according to a Sinc function with the first null at 1/T where T is the length of the window. (Or for any arbitrary window based on the Fourier transform of the window itself).

Also as I mentioned in the comment under the original posting, I believe fred harris handles the mathematics well in describing coherent vs non-coherent gain, equivalent noise bandwidth etc in windowed systems in this classic paper that I reference often: https://www.utdallas.edu/~cpb021000/EE%204361/Great%20DSP%20Papers/Harris%20on%20Windows.pdf


Amplitude not attenuated after bandpass filtering? I want to make a bandpass filter for an IQ signal recorded from an Software Defined Radio with center frequency "Fc" & sampling frequency "Fs". For futhur analysis I wanted to focus upon signal received between a fixed frequency range & filter out rest frequencies.

I generated "sinc" funtion using Python Scipy firwin then I convoluted it with my signal using Python fftconvolve. But as you see in the plots amplitude of the frequency bins outside the pass band of my are still not attenuated.

Approach

Read the data,convert it from uint8 & strore it as complex notation.
Generate a "sinc" function of width half of difference of cutoff frequencies eg. for 20MHz to 40MHz width of filter will be 10MHz. Since "sinc" is symmetric it will filter frequencies from -10MHz to +10MHz.
Move the centre of filter to mean of cutoff frequencies using frequency shift property.
Convolve the filter with the data.

Relevant Code

data_raw=signal.flatten()
data_raw=data_raw-127.5

data = np.empty(data_raw.shape[0]//2, dtype=np.complex128)
data.real = data_raw[::2]
data.imag = data_raw[1::2]

nyq = 0.5 * fs

# Normalized frequency to be used for frequency shift
centre=(highcut+lowcut)/2-fc

width=np.ceil((highcut-lowcut)/2)
h = firwin(513,width,nyq=nyq,window='hanning')
h = np.append(h,np.zeros(1024-len(h)))

t=np.arange(len(h))
h=h*(np.exp(1j*2*np.pi*t*(centre/fs)))
filtered_data = fftconvolve(data, h)

Details:

Fcentre = 137.65 MHz
Fsample = 2 MHz
Lowcut = 137.65 MHz
Hightcut = 138.15 MHz

Plots:

Firwin filter linear plot

Firwin filter semilog plot

Output after applying filter

You seem to be aware of the Fourier Transform pair that relates a sinc function in the time domain to a rectangle function in the frequency domain, but you don't seem to be able to correctly apply this to your situation. I am going to take one more stab at answering you. Let me know if you have questions.

A sinc filter can be designed using the relationship
\begin{equation}
\DeclareMathOperator{sin}{sin}
\DeclareMathOperator{sinc}{sinc}
\DeclareMathOperator{rect}{rect}
2 B \sinc\left(2 B t\right) \iff \rect\left(\frac{f}{2 B}\right)
\end{equation}
where $\sinc(x) = \frac{\sin(\pi x)}{\pi x}$. This shows us that the desired low pass bandwidth $B$ must be used in our expression to generate the sinc samples. Note that a modulated sinc has a bandwidth equal to $2B$ whereas a sinc centered at DC has a bandwidth of $B$.

In your example, you define several quantities. I am going to give these quantities symbols so I can do some math on them.

  Fcentre = 137.65 MHz

I'll call this $f_c$.

  Fsample = 2 MHz

I'll call this $f_s$.

  Lowcut = 137.65 MHz

I'll call this $f_l$.

  Hightcut = 138.15 MHz

I'll call this $f_h$.

For you to generate the sinc with the correct bandwidth, you would sample the sinc function above with
\begin{equation}
B = \frac{f_h - f_l}{2}.
\end{equation}

This means that your bandpass filter is given as
\begin{equation}
h[n] = \left(f_h - f_l\right) \sinc\left(\frac{f_h - f_l}{f_s} n\right) \exp\left(j \pi \frac{f_h + f_l - 2 f_c}{f_s} n \right)
\end{equation}
where I have ignored the required time shift to make the filter causal (implementable). The amount you shift it by will be determined by the length of the filter. The length of the filter will be determined by how much stopband attenuation you want.

Note that no finite length filter will yield zero gain in the stopband and all filters have a transition band (a band where the amplitude gradually goes from the desired passband value to the designed stopband value). You can't avoid this.

As a side note, you could have just as easily used the Remez Exchange algorithm (a.k.a. Parks-McClellan) to obtain an optimal filter. The window method is computationally simple, but it is in no sense optimal.


Can a causal filter without phase shifts exist? When I was studying dispersion of refraction index in semiconductors and dielectrics, my professor tried to explain that if a filter (like a dielectric absorbing some light frequencies, or an electric RC-filter) removes some frequencies, then the remaining ones must be phase shifted to compensate for those frequencies (which are infinitely spread in time as usual monochromatic signals) being subtracted from the whole signal, to preserve causality.

I intuitively understand what he was talking about, but what I'm not sure of is whether his argument is really justified - i.e. whether there can exist a non-trivial filter, which absorbs some frequencies and leaves remaining ones not shifted, but still preserving causality. I can't seem to construct one, but can't prove it doesn't exist as well.

So the question is: how can it be (dis)proved that a causal filter must shift phases of frequencies relative to each other?

Suppose that a linear filter has impulse response $h(t)$ and frequency response/transfer function $H(f) = \mathcal F [h(t)]$, where $H(f)$ has the property that $H(-f) = H^*(f)$ (conjugacy constraint).

Now, the response of this
filter to complex exponential input $x(t) = e^{j2\pi f t}$ is
$$y(t) = H(f)e^{j2\pi f t} = |H(f)|e^{j(2\pi f t + \angle H(f))}$$ and if
we want this filter to cause no phase shift, it must be that $\angle H(f) = 0$
for all $f$.

How about if, instead of no phase shift, we are willing to allow
a fixed constant phase shift for all frequencies? That is, $\angle H(f) = \theta$ for all $f$ is acceptable to us where $\theta$ need not be $0$? The extra
latitude does not help very much, because $\angle H(-f) = -\angle H(f)$, and
so $\angle H(f)$ cannot have fixed constant value for all $f$ unless that value
is $0$.

We conclude that if a filter does not change the phase at all, then
$H(f)$ is a real-valued function, and because of the conjugacy constraint, it is
also an even function of $f$.  But then its Fourier transform $h(t)$ is a
an even function of time, and thus the filter cannot be causal (except in trivial cases): if its impulse
response is nonzero for any particular $t > 0$, then it is also nonzero for $-t$ (where $-t < 0$).

Note that the filter need not be doing any
frequency suppression, that is, we did not need the assumption that
some frequencies are "removed" by the filter (as the OP's
professor's filter does) to prove the claim that zero phase shift
is not possible with a causal filter, frequency suppressor or not.


Proof of Forward Euler for discretizing a transfer function In Levine book "The control handbook" it is shown that, for discretizing
a transfer function $\frac{1}{s}$ using Forward Euler i simply have to replace s
with $\frac{z-1}{T}$. How can extend the proof and show that the substitution works for
every possible transfer function $G(s)$ ?

Another way to see how the forward Euler method approximates a continuous-time system is by considering the "ideal" mapping of the $s$-plane to the $z$-plane (why?):

$$z=e^{sT}\tag{1}$$

For frequencies that are much smaller than the sampling frequency (i.e., $|s|T\ll 1$) we can approximate $e^{sT}$ by its first order Taylor series:

$$z\approx 1+sT\tag{2}$$

which is exactly the forward Euler mapping, resulting in

$$s=\frac{z-1}{T}\tag{3}$$


Remove square-wave pulse artifact How can I remove a biphasic, square-wave-like artifacts from the signal, so I'll keep as much as possible information "under" the artifact? Below is the signal with artifacts:

The artifacts are those two big spike-like pulses. They are caused by square-wave pulses applied to the skin in neighborhood of recording electrodes.

Maybe you find this useful:

N=500;                          %signal length
t=0:N-1;
s=20*randn(size(t));
x=s;
p=200*[ones(1,5) -ones(1,5)];   %square artifact template
x(:,101:110)=x(:,101:110)+p;    %add artifact to signal at two random points
x(:,401:410)=x(:,401:410)+p;
figure;plot(x)
threshold = 1e5;                %peak detection threshold
c=xcorr(p,x);                   %cross-correlation
[pk l]=findpeaks(c,'MinPeakHeight',threshold);
l=N-l+1;                        %find the actual lag
y=x;
for j=1:numel(l)                %remove artifact
y(:,l(j):(l(j)+length(p)-1))=y(:,l(j):(l(j)+length(p)-1))-p;
end
hold on;
plot(y,'r--')
legend('corrupted signal','recovered signal')


Why is "$1/\mathrm{Hz}$" unchanged while converting$\mathrm{V}^2/\mathrm{Hz}$ to $\mathrm{dB}/\mathrm{Hz}$? When computing power spectral density via Fourier transform and Parseval's theorem, one gets units of $V^2/Hz$ (e.g. Wiki). As a popular next step, the following unit conversion is widely used: $V^2/Hz \to dB/Hz$, and it's done by taking the famous $10.\log10(\cdot)$.

Question: while $\mathrm{dB}$ is a log of a ratio of two $\mathrm{V}^2$'s (second $\mathrm{V}^2$ is a fictive one) and I get where $\mathrm{dB}$ comes from, I'm struggling to understand why $1/\mathrm{Hz}$ remains unchanged after taking the logarithm?

Alex

P.S. Some obscure explanation is in the last paragraph here (and I'm afraid the word 'loosely' conceals all the real stuff here):

  Since the argument of the logarithm is in units of $\mathrm{Hz}^{-1}$, this spectral measure can loosely be said to be in units of `dB/Hz'

Say you have a measurement in unit $U$ (V$^2$) per unit $u$ (Hz). I do claim (although not in public) that the logarithm is a kind of unitless transformation, see for instance What is the logarithm of a kilometer? Is it a dimensionless number? In other words, it is somehow $0$-homogeneous, while a power is $2$-homogeneous: multiplying by a constant yields a constant to some power.

First, the logarithm is  kind of homogeneous to a $0$-th power: it is the integral of $1/x$ (power $-1$). The derivative of an affine function is a constant. Derivatives decrease powers by one unit, integrals increase them of the same amount. The logarithm somehow stands with constants between $x^{-1}$ and $x^{1}$. Additionally, the coefficient of variation of $x^p$ when $p\to 0$ is about:

$$ \frac{x^p-x^0}{p-0}\approx \log x\,.$$

Second, a dB scale is often relative: to some signal measurement, to some "absolute" value. The $\log$ argument is better unitless: something in U divided by some unit $1$ U, called U$_0$. Thus, computing $\log(\mathrm{U}/\mathrm{U}_0)$ could be considered.

So  converting U/Hz into dB$_\text{U}$/Hz, to me, only affects the numerator.


Impulsive autocorrelation sequence I am looking for a proof that a finite length sequence can only have an impulsive autocorrelation sequence if the sequence is a single impulse itself, ie.:

  assume $s[n]$ is nonzero for $n=0,1,2,3$ only and compute
  $$r_s[k]=\sum_{n=0}^{3-|k|}s[n]s[n+|k|]$$
  for $k=0,1,2,3$.

What are all the solutions for
$s[n]$ that make $r_s[k]=1$ for $k=0$ and $r_s[k]=0$ otherwise?

Here is a detour over the frequency domain. We'll use index $k$ for frequency and index $n$ for time. Capital letters for Frequency Domain variables

The Fourier Transform (FT) of an impulse is $R(k)=1$
FT of the autocorrelation is the magnitude squared of the FT of the original signal, i.e $R(k) = |S(k)|^2$
Combining 1 and 2 we get $|S(k)| = 1 $,i.e. $S(k)$ must be an allpass filter
An allpass filter has zeros that are inverse of the poles
Any finite length sequence can be viewed as the impulse response of an FIR filter
An FIR filter has all the poles at $z=0$
An FIR allpass filter can therefore have only zeros at $z = \infty$
Zeros at infinity correspond to pure integer delays.
Pure delays have only one non-zero sample

And that's basically it. This says if $s[n]$ is a finite length sequence, the only way to get a unit impulse as auto correlation is to have only one non-zero sample. Any sample will work though, it doesn't have to be at $n=0$. For a length of 3 you get 4 solutions: $[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]$


Multiple cheap IMUs vs one expensive IMUs In my project, there is a possibility of either buying one expensive IMU or 5 cheaper IMUs. They have to be cheaper due to budget limitations, and the idea is to use them in order to be able to do "voting" - i.e. vote out outlier measurements of one IMU if 4 other IMUs show a different reading - and also for redundancy. Question: does anybody know if 5 cheaper IMUs (on the order to 300 dollars each) can create a signal at least as smooth, via sensor fusion, as one expensive IMU (on the order of 1500 dollars)?

That really will depend on the variance (standard deviation) of the measurements you get from the five vs the one.

Suppose your expensive IMU gives a measurement distributed as:
$$
x_\mbox{expensive}  \sim  N\left(x_\mbox{truth}, \sigma^2_\mbox{expensive}\right)
$$
and suppose your cheaper ones give measurements distributed as:
$$
x_\mbox{cheap}^i  \sim  N\left(x_\mbox{truth}, \sigma^2_\mbox{cheap}\right)
$$
for $i=1\ldots 5$ and that you form the average of these cheaper measurements to get the "real" measurement:
$$
\hat{x} = \frac{1}{5} \sum_{i=1}^5 x_\mbox{cheap}^i
$$
then $\hat{x}$ is distributed as:
$$
\hat{x} \sim  N\left(x_\mbox{truth}, \frac{\sigma^2_\mbox{cheap}}{5}\right)
$$
So if
$$
\frac{\sigma^2_\mbox{cheap}}{5} \le \sigma^2_\mbox{expensive}
$$
you should be OK with the cheaper sensors.

Of course, there may be other problems with your cheaper sensors: they may be biassed (so that the mean of their distributions will all be different from $x_\mbox{truth}$ ).

And I'm not quite sure what you mean by "voting"? Do you vote off the one furthest from the mean of the 5?  What if the mean changes for the remaining four?  I'll do a bit more digging to see if there's a viable algorithm there, but the direct mean of all five might be the easiest.