Offset after Butterworth bandpass filtering guess I'm stuck with this filtering problem and I'm thinking it has something to do with my raw data but I'm not so sure. So the problem is the following, please refer to the figures below, basically, all I'm trying to do is implement a bandpass Butterworth filter to filter the raw data trace (green). The raw data consists of 2501 points and was sampled at 20 kHz (100-ms). I'd like to bandpass filter this trace between 3 Hz and 170 Hz using the Matlab function 'butter' followed by either 'filter' (top figure) or 'filtfilt' (bottom figure), preferably the latter.
My Matlab code is as follows:
A = Average(1,:); % raw data
fs_Hz = 20000; % sampling rate (Hz)
order = 3;
fcutlow = 3; % Hz
fcuthigh = 170; % Hz
[b,a]  = butter(order,[fcutlow,fcuthigh]/(fs_Hz/2), 'bandpass');
x = filtfilt(b,a,A);
plot(x);
hold on
plot(A,'color','g');
legend('filtered data','raw data');

As you might guess the problem is the offset that I observe after both filtering functions. These 2501 point traces actually stem from a very long data-trace, should I possibly filter this long trace prior to extracting these smaller traces? Is the reason for the offset the fluctuating pattern of the trace? I can't detrend the data as I need the y-axis information for the subsequent analysis on this data. Should I be using a completely different type of filtering method; FIR instead of IIR??
Any advice on my problem would be greatly appreciated.
IljaPilja

Your data have a non-zero DC component which is filtered out by the bandpass filter. So the filter output has approximately zero DC, which causes the corresponding vertical shift of the signal. If you want to retain the DC information in the signal you must use a lowpass filter instead of a bandpass filter.


What does $I_{seg}$ in 3GPP TS 38.212 version 15.2.0 Release 15 represents? I'm working on 5G NR, whereby on reading Code block segmentation and code block CRC attachment given in Subclause 5.2.1 of 3GPP TS 38.212 version 15.2.0 Release 15.

What does $I_{seg}$ represent?

$I_{\mbox{seg}} = 1$.

if payload $A>1013$ or $A\ge 360$ and $G>1077$. Divided two block segment
PUCCH $n_{\mbox{max}}=10$;
PDCCH $n_{\mbox{max}}=9$;


changing circular shift in hologram reconstruction I'm trying to write a very simple code for the reconstruction of an in-line hologram of particles in a micro channel (digital holographic microscopy), but I'm getting some results that are hard for me to understand.

clear;

lambda = 632.8e-9/1.333; %wavelength
pix_size = 7.4e-6; %camera's pixel size
dx = pix_size;
dy = pix_size;

hologram = imread([file_path file_name]);
N = size(hologram);

tukey_win = tukeywin(N(1),0.075); %tukey window
% multiply image (single precision) with the tukey window to darken edges
trimmed_holo = single(hologram) .* (tukey_win * tukey_win');

prop_distance = 0.01:0.005:0.2; % propagation distances vector
fft_hologram = fft2(trimmed_holo); %Fourier transform of hologram
fft_shifted_hologram = fftshift(fft_hologram); %Fourier shift of hologram

for k = 1:length(prop_distance) % do multiple propagation distances
% obtain transfer function
G = trans_func(lambda, N(1), dx, dy, prop_distance(k));

%multiply FT of impulse response times FT of hologram
psi = G .* fft_shifted_hologram;

% get slice of reconstructed optical field through inverse FT
recon_image = ifftshift(ifft2(psi));

%write image to file
imwrite(uint8(abs(recon_image)), [file_path_save 'out' num2str(k,'%2.2d') '.tif']);

end

function G = trans_func(lambda, N, Delta_xi, Delta_eta, z)

M = N;
[m, n] = meshgrid(-N/2:N/2-1,-M/2:M/2-1);

first_term  = (lambda^2 * (n + (N^2 * Delta_xi^2)  ./ (2 * z * lambda)).^2) ...
    ./ (N^2 * Delta_xi^2);
second_term = (lambda^2 * (m + (M^2 * Delta_eta^2) ./ (2 * z * lambda)).^2) ...
    ./ (M^2 * Delta_eta^2);
arg_root = 1 - first_term - second_term;
arg_root(arg_root<0) = 0;

G = exp(-2*pi*1i*z/lambda * sqrt(arg_root));

The actual reconstruction appears to be correct, but I am getting a circular shift that varies with the reconstruction (propagation distance).
Here's the original hologram:

Here's z = 0.01, which (just) looks like I missed the ifftshift():

Here's at z = 0.0105:

At z = 0.011:

The circular shift keeps getting smaller, and by the last propagation distance (z = 0.2) it is practically gone:

I've looked at the code for quite some time and can't figure out what's causing the circular shift, so I would greatly appreciate it if anyone could explain to me what is causing/how to correct it. Thanks.

The problem is this line:

recon_image = ifftshift(ifft2(psi));

It should be

recon_image = ifft2(ifftshift(psi));

You need to shift the 0 frequency back to the corners before doing ifft2

I find fftshift / ifftshift and creating the frequency coordinates causes so many programming errors. This is my code for creating the coordinates (with 0 frequency in the corner):

function [ux,uy,r,th] = fftMesh(sizeV)
%FFTMESH - Create DFT cartesian and polar coordinates
%
% Inputs:
%
% sizeV     [width,height] of spectrum (same as image size)
%
% Outputs:
%
% ux        x in cartesian coords
% uy        y in cartesian coords
% r         radius in polar coords
% th        angle (rad CCW) in polar coords

rows = sizeV(1);
cols = sizeV(2);

[ux, uy] = meshgrid(...
    ([1:cols]-(fix(cols/2)+1))/(cols-mod(cols,2)), ...
    ([1:rows]-(fix(rows/2)+1))/(rows-mod(rows,2))...
    );
ux = ifftshift(ux);   % Quadrant shift to put 0 frequency at the corners
uy = ifftshift(uy);

% Convert to polar coordinates
th = atan2(uy,ux);
r = sqrt(ux.^2 + uy.^2);

end

If you use this, you dont have to worry about using fftshift / ifftshift. Alternatively, just remove the ifftshift lines in the function to keep the 0 freq in the middle.

Updated code:

clear;

lambda = 632.8e-9/1.333; %wavelength
pix_size = 7.4e-6; %camera's pixel size
dx = pix_size;
dy = pix_size;

hologram = double(rgb2gray(imread('test.jpg')));
N = size(hologram);

tukey_win = tukeywin(N(1),0.075); %tukey window
% multiply image (single precision) with the tukey window to darken edges
trimmed_holo = single(hologram) .* (tukey_win * tukey_win');

prop_distance = 0.01:0.005:0.2; % propagation distances vector
fft_hologram = fft2(trimmed_holo); %Fourier transform of hologram

for k = 1:length(prop_distance) % do multiple propagation distances
% obtain transfer function
G = trans_func(lambda, N(1), dx, dy, prop_distance(k));

%multiply FT of impulse response times FT of hologram
psi = G .* fft_hologram;

% get slice of reconstructed optical field through inverse FT
recon_image = ifft2(psi);

%write image to file
imagesc(real(recon_image)); pause;
imagesc(imag(recon_image)); pause;
imagesc(abs(recon_image))

end

function G = trans_func(lambda, N, Delta_xi, Delta_eta, z)

M = N;
[m, n] = meshgrid(...
    ([1:M]-(fix(M/2)+1))/(M-mod(M,2)), ...
    ([1:N]-(fix(N/2)+1))/(N-mod(N,2))...
    );
m = ifftshift(m);   % Quadrant shift to put 0 frequency at the corners
n = ifftshift(n);

first_term  = (lambda^2 * (n + (N^2 * Delta_xi^2)  ./ (2 * z * lambda)).^2) ...
    ./ (N^2 * Delta_xi^2);
second_term = (lambda^2 * (m + (M^2 * Delta_eta^2) ./ (2 * z * lambda)).^2) ...
    ./ (M^2 * Delta_eta^2);
arg_root = 1 - first_term - second_term;
arg_root(arg_root<0) = 0;


What are the most robust implementations of passive/self-localization assuming high noise and reverberation? I am doing research on a related subject that requires the use of TDOA estimation of two signals in an environment with a great deal of noise and reverberation. If anyone had any expertise on (or links to research regarding) robust passive localization/self-localization, it would be greatly appreciated - thanks!!

It is well known that the GCC-PHAT algorithm (based on signals correlation) is robust in the case of reverberation. You can also use preprocessing algorithms to reduce noise.


Which IIR filters approximate a Gaussian filter? So it recently dawned on me that Bessel filters, despite being listed along with the other common types, are really an oddball that belongs in a different "class", and I'm trying to learn more about it.

  The rectangular magnitude response represents the ideal frequency domain response, for the transition band is zero and the stopband has infinite attenuation. The Gaussian magnitude response, on the other hand, represents the ideal time-domain response, in that no overshoots occur in the impulse response and the step response. Many of the responses attained in practice are approximations to these ideal ones  source

So a brickwall filter is convolution with a sinc function, and has these frequency domain properties:

Flat passband
Zero stopband
Infinite roll-off rate/no transition band

It's non-causal and unrealizable because of the infinite tails in both directions.  It is approximated by these IIR filters, with the approximation improving as order increases:

Butterworth (maximally flat passband)
Chebyshev (maximum roll-off rate with stopband or passband ripple)
Elliptic (maximum roll-off rate with stopband and passband ripple)
Legendre  (maximum roll-off rate with monotonic passband)

The Gaussian filter is convolution with a Gaussian function, and has these time domain properties:

Zero overshoot
Minimal rise and fall time
Minimal group delay

It's unrealizable for the same reasons as the sinc function, and can be approximated by these IIR filters, more closely as order increases:

Bessel (maximally flat group delay)  according to 1 and 2

Here's Bessel filters of increasing order along with a Gaussian dashed line I picked merely because it seemed to fit the trend ($e^{-{1 \over 2}(\pi \omega)^2}$):

So my questions are:

Is everything right so far?  If so, are there other IIR filters that approximate the Gaussian?  What are they optimized for?  Maybe one that minimizes overshoot?

If you search for "IIR Gaussian" you can find a few things (Deriche? van Vliet?), but I don't know if they're really the same as a Bessel or if they optimize for some other property, etc.

The Deriche and van Vliet filters are heuristics.  In both cases they choose the locations of poles and zeros to minimize either the RMS difference or the maximum difference of the filter's impulse response from a Gaussian.

Both filters are causal-anti-causal pairs.  So I think they have no phase error or group delay, but you need to be able to run them backwards on the data as well as forwards.  This makes them popular in image processing, but perhaps restricts their applicability elsewhere.

That they are heuristics is evidenced by the wealth of papers that tweak them.  For example a google search (while I was looking for the link to the Deriche paper) turned up this one that tries to solve the problem that the Deriche derivative-of-gaussian filter doesn't have an exactly 0 DC response.  There are also some interesting issues about initializing the boundary conditions correctly.

I have found the following overview a good resource:
Dave Hale, Recursive Gaussian ﬁlters, Colorado School of Mines Center for Wave Phenomena report CWP-546.


How do I implement LIBSVM for multiclass classification using java? I am currently working on an image classification project ,where i need to use bag of features model with libsvm(HIK kernel) for classification.I have extracted sift feature and created a visual vocabulary,my fetaures are double values.But LIBSVM works with only integer values.Can anyone please help me with this?

libsvm doesn't work with only integer values -- in fact, you're supposed to scale your features to lie within [-1,1] or [0,1] before applying SVM's typically. Read their guide. They also have some tricks for dealing with categorical features as well.

If you want to get finitely many values which the feature can take, use quantization.


How to calculate to structure tensor? Structure Tensor is a matrix in form:

$S=\begin{pmatrix} W \ast I_x^2 & W \ast (I_xI_y)\\ W \ast (I_xI_y) & W \ast I_y^2 \end{pmatrix}$

where $W$ is a smoothing kernel (e.g a Gaussian kernel) and $I_x$ is gradient in the direction of $x$ and so on.

Therefore,size of structure tensor is $2N \times 2M$ (were $N$ is the image height and $M$ is its width). However it is supposed to be $2\times2$ matrix to decompose eigenvalues such that we obtain $\lambda_1$ and $\lambda_2$ as it is mentioned in many papers.

So, how to calculate $S$ matrix?

The $S$ matrix is finally of $2\times2$ the size of the initial image. Gathering it is a bigger image is somehow misleading. You effectively get a $2\times 2$ matrix at each pixel of the original image, by gathering values from the $4$ submatrices at corresponding locations.

Structure tensor provides local information. Do not forget to use smoothed directional derivatives as well when computing $I_x$ and $I_y$. You can for instance compute them with Difference of Gaussians (DoG) operators.
More discussion in  How do I implement the Structure tensor (in Matlab)? or Structure tensor tutorial.


Amplitude of an Image If I take a two dimensional image and conduct a fourier transform on it, I would get a two dimensional matrix of complex values. If I want to find the amplitude of each value, is that the same as calculating the magnitude of each component, ie Mag(re,im) = sqrt(re^2+im^2)?

I don't see the question here.

But you're right.
On the DFT map of the image, the magnitude / amplitude, etc. are powers of the length of the vector of the complex number.

On more useful visualization to help working with large Dynamic Range is using the log(1 + abs()) on the DFT.


Filter design / Analysis apps There are a plethora of tools, both commercial and free, which I have found online for designing filters. The ones I have tried (so far) prompt for frequency response, number of steps (for FIR), then generate coefficients and a frequency response plot.
Question: are there any particular sites or programs, which allow you to input or modify the coefficients? Free would be preferable, but (especially for an FIR with a large number of taps) commercial will do in a pinch.
(Matlab is too expensive for me, but I've heard of freeware Octave: if anyone recommends it for filter design I'd be willing to try it out)

I believe that Octave and Python/Scipy are the most commonly used free software packages for digital signal processing (including filter design). They offer of course much more than you might be looking for. On the other hand, if you want to play around with the coefficients and get some quick feedback, it might be very helpful to have one of these (or similar) packages because then you can easily modify the filter, plot frequency responses, and even write your own little programs if necessary. I've used Octave for that purpose and for me it worked fine. If you use Octave, you'll need to install the signal processing toolbox.

A totally different approach would be to choose your favorite programming language and write some simple routines yourself. This is perfectly feasible for some common FIR design methods (such as windowing). In combination with some plotting software this may be all you need. In this case I would consult this document on digital filtering. This latter approach is of course not the most efficient one because you'll be re-inventing the wheel, but in this way you'll learn how it all works and you will not just be using a black box design method.


Add noise/signal to mel log spectrogram I have a predicted mel log spectrogram $S$ and I want to add a noise signal $\eta$ to it.

S = log(mel-filterbank(abs/real part of(stft(signal)))

The naive way would be to approximate the linear spectrogram, add the noise signal and then apply the log mel transform again. However, it seems to me that this approximation adds unwanted artifacts which would not show up in a natural signal.

Is there a way to combine the signal without doing the approximation, given that I want the combined signal to be a mel log spectrogram too?

Except Logarithmic operation, other operations could be done on combined signal or separate signals both will give same result.

$STFT$($n$) $+$ $STFT$($S$) $=$ $STFT$($n$ $+$ $S$)
$mel$$(STFT$($n$) $+$ $STFT$($S$)) = $mel$($STFT$($n$ $+$ $S$)))

so just removing effect of logarithm and adding should do.
though I am not sure that artifacts will go or not, but you need not do reverse Mel(I am assuming that mel filter banks central frequencies and bandwidths are same for both $S$ and $n$).


Whats the difference between spatial and temporal resolution? Iam trying to understand how super-resolution works. But i think i have not understood correctly the difference between the optical resolution (spatial resolution?) and the resolution i know from a simple signal. I mean, i can have a higher sampling rate to improve resolution. Does this also effect images in the same way? And what about aliasing? Is aliasing in images difference than in a 1 D signal?
Here is a nice paper from google where they apply super-resolution on google pixel camera. They write:
Super-resolution techniques reconstruct a high-resolution signal from
multiple lower resolution representations.
They use multiple pictures taken at once, and use pixel shifts resulting from hand motion ... that pixel shifts and aliasing in the signal is used with kernel regression and other algorithms to reconstruct a higher resolution."The input must contain multiple aliased images, sampled
at different subpixel offsets. This will manifest as different
phases of false low frequencies in the input frames"
Is there a difference in applying this technique with Aliasing and multiple frames on 1d signals? How can i understand subpixel shifts applied on 1d signal.

A good 1d example of this is the foundation of the FFT algorithm in how an $N$ length DFT can be created from two $N/2$ length DFTs. If you look under the hood of this, we are increasing the resolution through multiple copies of a time domain signal each sampled at a different offset, and resulting in each signal containing the low frequency content as well as the aliasing of the high frequencies. The beauty is in the combining such that we can recover the low frequencies by adding the two FFT's and the high frequencies by subtracting the two FFT's (with an appropriate phase adjustment in frequency of one of the two before combining to compensate for the 1 sample shift in the time domain.

Let me demonstrate with formulas and graphics:

Given the general formula for an N length DFT:

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

Where $W_N^{nk}$ are the "roots of unity" phase rotations on a unit circle as $e^{-j2\pi nk/N} $

As further detailed in Cooley and Tukey's famous 1965 paper https://www.ams.org/journals/mcom/1965-19-090/S0025-5718-1965-0178586-1/ the equation above for an $N$ point DFT can be calculated from two $N/2$ point DFT's as:

$$ X[k] = \sum_{r=0}^{N/2-1}x[2r]W_{N/2}^{rk} + W_N^k\sum_{r=0}^{N/2-1}x[2r+1]W_{N/2}^{rk}$$

Here is what you can observe about each of teh two DFT's and how it relates to the OP's question:

Point 1: The first DFT is a DFT of all the even samples in x[n]

Point 2: The second DFT is a DFT of all the odd samples in x[n]

Point 3: The frequency response of $W_N^k$ has a magnitude of 1 for all frequencies and a phase that increases negatively from 0 to $2\pi$ as the frequency goes from 0 to $F_s$ where $F_s$ is the sampling rate (normalized radian frequency of $2\pi$ or normalized frequency in cycles/sample of 1).  This is exactly the same as the frequency response of $z^{-1}$, a unit delay of one sample (at the sample rate of x[n])!

First considering point 1 and 2 if we aligned the resulting even and odd samples in time. What we have done is decimated the sequence by two for the case of the even, and for the case of the odd we have (non-causally) advanced x[n] one sample and repeated the same decimate by two operation.

The transfer function of $z^{+1}$ is a magnitude of 1 for all frequencies but notably it advances the phase linearly from $0$ to $2\pi$ as we advance through all frequencies up to the sampling rate.

Next consider a digital spectrum so we can see how aliasing is handled in this case. The graphic below depicts a real spectrum where different symbols are used to differential the low frequency from the high frequency components in the first Nyquist zone extending from $0$ to $F_s/2$. The DFT of $N$ samples would return the block extending from $0$ to $F_s$ (with $mF_s$ actually cyclically repeating as bin 0 for all integers $m$).

Note what occurs to this spectrum when we compare the spectrum directly to the spectrum after it has gone through the $z^{+1}$ operation (we could equally say the bottom is the direct path and the top one goes through $z^{-1}$ since that could be actually implemented, but this is consistent with us aligning the output of the DFT result with the input without regard to processing delay and then will be consistent with the final formula, so at this point it is just math). The bottom plot doesn't quite show it since I couldn't draw a 3d spiral, but the phase shift is such that at the halfway point it will be $\pi$ or 180° representing a complete inversion of the spectrum such that it will be completely out of phase with the upper one, and then as it extends to the upper end it has spun around 360° so the spectrum at that point is back in perfect phase alignment with the upper one.

Next we see what happens when we decimate by 2 and how aliasing is created. When we sample any signal, all the spectrum around $mF_s$ for any integer $m$ is mapped via aliasing to $F=0$. So if we resample (which occurs when we select every other sample or decimate by 2), the same thing occurs, we have just created a new sampling rate. Most of the time when we decimate properly we are sure to low pass filter the signal first to eliminate anything in the middle of the spectrum that would alias (so true decimation is low pass filtering and down-sampling; here we are only down-sampling). In our case we want those images as we will be able to separate them with proper recombining:

Below shows the recombining where we recover the low frequency portion by summing the two and the high frequency portion by subtracting. The phase rotator $W_N^k$ is the summation undoes the opposite rotation that was innate in the even/odd FFT processing, and because of the difference in the rotation of the upper spectrum the aliasing can be isolated through adding and subtracting as depicted in the graphics. The subtraction occurs since each DFT here is only $N/2$ long and $W_N^k = -W_N^{k+N/2}$:

So note specific to your question how here was a 1D example of two lower resolution samples of the same data set with an offset in the sampling rate, and with it we are able to create the higher resolution data set (simply by interleaving the even and odd samples of course, but I believe this FFT view helps us see how the aliases are impacted and used to help create the higher frequency components).


How to convert ofdm symbol to signal? Hello may I know how to convert ofdm symbol generated from ifft block to an OFDM signal for the purpose of transmitting it?

I think if you added more details about your signal and the specifications of your hardware, that would be helpful to help you in that issue.
What you basically need is to multiply it by carrier signal and take the real part of the resulted signal, after that you need to do the inverse at the receiver end.  You might need also to upsample and use such filter equivalent to your system specifications.
Again, I think if you added more details about your system, we can also give you more information which may guide you.
Good luck


Why not use zero-signal in digital transmission? This is something that bothers me from time to time and has always been hand-waved away with some plausible but non-rigorous explanation. For clarity, consider the QPSK case where we have the constellation points $(0,1), (1,0),(-1,0),(0,-1)$. What is the fundamental flaw in using the origin, and then three equidistant points separated by $\frac{2\pi}{3}$ radians? We would have the same minimum distance but with smaller average energy.

I understand the added annoyance in synchronization and differentiating between "no message" and "the zero message". But these should be solvable with preambles and whatnot. I'm curious about the flaws that are more than just annoyances, the ones that makes this constellation a fundamentally bad idea.

EDIT: I spoke without thinking, clearly the minimum distance would not be the same, it would be one in the zero-signal case and $\sqrt{2}$ in the standard case. The number of pairs at min. distance would be... uh... six...?

All implementation aspects aside, the constellation you propose performs worse than QPSK in an additve white gaussian noise (AWGN) channel.

I claim this based on simulations that I have run with Matlab calculating the symbol error rate (SER) as a function of signal-to-noise ratio (SNR). Here is the result:

As you can see, for a given SNR, the proposed constellation has an higher SER than a QPSK constellation. Additionally, I expect the proposed scheme to perform even worse in terms of bit error rate (BER) as -- in contrast to QPSK -- Gray labelling is not possible. Some details concerning the simulation:

1 million bits per SNR value
maximum likelihood decision
additive, Gaussian distributed noise
memoryless channel (no inter-symbol interference)

It might well be that I made a mistake writing the simulation script, so I uploaded it in case you would like to check my implementation.

This might not be the mathematically rigorous explanation that you have been hoping for. To answer your question theoretically, one would have to calculate the SER analytically. For that, one should first find the decision thresholds, which is the main problem (for me, at least). It's far from obvious, what the optimal decision thresholds would be. They are not straight lines as with QAM constellations but probably a kind of zero-centered disc with three "rays" emerging from its center. Calculation of the SER should probably involve a transformation of the problem to polar coordinates.

With that kind of irregular constellation it also doesn't help much to compare the minimum distance, i.e. the length of a straight line between two given constellation points. It's true, that for a given mean signal energy, the outer points of the proposed contellation have greater distance to the inner point than any two points of the QPSK constellation. However, the inner point of the proposed contellation is surrounded by three other points which  decreases its decision area and consequently increases the symbol error probability for that point. It's not obvious what is better. According to the simulation, QPSK is better.

Update

The areas of the maximum likelihood decision can also be simulated. In the following image, each color corresponds to a distinct, decided symbol. The input to the decision device is a very noisy (0 dB SNR) signal with the proposed constellation.

We can clearly see how the inner point is "squeezed" between the outer points. In my opinion, that is the reason why the proposed constellation performs worse than QPSK. Now that I see this plot, I find the SER calculation not so difficult anymore. With a little time it could be done.


Detecting multi-path propagation when one acoustic signal is dominant I am sending an acoustic signal from a source to a receiver where there are limited number of paths the signal can travel. I want to detect if there has been multipath propagation or not when one path is the dominant line of sight path. Based on this link - Cross-correlation peaks in acoustic multi-path conditions , it seems like cross-correlation and checking for multiple peaks is an option, but it does not work when there is Fading (one dominant signal and others being attenuated).

It seems like this sort of Fading is called Rician Fading - however, I am not sure how to use this to detect if their has been multipath propagation or not. How to use the source and receiver signals to know if multipath propagation happened or not. (Also, I cannot use time of arrival of chirp pulses since the difference between time of arrivals of different paths is less than the chirp width)

It seems to me that, in the conditions you describe, the most you will see is flat fading. One indication that flat fading is occuring is that a change in the positions of transmitter and receiver, and/or the reflectors causing the multipath, produce a significant change in the received power. Usually the change has to be of the order of $\lambda/2$ meters, so this may be unfeasible in your case.


quick conversion of dB to $\log_2$ I was listening to this lecture

https://ocw.mit.edu/courses/electrical-engineering-and-computer-science/6-451-principles-of-digital-communication-ii-spring-2005/video-lectures/lecture-1-introduction-sampling-theorem

and the prof computed $\log_2(37dB)$ in no time. I tried to follow his logic but couldn't get it. Can someone please explain how can one do such a computation in head. He used $\log_{10}(2)=0.3$ as a building block for his computation


Suggest an image segmentation technique to find the paste radii? I am new to MATLAB/Digital Image Processing. So pardon me for any typing errors or wrong use of jargon.

1) I am trying to find the paste diameters for a total of 100 images in sequence. I am posting only one of the RGB images here.

2) I have attempted using imfindcircles and regionprops to find the diameters. imfindcircles demands a min and max radius range and since I have evolving paste diameters (for 100 images) it does not seem to be applicable. After some image segmentation and applying regionprops, the diameters were also not accurate.

3) What I have achieved so far using the image segmentation techniques is:

using the following code:

A = imread('ref.jpg');
I = imsharpen(A);
I = rgb2gray(I);
I = imadjust(I);
I = im2bw(I);
I = imcomplement(I);
I = im2bw(I, .8);
I = bwareaopen(I,2000);
nhood = true(9);
I = imclose(I,nhood);
I = imfill(I,'holes');
I = imcomplement(I);
I = bwareafilt(I,1,'largest');
imshow(I)

4) But I actually require is this (A perfect circle). With this I can find the pixel area and then the diameters:

5) Please suggest an alternate way to measure the diameter of the paste or help me reach step 4. I need to segment it from the mould (noise) above. Also I would like to reiterate that I have a bundle of images so any user-defined technique such as imdistline wont work.

Please let me know if you have any follow-up questions

Thankyou

This problem can probably be solved easier if there is a known relationship between the camera position and the metallic device that is over the paste. If the device and the camera are both fixed, for example, it would be easy to use a picture of the device without the paste to eliminate it from the other pictures.

But here is one suggestion for measuring the diameter of the paste using only this given image. I'm using DIPimage, because I know it best. I'm sure the same can be accomplished with the Image Processing Toolbox.

% Read in image:
img = readim('https://i.stack.imgur.com/5FLCs.jpg');
img = img{2}; % Use the green channel
% Alternating sequential filters to keep only the dark big circular shapes
paste = opening(img,15,'elliptic','add max');
paste = closing(paste,30,'elliptic','add max');
paste = opening(paste,60,'elliptic','add max');
paste = closing(paste,90,'elliptic','add max');
paste = opening(paste,120,'elliptic','add max');

The sequence of opening and closing are meant to remove specular reflections and much of the metallic device. Here is what the image looks like after this step:

There are some bulges to the left still, I didn't put any effort into trying to remove those, because there is sufficient outline for correct measurement.

% Find paste outline:
labs = label(~threshold(paste,'otsu')); % threshold and label
D = traceobjects(labs,[],2,'convex hull'); % get object outlines
[~,index] = max(cellfun(@length,D)); % find the longest outline
D = D{index}; % keep only the longest outline
% Display:
dipshow(img)
hold on
plot(D(:,1),D(:,2),'r')

This shows the convex hull of the outline over the original image, to show that the filtering didn't move the edges of the paste:

Next, I'll use code to list antipodal pairs that I get from this blog post. These are points on opposite side of the convex hull. This is the reason I used traceobjects with the 'convex hull' option, the algorithm to list antipodal pairs only works with a convex hull.

% Signed area of triangle
area = @(pt1,pt2,pt3) ( (pt2(1)-pt1(1))*(pt3(2)-pt1(2)) - ...
                        (pt2(2)-pt1(2))*(pt3(1)-pt1(1)) ) / 2;
% Height of triangle, pt3 is the top vertex
height = @(pt1,pt2,pt3) ( (pt2(1)-pt1(1))*(pt3(2)-pt1(2)) - ...
                          (pt2(2)-pt1(2))*(pt3(1)-pt1(1)) ) / ...
                                               norm(pt2-pt1);
% Next point along the polygon
next = @(p,N) mod(p,N)+1;

diameters = [];
N = size(D,1);
p0 = N;
p = 1;
q = 2;
while area(D(p,:),D(next(p,N),:),D(next(q,N),:)) > ...
      area(D(p,:),D(next(p,N),:),D(q,:))
   q = next(q,N);
end
q0 = q;
while p~=p0
   p = next(p,N);
   listq = q;                       % p,q is an antipodal pair
   while area(D(p,:),D(next(p,N),:),D(next(q,N),:)) > ...
         area(D(p,:),D(next(p,N),:),D(q,:))
      q = next(q,N);
      listq(end+1) = q;             % p,q is an antipodal pair
   end
   if area(D(p,:),D(next(p,N),:),D(next(q,N),:)) == ...
      area(D(p,:),D(next(p,N),:),D(q,:))
      listq(end+1) = next(q,N);     % p,q+1 is an antipodal pair
   end
   disp(listq)
   for ii=1:length(listq)
      v = D(p,:)-D(q,:);
      d = norm(v);
      diameters(end+1) = d;
   end
end

Finally, we can use these distances between antipodal pairs to find the average diameter of the paste. We'll assume that about 1/3 of the boundary is affected by the bulges on the left that I didn't completely remove. Therefore we'll use the 2/3 smallest values, and take the mean of those:

diameters = sort(diameters);
diameters = diameters(1:round(length(diameters)*2/3));
disp(mean(diameters))

This gives me 773.25 pixels. You will need to know the size of a pixel to translate this value to physical units.


MFCC for very short signals I have very short signals (usually 100 samples from a 44100 hz audio recording ). These are 'bang' noises that a hammer makes on a piece of metal. I want to use these short signals as input to some analysis (and a machine learning algorithm). I am not experienced in sound processing but from what I've gathered one of the most informative features you can use from a sound recording is the MFCC. I am using python's scikit talkbox package to compute the mfcc of my signals. The problem is that this package requires the signal to be at least 300 samples long (and as i've said mine are about 100 samples long and sometimes shorter).

I wanted to know if there is a different process I should be doing in this case to calculate the MFCC or if there are other suggestions for other methods, besides the MFCC, that I should be using.

Have you tried zero-padding your signals? By adding zeroes to the front and/or end of the signals you could make them longer than 300 samples without affecting the outcome of your analysis.

Bear in mind that with very short signals your frequency-domain resolution is going to be pretty limited, and an impulsive sound such as that you describe is likely to have a very broad spectrum. The impact this has will depend on what you plan to do with your MFCCs, but essentially, there'll be fewer of them!

Using an existing package to perform the signal processing is generally a great idea, as long as you're happy to trust it! I'd certainly encourage you to make sure you understand all the processing being carried out on your data - this will affect your ability to demonstrate that your work is trustworthy.


is there a second order approximation to perceived loudness Quoted from this

  As a first order approximation it takes ten times the energy to double the perceived loudness

So is there a 2nd order approximation. If so then how much change in RMS value doubles perceived loudness? Is it governed by any system of equations?

Thanks

This depends highly on the signal and it's content. For narrow band signal, the loudness can fairly well be estimated through the equal loudness curves as published in ISO 226 (see for example) http://en.wikipedia.org/wiki/Equal-loudness_contour

For wide band signal, things are more complicated. If the signals are stationary, you apply frequency weighting curves for the appropriate overall loudness (such as "A" or "C" weighting), which gives in some cases a reasonable estimate but can also be quite off by  a fair amount.

Next step up in complexity would be a loudness calculation based on on a detailed perceptual model. So you split the signal into "critical bands", apply forward, backward, simultaneously masking thresholds, apply frequency weightings and than integrate this all up. The resulting loudness is often measure in "Sone" which intended to be a "ratio scale", meaning it has an absolute zero value (threshold of audibility) and 6 Sones is twice as loud as 3 Sone.

So overall it's possible to get very accurate estimates of perceived loudness even for non stationary signals, but it requires deploying very complicated perceptual models. These models basically the same ones as used in perceptual codecs such as MP3, AAC, Dolby Digital etc.


Performing a Deconvolution So I have the following data:

[[-344    13771   4600 ]
 [-275.2  12478   6410 ]
 [-206.4  19443    830 ]
 [-137.6  69392   3830 ]
 [ -68.8  143737  3780 ]
 [   0    189278 16870 ]
 [  68.8  184486  5090 ]
 [ 137.6  188466  9380 ]
 [ 206.4  185023 21680 ]
 [ 275.2  128133  1460 ]
 [ 344    51288   1950 ]
 [ 412.8  10854   4290 ]]

First column is the x value (position in microns).  Second value is the recorded data point (y), and the third column is the error (noise).

I am wondering how I can find the response of my device if the data I record is known to be a convolution of the response and a slit of width 83.6666 (plus noise).

My attempt was to use the convolution theorem of fourier transforms and use MatLab's fft() to solve for the desired function, but I could not figure out how to get everything the same length.  I also thought about being blunt and using deconv() but that gave me something very strange, (with my rect function defined as

>> x=-344:1:412;
>> h=zeros(size(x));
>> h(300:386)=1;

So I used deconv(h,y) and got a very strange looking plot that was not a function.

So I was hoping someone could give me some advice on how to perform this mathematically (code or no code).  Especially with the addition of noise.

Thanks!

en.wikipedia.org/wiki/Wiener_filter

what you are looking for is a Weiner  filter  that takes into account the standard deviations of your data. See http://www.astroml.org/book_figures/chapter10/fig_wiener_filter.html
 and the description in Wikipedia  for implementation.  Basically you will do what you were suggesting,  but with an additional term in the denominator of the ratio of fft's  that is the fft  of the noise.


Delay filters in Python using for loop and $\tt lfilter$ I have these two filters as given in recursive form:
\begin{align}
y_1[i] &= y[i] - y[i-d] + 1.2y_1[i-d] - 0.8y_1[i-2d]\\
y_2[i] &= y[i] - 2y[i-d] + y[i-2d]
\end{align}

where $d$ is an integer for setting the delay. I wanted to have a look and see how such filters change input signals depending on the value of $d$, but I faced the problem with implementing the actual filters in Python: I tried iterations and scipy.signal.lfilter, and each method was giving me different results. Let me show you the code:

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import lfilter
# The test signal is just a sine of the length l, constant amplitude of 1
# and varying frequency w
l = 1024
x = np.array(range(l))
w = x/l
y = np.sin(w*x)
d = 4

# Let's caclulate the signals filtered with IIR filter
## First using forward iteration
y11 = y.copy()
for i in range(2*d, l):
    y11[i] = y[i] - y[i-d] + 1.2*y11[i-d] - 0.8*y11[i-2*d]
## Second using backward iteration
y12 = y.copy()
for i in reversed(range(2*d, l)):
    y12[i] = y[i] - y[i-d] + 1.2*y12[i-d] - 0.8*y12[i-2*d]
## Third using lfilter
y13 = lfilter([1, 0, 0, 0, 1], [1, 0, 0, 0, -1.2, 0, 0, 0, 0.8], y)

# Let's caclulate the signals filtered with FIR filter
## First using forward iteration
y21 = y.copy()
for i in range(2*d, l):
    y21[i] = y[i] - 2*y[i-d] + y[i-2*d]
## Second using backward iteration
y22 = y.copy()
for i in reversed(range(2*d, l)):
    y22[i] = y[i] - 2*y[i-d] + y[i-2*d]
## Third using lfilter
y23 = lfilter([1, 0, 0, 0, 2, 0, 0, 0, -1], [1], y)

# Let's start the figure and plot the signal in the top row, then the filtered
# signals below it
fig = plt.figure()
ax1 = plt.subplot(421)
plt.title('Test signal')
plt.plot(x, y)
ax2 = plt.subplot(422)
plt.title('Test signal')
plt.plot(x, y)

# The left column is for the IIR-filtered signals
plt.subplot(423, sharex=ax1)
plt.title('Forward iter.')
plt.plot(x, y11)
plt.subplot(425, sharex=ax1)
plt.title('Backward iter.')
plt.plot(x, y12)
plt.subplot(427, sharex=ax1)
plt.title('lfilter')
plt.plot(x, y13)

# The right column is for the FIR-filtered signals
plt.subplot(424, sharex=ax2)
plt.title('Forward iter.')
plt.plot(x, y21)
plt.subplot(426, sharex=ax2)
plt.title('Backward iter.')
plt.plot(x, y22)
plt.subplot(428, sharex=ax2)
plt.title('lfilter')
plt.plot(x, y23)

plt.tight_layout()
plt.show()

And here's the resulting figure:

Judging by the look of the plots in rows 2 to 4 one can see that each calculation produced different results except the 2nd and 3rd plot in the right column where backward and forward iteration worked out similarly. Normally I would rely on the results from Scipy's lfilter, but there I'm not sure if I specified the coefficients correctly. Somebody please take look and point out the mistakes.

After a trial-and-error session I figured what was wrong in my code:

The coefficients I fed to lfilter had incorrect signs, because only
the signs of the coeffs from the a array have to be reversed and
those in the array b stays the same.
The backward iteration did not work properly for the IIR filter,
because when running backwards the delayed values were not yet
filtered and in that case equal to the input signal.

Here is the fixed code:

l = 1024
x = np.array(range(l))
w = x/l
y = np.sin(w*x)
d = 4

# Let's caclulate the signals filtered with IIR filter
## First using forward iteration
y11 = y.copy()
for i in range(2*d, l):
    y11[i] = y[i] - y[i-d] + 1.2*y11[i-d] - 0.8*y11[i-2*d]
## Second using backward iteration
y12 = y.copy()
for i in reversed(range(2*d, l)):
    y12[i] = y[i] - y[i-d] + 1.2*y12[i-d] - 0.8*y12[i-2*d]
## Third using lfilter
# was:
#y13 = lfilter([1, 0, 0, 0, 1], [1, 0, 0, 0, -1.2, 0, 0, 0, 0.8], y)
y13 = lfilter([1, 0, 0, 0, -1], [1, 0, 0, 0, -1.2, 0, 0, 0, 0.8], y)

# Let's caclulate the signals filtered with FIR filter
## First using forward iteration
y21 = y.copy()
for i in range(2*d, l):
    y21[i] = y[i] - 2*y[i-d] + y[i-2*d]
## Second using backward iteration
y22 = y.copy()
for i in reversed(range(2*d, l)):
    y22[i] = y[i] - 2*y[i-d] + y[i-2*d]
## Third using lfilter
# was:
#y23 = lfilter([1, 0, 0, 0, 2, 0, 0, 0, -1], [1], y)
y23 = lfilter([1, 0, 0, 0, -2, 0, 0, 0, 1], [1], y)

# Let's start the figure and plot the signal in the top row, then the filtered
# signals below it
fig = plt.figure()
ax1 = plt.subplot(421)
plt.title('Test signal')
plt.plot(x, y)
ax2 = plt.subplot(422)
plt.title('Test signal')
plt.plot(x, y)

# The left column is for the IIR-filtered signals
plt.subplot(423, sharex=ax1)
plt.title('Forward iter.')
plt.plot(x, y11)
plt.subplot(425, sharex=ax1)
plt.title('Backward iter.')
plt.plot(x, y12)
plt.subplot(427, sharex=ax1)
plt.title('lfilter')
plt.plot(x, y13)

# The right column is for the FIR-filtered signals
plt.subplot(424, sharex=ax2)
plt.title('Forward iter.')
plt.plot(x, y21)
plt.subplot(426, sharex=ax2)
plt.title('Backward iter.')
plt.plot(x, y22)
plt.subplot(428, sharex=ax2)
plt.title('lfilter')
plt.plot(x, y23)

plt.tight_layout()
plt.show()

Now on the plots it is clearly visible that iterations produce the same results as scipy's lfilter except for the backward iteration attempt at the IIR filter as it was explained above.


Accelerometer sensor's data analysis Here is a sample data that I got from the accelerometer sensor:

The first selected area shows the acceleration of the IMU device to right side and the second shows the movement to left side.
I am looking for an algorithm to detect selected areas (as an array of floats) and provide the direction.
I almost don't know anything in this area (signal processing), So I tried to avoid using incorrect keywords, feel free to edit this post.

First, learn about digital differentiators and how to implement them with software. Then apply your accelerometer output signal to a digital differentiator. Next, use software to detect a large positive-valued sample followed by a large negative-valued sample in the differentiator's output sequence. That sample will correspond to your "first selected area". Finally, use software to detect a large negative-valued sample followed by a large positive-valued sample in the differentiator's output sequence. That sample will correspond to your "second selected area".


Array Signal Processing: What are snapshots in Direction of Arrival estimation? I do have an overview of basic radar signal processing chain [Data cube, 3 stages of FFT, Range, Range-Doppler map, coherent/non-coherent integration, coarse DoA etc.].

I am now trying to understand an advanced DoA estimation technique-MUSIC and one of the road blocks is the snapshot concept. Nearly every paper retrieved from google assumes an prior understanding of snapshots and am yet to find one that tells me what exactly is it.

So, what are snapshots in DoA estimation?

A snapshot is simply a data capture (sample window) of ADC samples from all receivers, that we know holds all the information to resolve what we need (distance and/or direction) given the various dimensions involved (transmit pulse length, nearest possible target distance, farthest possible distance, maximum array dimension, and the echo wavelet length). This defines the delay (from start of transmit time) of the window and the length of the window. The literature refers to a "snapshot" because many systems can repeatedly take snapshots to integrate-out noise, etc., but other cheaper systems can't.


What statistic is used to determine presence of a signal in noise? This is a detector problem I believe:

I am being stumped by what appears to be a simple problem. Basically, I have a band of interest. If signal energies exist within this band of interest, then I perform operation X on my signal.

My problem is that I am not sure exactly how to go about 'deciding' if a signal exists or not. In that, after I perform an FFT, I can look for peaks.

But now what?

Is the statistic used usually comparing this peak to the surrounding mean of the spectrum? Or is it some other statistic?
What statistical measure do I use to simply determine if a signal is present, and go from there?
How do I set this value? Simple thresholding?

EDIT Based on feedback:

For this simple case, I am assuming a tone, in white gaussian noise. What I am trying to get  a handle on are:

How exactly does one generated a ROC curve. Does one have to go and label all the data first, and then get the true-positive and false-positive rates for a multitude of thresholds?
How does decreasing SNR affect the ROC curve? Move it towards the diagonal?
What is adaptive thesholding doing to a given ROC curve that was otherwise generated without an adaptive threshold?

3a. What are some common adaptive threshold techniques I can look at that are common?

This is one of the oldest signal processing problems, and a simple form is likely to be encountered in an introduction to detection theory. There are theoretical and practical approaches to solving such a problem, which may or may not overlap depending upon the specific application.

A first step toward understanding the approaches to the problem is understanding how you would measure the performance of your signal presence detector. There are two important and related metrics used to quantitatively measure how good a detector is: its probability of detection $P_d$ and its probability of false alarm $P_{fa}$.

$P_d$ is specified as the probability that your detector will indicate the presence of the signal of interest, given that the signal is actually there. Conversely, $P_{fa}$ is the probability that your detector will indicate the presence of the signal of interest, given that the signal is not there. As you might expect, then, in a perfect world, we would design a system that yields $P_d = 1$ and $P_{fa} = 0$ and call it a day. As you might also expect, it's not that easy. There is an inherent tradeoff between the two metrics; typically if you do something that improves one, you will observe some degradation in the other.

A simple example: if you are looking for the presence of a pulse against a background of noise, you might decide to set a threshold somewhere above the "typical" noise level and decide to indicate presence of the signal of interest if your detection statistic breaks above threshold. Want a really low false-alarm probability? Set the threshold high. But then, the probability of detection might decrease significantly if the elevated threshold is at or above the expected signal power level!

To visualize the $P_d$ / $P_{fa}$ relationship, the two quantities are often plotted against one another on a receiver operating characteristic curve. Here's an example from Wikipedia:

An ideal detector would have a ROC curve that hugs the top of the plot; that is, it could provide guaranteed detection for any false alarm rate. In reality, a detector will have a characteristic that looks like those plotted above; increasing the probability of detection will also increase the false alarm rate, and vice versa.

From a theoretical perspective, therefore, these types of problems boil down to selecting some balance between detection performance and false-alarm probability. How that balance is described mathematically depends upon your statistical model for the random process that the detector observes. The model will typically have two states, or hypotheses:

$$
H_0: \text{no signal is present}
$$
$$
H_1: \text{signal is present}
$$

Typically, the statistic that the detector observes would have one of two distributions, according to which hypothesis is true. The detector then applies some sort of test that is used to determine the true hypothesis and therefore whether the signal is present or not. The distributions of the detection statistic is a function of the signal model that you choose as appropriate for your application.

Common signal models are the detection of a pulse-amplitude-modulated signal against a background of additive white Gaussian noise (AWGN). While that description is somewhat specific to digital communications, many problems can be mapped to that or a similar model. Specifically, if you are looking for a constant-valued tone localized in time against a background of AWGN, and the detector observes the signal magnitude, that statistic will have a Rayleigh distribution if no tone is present and a Rician distribution if one is present.

Once a statistical model has been developed, the detector's decision rule must be specified. This can be as complicated as you desire, based on what makes sense for your application. Ideally, you would want to make a decision that is optimal in some sense, based on your knowledge of the distribution of the detection statistic under both hypotheses, the probability of each hypothesis being true, and the relative cost of being wrong about either hypothesis (which I'll talk more about in a bit). Bayesian decision theory can be used as a framework for approaching this aspect of the problem from a theoretical perspective.

In the simplest practical case, the detector might trigger a detection if the detection statistic exceeds a fixed threshold $T$. In a more complicated and practical case, the detector might have some criteria for setting an adaptive threshold $T(t)$ and trigger a detection at time $t$ if the detection statistic breaks the threshold value at that instant. In your description of your problem, you hit on one common method for setting such an adaptive threshold: calculate the neighborhood mean to estimate the "background level", then set a detection threshold some amount above that mean. This can work for some applications, and there are many other ways to arrive at such a threshold.

Given a statistical model for the detector's input and the decision rule used to map that statistic to detection conclusions, one can then calculate the detector's theoretical performance metrics. In the design phase, you would typically calculate these metrics as functions of the free design parameters that you have (for instance, the threshold $T$ above). You can then evaluate the inherent tradeoffs: "if I set $T=5$, then I get $P_d = 0.9999$, but $P_{fa} = 0.01$. That's too high of a false alarm rate, so I better increase the threshold."

Where you eventually decide to sit on the performance curve is up to you, and is an important design parameter. The right performance point to choose depends upon the relative cost of the two types of possible failures: is it worse for your detector to miss an occurrence of the signal when it happens or to register an occurrence of the signal when it hasn't happened? An example: a fictitious ballistic-missile-detector-with-automatic-strikeback-capability would be best served to have a very false alarm rate; starting a world war because of a spurious detection would be unfortunate. An example of the converse situation would be a communication receiver used for safety-of-life applications; if you want to have maximum confidence that it doesn't fail to receive any distress messages, it should have a very high probability of detection of the transmitted signal.


First order low pass filter I am trying to better understand the first-order low pass filter:

Summary:

Per wikipedia,  a first order low pass filter yields the following in discrete time:
$$
\frac{Y(s)}{U(s)}= \frac{\omega_{c}}{s+\omega_{c}}
$$
yields
$$
y[k] = \left(\frac{\omega_c T_{s}}{1+\omega_{c} T_{s}}
\right) u[k]+\left(\frac{1}{1+\omega_c T_{s}}\right) y[k-1]
$$
or
$$
y[k] =\alpha \, u[k] \ + \ (1 - \alpha)y[k-1]
$$

where

$
\begin{array}{cclc}
\omega_{c}
&:&
\text{Cutoff angular frequency of filter}
&
[\frac{rad}{s}]
\\
T_{s}
&:&
\text{Sampling period}
&
[s]
\\
\end{array}
$

Question 1:

Even though this filter is in discrete time, does it still model an analog ($s$-plane) filter?

If I wanted to use a discrete computing system to filter in real-time,
would I need to use the digital ($z$-plane) equivalent?
If so, what is the general process for performing this?
My best guess is:

Determine digital cutoff frequency $\omega_d$.
Convert to analog cutoff frequency $\omega_a =
      \frac{2}{T} \cdot \mathrm{tan}(\omega_d \cdot \frac{T}{2})$.
Determine transfer function for (first-order low-pass) analog filter
using analog cutoff frequency $\omega_a$.
Transform to transfer function for digital filter
using bilinear transform $z = \frac{2}{T} \cdot \frac{z-1}{z+1}$

Relation to exponential smoothing:

On the same page, exponential smoothing is referenced.
The exponential smoothing page describes an exponential weighted average as:
$$
y[k] =\alpha u[k]+\left(1 - \alpha\right)y[k-1]\quad\text{where}\quad\alpha = 1 - e^{-\omega_{c} \cdot T_{s}}$$

Question 2:

How is it possible to relate the first-order low-pass filter alpha
with the exponential smoothing alpha?

To answer your first question, yes, you will need to convert the signal into Z-plane. Bi-linear transformation is one way to achieve the desired result. You can even try Impulse Invariant Transform to convert the analog filter to a digital one.


Evaluation of Autocorrelation and Power Spectral Density of white noise through a filter So say there's a filter with an impulse response of $h(t) = (0.8)^t u(t)$. I'd like to pass white noise through this and figure out the autocorrelation and power spectral density of the output.

I'm having a bit of trouble figuring out exactly how to do this. I know if you pass a signal through a filter, you simply take the Fourier transforms of the signal and the impulse response, multiply them together, and that's the output.

However, not being any mathematical expression for the noise, I'm unsure how to proceed. Any tips would be appreciated.

White noise is modelled as a process having spectral density $G_n(f) = k$, a constant.

After the filter, the resulting process $y$ will have spectral density $G_y(f) = G_n(f) |H(f)|^2$. And the autocorrelation will be the inverse Fourier transform of the spectral density.

This is a general result for any stable filter and finite power, wide-sense stationary input process (in this case we are pushing the boundaries because white noise has infinite power as defined above).


FIR filter design with flat passband but equi-ripple stop-band I am looking for a lowpass FIR filter with flat passband but equi-ripple stopband. In other words, it likes Chebyshev_filter Type II except that it is FIR instead of IIR. Linear-phase is preferred.

Thanks

I suggest you the Parks-McClellan method. In Matlab you can use $\tt{firpm}$.

Matlab's (now obsolete function) $\tt{remez}$ also uses this scheme.

The FIR filter is optimally designed to approximate e.g. a Chebyshev filter such that the maximum error between the filter's response and the desired response is minimized.

For more information see this and the examples therein.


Why is interpolation a time varying system I was reading about interpolation (Interpolation and Decimation of Digital Signals - A tutorial Review, Ronald E. Crochiere) and found that Interpolation filter is a time varying system. Can someone give an simple example why is it a time varying system. Because, as far as I understand, interpolation just stuffs zeroes between the samples to increase the sample rate by appropriate amount, even if we apply the input at a later time it should give the same output.
Also, even though the process is time varying we can represent it by using transfer functions. For instance, Zero order hold is represented by sinc filter. Why is this possible?

Considering a discrete time framework, the interpolation operation has two stages; first expand the signal by zero stuffing between the existing samples, and then lowpass filter the expanded signal to get the interpolated samples. The second operation is LTI, hence it's the first stage which is not LTI.

The expansion block is shown like $$x[n]\longrightarrow \boxed{ \uparrow L } \longrightarrow y[n]$$ and mathematically defined as

$$ y[n] = \begin{cases}{ x[\frac{n}{L}] = x[m] ~~~ , ~~~ \text{ for } n = m L~~ ,~~  m=0,\pm 1, \pm 2,... \\  ~~~~~ ~~~~ ~~~ 0  ~~~~~ ~~~~ ~ , ~~~ \text{ otherwise } }\end{cases} $$

Now apply the test for time-invariance  on $y[n]$ :

First shift the ipnut. Let $x_d[n] = x[n-d]$, then call $y[n]$ as $y_2[n]$ as:

$$ y_2[n] = \begin{cases}{ x_d[\frac{n}{L}] = x_d[m] ~~~ , ~~~ \text{ for } n = m L~~ ,~~  m=0,\pm 1, \pm 2,... \\  ~~~~~ ~~~~ ~~~ 0  ~~~ ~~~~~~~~~~~~ ~, ~~~ \text{ otherwise }}\end{cases} $$

which becomes:
$$ y_2[n] = \begin{cases}{ x[\frac{n}{L} - d]  = x[m-d] ~~~ , ~~~ \text{ for } n = m L~~ ,~~  m=0,\pm 1, \pm 2,... \\  ~~~~~ ~~~~ ~~~ 0   ~~~ ~~~~~~~~~~~~~~~~~ ~ , ~~~ \text{ otherwise } }\end{cases} $$

Second, similarly shift the output $y[n]$ by $d$ which is denoted as $y_d[n] = y[n-d]$:

$$ y_d[n] = y[n-d] = \begin{cases}{ x[\frac{n-d}{L}]=x[m] ~~~ , ~~~ \text{ for } n-d = m L~~ ,~~  m=0,\pm 1, \pm 2,... \\  ~~~~~ ~ ~~~~ ~~~~ 0   ~~~ ~~~~~~~~~~ ~, ~~~ \text{ otherwise }  }\end{cases} $$

which becomes:
$$ y_d[n] = \begin{cases}{ x[\frac{n-d}{L}] = x[m] ~~~ , ~~~ \text{ for } n = m L + d~~ ,~~  m=0,\pm 1, \pm 2,... \\  ~~~~~ ~~~~ ~~~ 0  ~~~ ~~~~~~~~~~ ~, ~~~ \text{ otherwise } }\end{cases} $$

As you can see, $y_2[n]$ and $y_d[n]$ are not the same signals. You can explicitly see this by evaluating them for specific values of $n$ and $d$, such as $d=1$ and $n = L$, then $y_d[L] = 0$ while $y_2[L] = x[0]$ and hence $y_d[L] \neq y_2[L]$. And therefore the conclusion that the expansion system (or as a consequence the interpolation system) is a Time-Varying and not LTI system.


KLT vs DCT basis restriction error I want to implement 8x8 DCT and 8x8 KLT matrices to the 8-length vectors of 256x256 Lena image (e.g. each row of image there is 32 vectors 256/8 = 32 and for all rows 256x32 = 8192). I do the eigen analysis for finding KLT matrix in MATLAB as:

co = corr(I); % I is Lena image
co = co(1:8,1:8); % 8x8 autocorrelatin matrix
[V,D] = eig(co); % eigen analysis
eigval = diag(D); % eigen values
[desc, ind] = sort(eigval,'descend'); % sorting eigen values for optimum KLT
eigvec = V(:,ind); %sorted eigenvectors = 8x8 KLT matrix
y = transpose(eigvec*transpose(I(t,m:(m+7)))); % KLT matrix*8-lenght vectors-- t=1:256 and m=1:8:256 iterations of for loop
xtilda = transpose(inv(eigvec)*transpose(y)); % inverse KLT matrix*basis restricted vectors
MSE(k) = MSE(k) + immse(I(t,m:(m+7)),xtilda); % MSE between restricted vector and original vector, k iteration of for loop
MSEfor32vec(g) = mean(MSE); % average error for one row=32 vectors, g iteration of for loop

In the same way I calculated MSEfor32vec of 8x8 DCT matrix using dctmtx command.

When I got the result, I was expecting KLT had less error; but DCT had less error. In the most of signal processing books, it is concluded that KLT must have less error.
What is my mistake or if my solution is true what is the point that I miss in the books?
I tried to explain my question in the best way, please ask me if it is not crystal clear..

Natural images can seem complicated. However, the set of humanly interpretable images is relatively limited, with respect to all possible images ($24^{x\textrm{millions of pixels}}$).
There is a belief that their structure can be better summarized. Indeed, lossless compaction can achieve 2-3 factor ratios, and lossy compression about 8-16  with little perceivable loss.

Most of the present image compression methods (that I  know of) are based on some kind of interpolation/extrapolation schemes. The rationale is to predict some samples, and to code the hopefully smaller prediction errors either on less bits (lossless), or with less precision (lossy). Those schemes can be applied either in the original space domain, or in a transformed space: either on fixed bases (Fourier, discrete cosine transforms, wavelets, etc.), on adaptive bases (PCA, NMF, KLT), learnt dictionaries, machine learning or artificial intelligence dabases.

Adaptivity sounds sound. One can thus adapt to image features. Globally on images? Why not, but the computational prize can be important: eigen-algorithms can be $O(N^3)$. One can reduce the size to smaller image patches. You have to compute eigenvectors for each patch, passed on its second-ordre statistics (covariance).
Then, less memory is required, and this can be parallelized easily. Adaptive bases are thought optimal in several cases to pack information. But they come with a price for compression. Since bases are different from patch to patch, one should code the eigenvectors and transmit them as well. For traditional rate/distortion rates, this is not favorable.

So is there a good-enough fixed-basis? Among the many choices, it turned out than one discrete cosine transform, the DCT-II, could do a good job, either theoretically or practically. This is celebrated in a 1974 paper by N. Ahmed et al.: Discrete cosine transform
It is better suited to small patches, and possesses fast algorithms. There are other options, adaptations, but let's keep it simple for now.

So, assuming that one is working on patches, and computing covariance, DCT should be close to KLT in energy, and somehow better as one does not have to code vectors.

Below are previous answers related to properties of the DCT

How to understand the de-correlation property of DCT? What does de-correlation mean?
Could a DCT be used for an audio magnitude spectrum rather than DFT?
Why do we use scaling factors in DCT? Are they important?
Number of cosine waves in DCT
How to get the DCT energy distribution image in MATLAB,or we may call it the DCT coefficients distribution image?


Gabor filters for digital image processing For the application of Gabor filters in feature-based (e.g. lines and edges in images) digital image processing, it is desired that the DC component is 0. Why is this so? What effect does it have mathematically if the DC component of the filter is not 0?  Thanks in advance!

In short: If a filter have DC component it means that for a constant region of the image, the filter response will depend on the grey value of that region. Constant region with higher values ==> higher response.

One usually wants higher responses only in case of edges/corners and not in the case of bright but constant regions.

As a picture worth 1000 words, here is an example. I took a simple image consisting of two horisontal lines, each one in different background. I took the Gabor and LogGabor filter. LogGabor is a variant that have 0 DC response.

I filtered the image with a Gabor filter (tuned to respond to horizontal signals) and then, the same image, with LogGabor filter (tuned with exactly the same parameters as Gabor filter).

Original image:

Gabor convoluted:

LogGabor convoluted:

Another example of Gabor filter, with a bandwidth of 0.7 octaves.

Now, let's ignore the edges for a moment and focus on the constant regions.
In the second image, you see different responses for the constant regions! Darker region have a lower response than the lighter one. For the LogGabor case (third image) you see that the response is the same regardless of the original intensity.

Note: The response values were scaled in 0-255 interval, so, neutral grey means 0 response, darker means negative values and lighter positive ones. The bandwidth is set to 2 octaves (except for the 4'th image)

EDIT

According to a paper by Boukerroui et al., for bandwidths around 0.7 the DC response is close to the computational errors ($10^{-6}$). The 4'th figure demonstrate this.

My personal advice to you is to use another bandpass filter. I strongly recommend Gauss Derivatives or Difference of Gaussians (in this order) if it is necessary for you to stay in spatial domain. If you don't have this restriction, I personally prefer LogGabor.

From my experience there is no net advantage or disadvantage in terms of end result between logGabor, Gaussian Derivatives, Difference of Gaussians or more exotic Cauchy filters. BUT this might vary, depending on your practical application.

As for the code, here is a snippet of my Gabor/LogGabor implementation.
The code is in Matlab, the convolutions are done in Fourier space. The complexity is $O(n log n)$ and it depends only on the image dimensions, not on the filter dimensions. The code follows the guidelines in Boukerroui et al.. It is developed by me and released under MIT licence (not shown here) so feel free to use it as you see fit.

Img must be a MxNx1 of doubles, (preferably between 0 and 1 but it should work with 0-255 too); P0 is the tuning period (wavelength), that is, $P_0=1/f$, where f is tuning frequency; Orient is the filter orientation, expressed in radians; FBW frequency bandwidth in octaves; ABW the angle bandwidth.

function [ ReConv,ImConv ] = GaborConv( Img, P0,Orient, FBW, ABW,varargin )
 ReConv=Img; ImConv=Img;
 ImgFFt=fft2(Img);
 lg2pi=sqrt(log(2)/pi);
 A=P0*lg2pi*(2^FBW+1)/(2^FBW-1); B=P0*lg2pi*(1/tan(0.5*ABW));
 K=1/(A*B); KF=2*pi/P0;
 [N,M]=size(Img); u0=1/P0*cos(Orient); v0=1/P0*sin(Orient);
 FFilt=zeros([M,N]);
 stx=1/(M-1); sty=1/(N-1);
 [U V]=meshgrid(-0.5:stx:0.5,-0.5:sty:0.5);
 Ur=(U-u0)*cos(Orient)+(V-v0)*sin(Orient);
 Vr=-(U-u0)*sin(Orient)+(V-v0)*cos(Orient);
 FFilt=exp(-pi*(Ur.*Ur.*A.*A+Vr.*Vr.*B.*B)) ;
 FFilt=ifftshift(FFilt);
 normFact=sum(sum((abs(ifft2(FFilt)))));
 RezFreq=FFilt.*ImgFFt/normFact; RezSpace=ifft2(RezFreq);
 ReConv=real(RezSpace); ImConv=imag(RezSpace);
end

Log Gabor method:

function [  ReConv,ImConv ] = LogGaborConv( Img, P0,Orient, FBW, ABW,useCos ).
 ReConv=Img;
 ImConv=Img;
 ImgFFt=fft2(Img);
 %Compute the k_beta according to the desired frequency bandwidth
 kBeta=exp(-1/4*sqrt(2*log(2))*FBW);
 gaussSD=sqrt(ABW^2/(2*log(2))); %compute Gauss SD to ensure that Gauss(ABW,0,SD)==0.5
 lnkB2=2*log(kBeta)^2;
 nc=exp(-1/8*lnkB2/2)*sqrt(-2*sqrt(pi)/(1/P0*log(kBeta)));
 [N M]=size(Img);
 stx=1/(M-1);
 sty=1/(N-1);
 [Ug Vg]=meshgrid(-0.5:stx:0.5,-0.5:sty:0.5);
 Radius=sqrt(Ug.^2+Vg.^2);
 Radius = ifftshift(Radius);
 Radius(1,1)=1;
 Radius = fftshift(Radius);
 %prepare a low pass filter, as proposed by Kovesi
 lowpass=1./(1+(Radius./0.45).^(30));
 %Compute the bandpass filter
 logRad=nc*exp(-log(Radius*P0).^2/lnkB2);
 RadialBP=logRad.*lowpass;
 %Compute the angle component
 Theta=atan2(Vg,Ug);
 sintheta = sin(Theta);
 costheta = cos(Theta);
 ds = sintheta * cos(Orient) - costheta * sin(Orient);
 dc = costheta * cos(Orient) + sintheta * sin(Orient);
 dtheta = abs(atan2(ds,dc));
 if ~exist('useCos','var')
     useCos=0;
 end
 if useCos==0
     %Use a gauss angle
     AngFilter=exp(-dtheta.^2/(2*gaussSD^2));
 else
     %Use a cosine angle
     AngFilter=cos(dtheta).^2.*heaviside(dtheta./pi+0.5).*heaviside(0.5-dtheta./pi);
 end
 logGaborFilt=RadialBP.*AngFilter;
 logGaborFilt=ifftshift(logGaborFilt);
 logGaborFilt(1,1)=0;
 normFact=sum(sum((abs(ifft2(logGaborFilt)))));
 RezFreq=logGaborFilt.*ImgFFt/normFact;
 RezSpace=ifft2(RezFreq);
 ReConv=real(RezSpace);
 ImConv=imag(RezSpace);
end

As for how the DC response varies with various parameters and why, a comprehensive discussion is in Boukerroui et al., equations (19) and (20) at pages 75, 76. Take a peek at Appendix C (pag 101). I would like to redirect you to the paper rather than just copy/paste the relevant paragraphs.

Hope it helps!


PSD of a sum of two stationary real signals I have 4 signals ($w_1(t),w_2(t),u_1(t),u_2(t)$ of which I know the power spectral densities (PSDs) $S_{w_1,w_1}(\omega),S_{w_2,w_2}(\omega),S_{u_1,u_1}(\omega),S_{u_2,u_2}(\omega)$ and the cross PSD $S_{w_1,w_2}(\omega),S_{w_1,u_1}(\omega),S_{w_1,u_2}$,... etc for all.

Is it possible to obtain the PSD of $S_{\Gamma,\Gamma}(\omega)$:

$\Gamma(t)=w_1(t)-w_2(t)-u_1(t)+u_2(t)$

using the PSDs that I know?

Yes, it is possible. Note that the PSD $S_{\Gamma,\Gamma}(\omega)$ is the Fourier transform of the auto-correlation function $r_{\Gamma,\Gamma}(\tau)=E\{\Gamma(t)\Gamma(t+\tau)\}$, so we can proceed by finding $r_{\Gamma,\Gamma}(\tau)$:

$$\begin{align}r_{\Gamma,\Gamma}(\tau)=&E\{[w_1(t)-w_2(t)-u_1(t)+u_2(t)]\\ &[w_1(t+\tau)-w_2(t+\tau)-u_1(t+\tau)+u_2(t+\tau)]\}\end{align}$$

Writing this out and exchanging the expectation operator with the summations gives you an expression for $r_{\Gamma,\Gamma}(\tau)$ as the sum of the auto-correlation functions of $w_1(t)$, $w_2(t)$, $u_1(t)$, and $u_2(t)$, and all possible cross-correlations between the different signals (some with positive and some with negative signs). Each cross-correlation shows up twice, e.g., you get a term $r_{u_2,w_1}(\tau)+r_{w_1,u_2}(\tau)$. Since

$$r_{y,x}(\tau)=r^*_{x,y}(-\tau)$$

you have

$$r_{x,y}(\tau)+r_{y,x}(\tau)=r_{x,y}(\tau)+r^*_{x,y}(-\tau)$$

and the corresponding cross-PSD term becomes

$$S_{x,y}(\omega)+S_{y,x}(\omega)=S_{x,y}(\omega)+S^*_{x,y}(\omega)=2\text{Re}\{S_{x,y}(\omega)\}$$

With this result you finally obtain for the PSD of $\Gamma(t)$

$$S_{\Gamma,\Gamma}(\omega)=S_{u_1,u_1}(\omega)+S_{u_2,u_2}(\omega)+S_{w_1,w_1}(\omega)+S_{w_2,w_2}(\omega)\\-2\text{Re}\{S_{u_1,u_2}(\omega)+S_{u_1,w_1}(\omega)+S_{u_2,w_2}(\omega)(\omega)+S_{w_1,w_2}(\omega)\}\\+2\text{Re}\{S_{u_1,w_2}(\omega)+S_{u_2,w_1}(\omega)\}$$


Multiplication of two impulse function $\delta(t)\cdot \delta(t)=?$ If a impulse function is mulplied with a function $f(x)$ then the formula will be apply
$$f(x)\delta(x-a)=f(a)\delta(x-a)$$ so from this formula mulplication of two impulse function will be

  $$\delta(t)\delta(t)=\delta(0)\delta(t)$$

But,what if i take two impulse function as the approximation of standard signal (rectangular or gaussian) and mulplied them directy then i will get again an impulse function of same property.

What is wrong with this two?

The property

$$f(t)\delta(t-t_0)=f(t_0)\delta(t-t_0)\tag{1}$$

is only valid for a function $f(t)$ that is continuous at $t=t_0$. Since the Dirac delta impulse $\delta(t)$ is not a function (it is a distribution) and since $\delta(t)$ is not continuous, property $(1)$ does not hold (and does not make sense) for $f(t)=\delta(t)$.

The product $\delta(t)\cdot \delta(t)$ is undefined. Also note that the quantity $\delta(0)$ is meaningless. Since $\delta(t)$ is not a function you cannot determine its value for any value of its argument.


How to correctly compute the EEG Frequency Bands with Python? so I am trying to compute the EEG(25 channels, 512 sampling rate, 248832/channel) bands (alpha, beta, gamma, etc.) with Python. I managed to do so by: firstly filtering the signal with a butterworth filter that looks like this:

def butter_bandpass_filter(data, lowcut, highcut, fs, order=2):
    nyq = 0.5 * fs
    low = lowcut /nyq
    high = highcut/nyq
    b, a = butter(order, [low, high], btype='band')
    #print(b,a)
    y = lfilter(b, a, data)
    return y

I used the filter with the low set on 0.1 and high on 80. Than I compute the fft of the signal and store it in fft1, on which I use again the butterwort filter to extract the frequencies of each band and it looks like this:

for i in np.arange(n):
    alpha1 = butter_bandpass_filter(fft1[i, :], 8.1, 12.0, 256)
    beta1 = butter_bandpass_filter(fft1[i, :], 16.0, 36.0, 256)
    gamma1 = butter_bandpass_filter(fft1[i, :], 36.1, 80, 256)
    delta1 = butter_bandpass_filter(fft1[i, :], 0.0, 4.0, 256)
    sigma1 = butter_bandpass_filter(fft1[i, :], 12.1, 16.0, 256)
    theta1 = butter_bandpass_filter(fft1[i, :], 4.1, 8.0, 256)
    sumalpha1 = sum(abs(alpha1))
    sumbeta1 = sum(abs(beta1))
    sumgamma1 = sum(abs(gamma1))
    sumdelta1 = sum(abs(delta1))
    sumsigma1 = sum(abs(sigma1))
    sumtheta1 = sum(abs(theta1))
    objects = [sumalpha1, sumbeta1, sumgamma1, sumdelta1, sumsigma1, sumtheta1]
    N = len(objects)
    ra = range(N)
    plt.title(signal_labels[i])
    plt.autoscale
    somestuffneeded = np.arange(6)
    ticks = ['alpha','beta','gamma','delta','sigma','theta']
    plt.xticks(somestuffneeded, ticks)
    plt.bar(ra, objects)
    plt.show()

The problem is that I get a high value for gamma( which represents a high cognitive activity, which the person clearly didn't have). The result for the first channel is:

Can anyone point out a better solution/ has any ideas about what is wrong/ can tell me what is wrong? Am I using the filter incorect?
Thanks!

I believe there is a much simpler way to do this with numpy.fft.rfft and numpy.fft.rfftfreq.

In the below example, I have two seconds of random data between 0.0 and 100.0 sampled at 512 Hz. I'm plotting the band amplitude below, but if you wanted power, it would just be the square of the values.

I think the reason you saw overly large values for Gamma, was that your gamma range is much larger than the others, and you're taking the sum of all the values in that range.  This is fine if you are going to compare one EEG set's gamma value to another (relative amp/power), but not if you're just looking at one set and want to compare gamma to say alpha.

If you just want to know what EEG state the person is in, you could use np.max instead of np.mean.  Then whichever band has the max will tell you about the person's state.

import numpy as np

fs = 512                                # Sampling rate (512 Hz)
data = np.random.uniform(0, 100, 1024)  # 2 sec of data b/w 0.0-100.0

# Get real amplitudes of FFT (only in postive frequencies)
fft_vals = np.absolute(np.fft.rfft(data))

# Get frequencies for amplitudes in Hz
fft_freq = np.fft.rfftfreq(len(data), 1.0/fs)

# Define EEG bands
eeg_bands = {'Delta': (0, 4),
             'Theta': (4, 8),
             'Alpha': (8, 12),
             'Beta': (12, 30),
             'Gamma': (30, 45)}

# Take the mean of the fft amplitude for each EEG band
eeg_band_fft = dict()
for band in eeg_bands:
    freq_ix = np.where((fft_freq >= eeg_bands[band][0]) &
                       (fft_freq <= eeg_bands[band][1]))[0]
    eeg_band_fft[band] = np.mean(fft_vals[freq_ix])

# Plot the data (using pandas here cause it's easy)
import pandas as pd
df = pd.DataFrame(columns=['band', 'val'])
df['band'] = eeg_bands.keys()
df['val'] = [eeg_band_fft[band] for band in eeg_bands]
ax = df.plot.bar(x='band', y='val', legend=False)
ax.set_xlabel("EEG band")
ax.set_ylabel("Mean band Amplitude")


Decoupling consonant phoneme This question follows from Designing a sound that localises well

I am developing an experimental system for pitch training.  I am giving each of the 12 pitch-classes {C CD D DE E F FG G GA A AB B } an associated consonant phoneme: {d b r ng m v sh z h l k t }

So I could sing a chromatic scale 'da ba ra nga ... etc'

My question is:  could I decouple the consonant from the vowel? So that the computer just gives the consonant as a cue?

'k', for example, would be easy: I can just 'k' into a microphone,  but others like 'ma' would be harder as it is a half vowel.

So, would there be any way to communicate "m" without suggesting pitch?

And for the less tricky consonants, like 'b',  how to separate the 'b'-ness from the 'a'-ness?  maybe I could sing 'ba' at 440Hz, and look at it on the frequency domain and just remove spectral components at 440, 880, 1320, ... ?

Or even better, I could take a slice half a second into the track which consist of pure 'aaa', convert it to the frequency domain,  and then take the beginning of the track and just subtract on the frequency domain.

Maybe I could even do this on the time domain?

Would this work with 'maaaa'?  I don't think so, I think for half-vowels I might need to follow Hilmar's answer to the linked question, and construct my own set of transients.

Maybe I could use a speech synthesiser to get my biphones 'da' 'ba' 'ra', and then work with that?

π

Consonants come in various types and classes. One major differentiation is voiced vs unvoiced (in German: "stimmhaft" vs. "stimmlos", not sure I have the proper English words here). Voiced means that the vocal cords are moving which usually also results in a pitch perception (such as "mmm"). Unvoiced means the vocal cords do not move and there is not pitch ("s", "t" for example).

It looks to me that most unvoiced consonants should do the trick. Not sure though whether your language has 12 of those.

For further classifications of consonants and vowels see for example: http://www.zompist.com/kitlong.html


Best temperature compensation equation? I'm looking for the correct temperature compensation equation to use on our project.
We are measuring the output of a detector who's signal is very sensitive to temperature drift.
Any external Temperature drift is reflected in the output.
A cycle due Day/night variation is typically what we see.
So, we apply a formula to the detector's output to try to compensate for the changes in temperature.
The formula we use is:
RITC = ((RT-AT)*TC)+I
Where:
RITC = temperature corrected output
RT = Reference Temperature
AT = Actual Temperature as Read on Instrument
TC = Temperature Coefficient
I = Input (Actual reading from the detector)
This equation works if the temperature variation is constant say (+/- 5 C) from day to day.
However, if the variation isn't constant and it changes from day to day the equation doesn't work well.
Also if there is sudden changes due to maybe air conditioning or Fans turning On/Off again the equation does not hold up.
So what is the solution?
How can the affects of temperature be removed from the detector signal?
What is the correct method to deal with sudden dynamic changes and also to deal with the slower variations?
I feel a better equation is needed.. just not sure what it is!
Best regards
Conor

Open loop temperature correction techniques are very challenging as they have to account for both static errors and dynamic temporal effects. The dynamic temporal effects refer to the thermal memory in your device under test and its interaction with the environment such that there will be response time involved from any change in the environment or its own self heating, and that response could be of multiple order.
If open loop is necessary then the static effects are best compensated with an actual measurement from the sense thermistor (with intent to compensate for any offset from that  measurement point from nominal). The dynamic effects can be best compensated with similar dynamic characterization to determine the dominant poles and create the compensator from that- the characterization can be done by sweeping a thermal sine wave (or ramp but that is mathematically more complicated) in a thermally controlled chamber starting with a very low frequency oscillation (cycle time over many minutes or more) to find the point where the output phase has shifted 45 degrees from the input (first pole) and then continue to 135 degrees to find the second pole (typically a thermal system will consist of very low frequency poles). Using a two pole open loop compensator should significantly improve an otherwise static solution.  This would be necessary if the temperature is still desired to be controlled, although open loop. If not then the direct measurement could be made and the error of the result vs temperature can be corrected from other characterization measurements- for this it is critical that the thermistor is thermally intimate with the temperature sensitive component with no additional thermal time constants in between- otherwise significant 3d finite element analysis modeling of the path between the core sensitive area and the closest measurement point may be required to resolve these same questions.
When temperature compensation is critical, a closed loop solution can be much more practical and will provide the best result, if you can spare the additional cost and power to temperature stabilize the device. For this, small ovens (a power FET or resistor can be used as a heat source) can be used when the device temperature can be maintained above the otherwise maximum temperature due to environment and self-heating, or when lower temperatures are desired, small thermo-electric coolers can be put to use. Finally where possible further insulation from the environment to make the heating or cooling most efficient can also be employed. A temperature sensor on the device combined with a PID controller makes for a highly effective solution in minimizing any temperature effects.


Uncertainty Relation between time and frequency Typically, humans can hear sound waves in the frequency range 20 Hz to 20 kHz. If one wants to make digital record of sound such that no audible information is lost, what is the longest time interval between samples that could be used?

Sampling a signal creates its spectrum replica as Nyquist Shannon theorem stated. As a consequence, to avoid aliasing, one must sample at a rate greater than twice the max frequency $B$ of the signal spectrum. Otherwise, you have aliasing as illustrated in the Figure below.

In your case, $B = 20$kHz hence sampling rate $f_s > 2B$ and the interval between samples $T_s = 1/f_s < 1/(2B)$. You can see that $T_s$ is upperbounded by $1/(2B) = 50$ microseconds.

Note that the longest $T_s$ is not 50 microseconds but $(50-\epsilon)$ microseconds  with $\epsilon$ is arbitrarily small positive quantity.


Is the discrete Gaussian kernel an eigenfunction of the DFT? So the Gaussian function is an eigenfunction of the Fourier transform because it transforms into itself, right?
But this isn't true for the sampled Gaussian in the DFT because the tails of the function are truncated, right?
Wikipedia describes a discrete Gaussian kernel here and here (solid lines), which is different from the discretely-sampled Gaussian (dashed lines):

the discrete counterpart of the continuous Gaussian in that it is the solution to the discrete diffusion equation (discrete space, continuous time), just as the continuous Gaussian is the solution to the continuous diffusion equation

Does that mean it also DFT transforms exactly into itself?  If not, is there a similar Gaussian-like function that does?

Since the DFT is representable by multiplication with the Fourier matrix, your question is equivalent to asking what the eigenvectors of the Fourier matrix are.

Actually, Wikipedia provides the answer (http://en.wikipedia.org/wiki/Discrete_Fourier_transform#Eigenvalues_and_eigenvectors).

However, since the eigenvalues ($1, -1, i, -i$) are not simple, the eigenvectors are not unique (i.e. linear combinations are also eigenvectors). Also no simple closed formula exists.

One formula for an eigenvector close to what you ask is provided by Wikipedia

$$F_m = \sum_{k = -\infty}^\infty \exp\left(-\frac{\pi\cdot(m+N\cdot k)^2}{N}\right) \quad\quad\quad m = 0,\ldots,N-1$$

So concluding, the Gaussian function itself is not an eigenvector, but an infinite sum of Gaussians.
The infinite sum can probably be interpreted as equivalent to the discretization of the frequency and time domain when we going from the FT to the DFT. So it is not as simple, as just truncating the discrete Gaussian.


Determine Frequency Content of Variable Rate Signal (Non Uniform Sampling of Signal) What is the proper way to calculate the frequency content or spectrum of a signal with a variable (but measured) sample rate?  I.e. get similar results as you would from an FFT, but for data not evenly spaced in time.

I know I can resample the data by interpolating on an evenly-spaced time vector then apply a FFT, but this seems less than ideal.

Basically what you have is data which isn't evenly spaced in time.
The conventional way to take care of the is using the Non Uniform Discrete Fourier Transform.

Another approach would be looking at the Fourier Transform as a weighted sum of the signal in time.
Since you know the time shift between the samples you can build the right "Weights" for this kind of sampling.
Once you have the weights, all you need is to calculate the weighted sum for each frequency you're interested in.


Harmonic estimation from a single note I have a recording of a single note played several times on a violin. I'm working with Matlab to estimate the amplitude of the first N harmonics. I already know which note is played but there can be fluctuations in frequency and noise (very little, actually). Can you suggest me the more robust method in this kind of situation? At the moment I am windowing the signal and using a peak detection algorithm on the FFT of each frame in order to find the amplitude of the harmonics. Is there anything better?

pitch detection and tracking.  use autocorrelation or perhaps AMDF (average-magnitude difference function) or ASDF (average-squared difference function).
either the "heterodyne oscillator" or, i would recommend, wavetable synthesis techniques to get the amplitude and phase of each harmonic.


How much noise gets through my raised cosine filter? I want to figure out how much of my total noise is getting through my filter. To see this, I generated a noise vector, applied the filter, and measured the output variance. However, the output variance changes if I scale my filter (as would be expected). I suspected that I needed to normalize the DC gain to be one, but I'm not sure if that is the right approach. Here's my code:

noiseVector = randn(1,1e6);
varOrig=var(noiseVector)

filterCoeff = rcosdesign(0.2,10,12);

filteredNoise1 = conv(noiseVector,filterCoeff);
var1=var(filteredNoise1)

filteredNoise2 = conv(noiseVector, filterCoeff/sum(filterCoeff));
var2=var(filteredNoise2)

There's a big difference between var1 ($\approx 1.00$) and var2 ($\approx 0.08$). Var1 would indicate that essentially all of the noise power was getting through the filter, but var2 would say that most of the noise is filtered out.

I also plotted the response of the unnormalized filter

 and of the normalized filter.

As expected, the DC gain of the normalized filter is 1 (0 dB).

For my purpose of characterizing how much noise is getting through the filter, should I be using the normalized one (and is that normalization factor the right one)?

For my purpose, the correct method should normalize the filter so that the sum of the filter coefficients is 1. The DC value is normalized because this is a low-pass filter. It would be different if the filter were something else, like pass-band.

Here's the code and plots that convinced me.

noiseVector = randn(1,1e5);
[Pxx_noiseOrig, W]=pwelch(noiseVector,hamming(100),[],[],60e6,'centered');

filterCoeff = rcosdesign(0.2,10,12);

filteredNoise1 = conv(noiseVector,filterCoeff);
[Pxx_noise1, ~]=pwelch(filteredNoise1,hamming(100),[],[],60e6,'centered');

filteredNoise2 = conv(noiseVector, filterCoeff/sum(filterCoeff));
[Pxx_noise2, ~]=pwelch(filteredNoise2,hamming(100),[],[],60e6,'centered');

figure('Units','inches')
plot(W/1e6,10*log10([Pxx_noiseOrig, Pxx_noise1, Pxx_noise2]),'Linewidth',2)
grid on
set(gca,'FontSize',12)
legend({'Unfiltered','No scaling','Scaled'},'FontSize',12,'Location','SouthEast')
xlabel('Frequency (MHz)','FontSize',12)
ylabel('Power (dBW)','FontSize',12)
title('Comparison of ouput PSDs','FontSize',14)
get(gcf,'Position');
set(gcf,'Position',[pos(1) pos(2) 5 3])
print('~/Downloads/filteredNoiseComparison','-dpng','-r400')

The plot shows that the original filter had a gain in the passband. After scaling, the noise power in the passband stays the same. My goal was to see how much noise was getting through the filter, so scaled version was the one to use.


Fine-grained evaluation of speech recognition results I know that usually ASR systems are evaluated using WER (word error rate), and I was wondering if it makes sense to use a more fine grained approach. E.g. transforming both the ASR hypothesis and the reference transcription to its phonetic representation (e.g. using an algorithm such as soundex or metaphone), and then taking the WER.

Some half-baked ideas for use cases:

perhaps this evaluation method can help to distinguish between errors that stem from a bad language model to those that stem from a bad acoustic model
when using third-party ASR APIs, where the there's no control over the language model, perhaps it can be used to pick the best hypothesis among the top results (especially when the speech contains non-standard names and other out-of-vocabulary) words

You have misconception that errors are caused by either language model or by acoustic model. There is no clear distinction because models are combined together in during search for best hypothesis. Interaction is complex and another major player is beam search pruning.

Most errors are originated from the acoustic model when it fails to give proper estimate for a certain frame, but the problem is that they propagate over next few words. For example, if you have an unknown word, it will not be recognized properly but the next several words will be wrong too because search state will not match the actual content. And phone-by-phone comparison will not be really meaningful.

So it's a bit more complex.


Group delay of $H(\omega)= 1- re^{j \theta}e^{ - j \omega} $ I'm studying chapter 5 of Discrete-Time Signal Processing 3rd edition by Alan Oppenheim and I'm having serious difficulties understanding how he obtained equation 5.57. For those who don't have this book I tell you that in this part it is analyzing the frequency, phase and group delay of

$$ 1- re^{j  \theta} e^{-j \omega} $$

which could either be a pole or a zero depending on whether this factor is in the denominator or numerator of the frequency response. Here r is a random magnitude variable and theta is a random phase variable

So far I have been able to understand how is the phase expression obtained as it is defined as $$\arctan  \frac {\Im}{\Re}$$

and

$$ 1- re^{j  \theta} e^{-j \omega} = 1- r[ \cos (\theta - \omega) + j \sin(\theta - \omega) ]
$$

and as cosine is even and sine is odd

$$ = 1- r[ \cos (\omega- \theta) - j \sin(\omega- \theta) ]
= 1- r\cos (\omega- \theta) + j r\sin(\omega- \theta)  $$

the resulting phase expression is:

$$ \arctan  \frac {r \sin(\omega- \theta)}{1- r \cos (\omega- \theta)}$$

which coincides with equation 5.56 in the book

But when it comes to finding the group delay (which is the negative derivative of this expression) I'm not obtaining what it says in the book. Moreover I introduced the expression in Matlab and I'm obtaining the following result:

According to the book the group delay is;

$$ \frac{r^2 - r\cos (\omega- \theta)}{1 + r^2 - 2r\cos (\omega- \theta)}$$

How did they get there? Can you help me?

The book's formula is right.

Let $$H(w) = 1 - r e^{j(\theta - w)} = [1-r \cos(\theta - w)] + j [-r \sin(\theta - w)]$$ Since the group delay $\tau$ is the negative of the derivative of the phase of $H(w)$, we first define the phase as:
$$\phi(w) = \tan^{-1}\left( \frac{-r \sin(\theta - w)}{1-r \cos(\theta - w)} \right)$$

Using the derivative rule for the inverse tangent as:
$$\frac {d}{dw} \left( \tan^{-1}(u)\right) = \frac{u'}{1+u^2}$$

We shall compute the derivative of the phase as;
$$ \phi^{'}(w) = \frac{ \left( \frac{-r \sin(\theta - w)}{1-r \cos(\theta - w)} \right)^{'} }{ 1 + \left( \frac{-r \sin(\theta - w)}{1-r \cos(\theta - w)} \right)^2 }$$

which continues as:

$$\begin{align}
\phi^{'}(w) &= \frac{ r \cos(\theta-w) \cdot[1 -r \cos(\theta-w)] + r \sin(\theta-w) \cdot [-r \sin(\theta-w)]}{(1-r \cos(\theta-w))^2 \cdot \left( 1 + \left( \frac{-r \sin(\theta - w)}{1-r \cos(\theta - w)} \right)^2 \right) }\\
\phi^{'}(w) &= \frac{ r \cos(\theta-w) -r^2 \cos^2(\theta-w) - r^2 \sin^2(\theta-w) }{(1-r \cos(\theta-w))^2 \cdot \left( 1 + \left( \frac{-r \sin(\theta - w)}{1-r \cos(\theta - w)} \right)^2 \right) }\\
\phi^{'}(w) &= \frac{ r \cos(\theta-w) -r^2 } { (1-r \cos(\theta-w))^2 \cdot \left( \frac{ (1-r \cos(\theta-w))^2 + r^2 \sin^2(\theta - w)} { (1-r \cos(\theta-w))^2} \right) }\\
\phi^{'}(w) &= \frac{ r \cos(\theta-w) -r^2 } { (1-r \cos(\theta-w))^2 + r^2 \sin^2(\theta - w)}\\
\phi^{'}(w) &= \frac{ r \cos(\theta-w) -r^2 } { 1 - 2r\cos(\theta-w) + r^2 \cos^2(\theta - w) + r^2 \sin^2(\theta - w)}\\
\end{align}
$$

Which finally simplifes to
$$ \phi^{'}(w) = \frac{ r \cos(\theta-w) -r^2 } { 1-2r \cos(\theta-w) + r^2}$$

Then since the group delay is the negative of this we finally have:
$$ \tau(w) = -\phi^{'}(w) = \frac{ r^2 - r \cos(\theta-w) } { 1 + r^2 -2r \cos(\theta-w) }$$

$$ \tau(w) = -\phi^{'}(w) = \frac{ r^2 - r \cos(w-\theta) } { 1 + r^2 -2r \cos(w -\theta) }$$

Which is the same as the book's formula.

Coming to Matlab's output, you may try simplifications. But they are known to produce more verbose outputs than possible with simplifications.


How to Realize the Sigma Point Sampling Function in Unscented Kalman Filter? Recently I'm learning the unscented kalman filter (UKF). When designing the unscented kalman filter, it involves a non-linear function to generate the sigma points and then use the system non-linear function to transform those sigma points.
However, I'm still confused about why we need to select sigma points to compute the estimate value and how the sigma points selecting function affects the efficiency of UKF.

Here is the pdf book I'm learning now:

https://drive.google.com/file/d/0By_SW19c1BfhSVFzNHc0SjduNzg/view

In this book, it mentions the sigma function in Chapter 10.3

In any Kalman Filter one need to calculate the 1st and 2nd moment of the data under the transformation.

The image above taken from The Unscented Kalman Filter for Nonlinear Estimation
by Eric A. Wan and Rudolph van der Merwe

The problem is that there are some transformations which their linearization (As done in the Extended Kalman Filter - EKF) yield the wrong evaluation of the mean and covariance.
Then the Unscented Transform comes to help and able, using only few data points from the data prior to the transformation to estimate those moments pretty good in many cases.

Why does it work?
One way to look at it is through the analysis done in Cubature Kalman Filters (CKF) which is a generalization of the UKF.
The problem there is solving high dimensional integrals weighted by Gaussian Kernels.
This problem can be solved using the Gaussian Quadrature method which is a numerical method to solve high dimensional integrals.
The problem is basically done by selecting pre defined points in the data and their weights and summing them to calculate the integral.
This is just what's done in the UKF (unscented Transform to be exact).


How to run Welch’s method on cross-spectral density I have a multivariate time series and I’m trying to figure out the cross spectral density of all pairs of variables efficiently. I’m unsure of how to implement Welch’s method in this context. Do I calculate the PSD of each series individually using Welch and then use those to compute the cross-spectral density or do I average the cross-spectral density of each sub-band?
If the answer is to average the cross-spectral density of each sub-bad, then how does one write an efficient algorithm to do this for a high dimension?
Thanks!

Let's define Cross-Spectral Density; you can do that via the usual "Fourier transform of the cross-correlation", but assuming WSS, Wiener-Khinchin allows us to state it as
$$R_{XY}(f) = \lim_{\Delta t \to \infty} {\large\mathbb E} \left\{\frac{X_{[t_0;t_o+\Delta t]}^*(f)Y_{[t_0;t_o+\Delta t]}(f)}{\Delta t}\right\}$$
(abusing the subscript here to denote this is a finite-duration "observation" from point $t_0$ to $t_0+\Delta t$).
Now, you could of course first estimate the cross-correlation and then estimate the spectrum of that, and end up with an estimate of the CSD; however, estimation of a cross-correlation is usually done by first transforming the both signals to frequency domain, (conjugate) multiplying there, and then transforming back.
Transforming that back into frequency domain seems like you've basically reverted the last step.
Small problem with using Welch here: it usually ends with squaring the transform of your averaged segments – that erases the phase information necessary in the formula above.