Tracking a sine wave with noise This will probably be an extremely simple question for some one with any background in signal processing(not my background)

Let say I have signal $$x(t)=A\sin(\omega t)$$ where A is known and $\omega$ is unknown. Let say I only observe a noisy signal of $x(t)$ at discrete time periods $t=1,2,3....$. My signal $y(n)$ has the following form
$$y(n)=x(n)+\epsilon$$ where $\epsilon\sim N(0,\sigma^2)$.
My question is , what would be the best estimate of $y(n+1)$ given the observations till time time period $n$.

That's actually not a trivial problem and the best approach really depends on your specific requirements.

My first approach would be a phase-locked loop and simply use the local oscillator as the predictor (with proper amplitude scaling). This can give you extreme accuracy if you are willing to pay with time.


Signals and Filters 101 I am a complete novice to the world of signal processing, and I am an environmental engineering grad by degree (i.e. zero background in signals/systems). We acquired a couple sensors and a data recorder at work and I am required to analyze the data we measure. Pertaining to that, I have two questions:

Is anyone here familiar with Engineering Data Management software by Crystal Instruments?
I just need to learn the basic nuts and bolts of signals without having to go in depth into the math of deriving the Fourier Transform or signal design etc. Is that a possibility, and what is a good resource for this?

I just need to learn the basic nuts and bolts of signals without having to go in depth into the math of deriving the Fourier Transform or signal design etc. Is that a possibility, and what is a good resource for this?

Sorry, I don't think it works this way. Without a proper understanding of the fundamentals, it's unlikely what you will every get useful results of your analysis. Check https://ccrma.stanford.edu/~jos/filters/ or http://www.dspguide.com/ for starters. This is good too, if you are up for it https://ocw.mit.edu/resources/res-6-007-signals-and-systems-spring-2011/introduction/


How can I plot a 3D graph of a given Laplace Transform of a function? Let's say I have a function called $f(t)$ in time domain as:

$$f(t) = \exp(-3t)\cos(5t)$$

And the Laplace transform of this function call it $F(s)$ becomes:

$$F(s)=\frac{(s + 3)}{(s + 3)^2 + 25}$$

I want to plot the 3D plot of $\lvert F(s)\rvert$ as a surface above the $s$-plane.

I couldn't find any script. I have the version R2014a of MATLAB.
How is this done in MATLAB? Something similar to this plot:


Why does raised cosine pulse through Butterworth filter introduce (significant) ISI? Here is a simple MATLAB experiment:

alpha = 1;
R = 10;
[b,a] = butter(7,alpha*1/R);
x = resample(sign(randn(100,1)),R,1,0); % zero-order-hold
y = filter(b,a,x);
eyediagram(y, 2*R);

With channel twice the bandwidth than the pulses (alpha=2) the eye diagram looks clean. The first thing I am concerned that the eye diagram at the critical bandwidth (alpha=1, i.e. channel bandwidth is bitrate/2) still looks reasonably clean (little bit ISI but not too much):

I would have expected significant ISI at the critical bandwidth.

Now instead of using rectangular pulses, I use a raised cosine filter:

 hcos = rcosdesign(0, 30, R);
 x = upfirdn(sign(randn(1000,1)), hcos, R);

Now the eye diagram is significantly worse:

I would expect the result to be significantly better. I played with all possible rolloff factors 0...1 and number of symbols. I also tried cutting the signal at beginning/end, thinking it may be the initial transient. No changes.

What is going on here?

(Additional info: I do not want to simulate TX/RX with raised cosine but assess a channel without having to use an equalizer - in simulation. One way to do this is using an equalizer or just using a pulse shape that works close to the theoretical bandwidth - that should be raised cosine with a low beta).

That's because Matlab's rcosdesign creates a root-raised-cos filter, not a raised cos filter. If you want no ISI (ie meet Nyquist's criterion), you need to do the match filtering of the reception part :

%%
alpha = 1;
R = 10;
hcos = rcosdesign(0.25, 30, R);
x = upfirdn(sign(randn(1000,1)), hcos, R);
y = upfirdn(x, hcos, 1, R);
eyediagram(x, 2*R);
hold on
eyediagram(y, 2*R);

eye of x

eye of y


Processing using Matlab I'm new to Matlab.As a beginner, I want to work on few signals which are used in day today life and extract useful information through it. Can someone please help me from where I could start or share a blog which might help me with this. Thank you.

Your question is a bit broad. Since you probably have Matlab (I hope), Matlab help is where to start, for instance Measurements and Feature Extraction. From there, you can start with signals available in Matlab, get some tutorial code, etc.


What Are Intuitive Explanations for Shrinkage in the Context of Image Denoising? I have occasionally come upon the term "shrinkage", mostly in the context of denoising methods.
My rough understanding is that it refers to the part where the real distribution might not be exactly like the theoretical one that's being modeled while still being relatively similar (e.g. mostly Gaussian with some outliers), and to the fact that "shrinkage factors" are introduced/estimated in statistical estimators to counter that and have a "better" estimation (for a given definition of "better" that's probably not always the same depending on your goal).
The Wikipedia page gives out the example of the sample/population variances and hints at the fact that using $\frac{1}{n-1}$ or $\frac{1}{n}$ depends on the shape of the distribution, but it's hard to have a good idea of why I should choose one or the other.
How can I put that more intuitively? Can you point me towards other examples that may help me form a bigger picture of this concept?

In the image processing field shrinkage means we reduce the power of some components of the signal.
The idea is to have a basis which the signal can represented in in a sparse manner.
For instance, if the signal is harmonic, efficient basis would be the Fourier.
So if we look at harmonic signals and noise in Fourier we expect to see dominant component and other components are the noise, which will be shrinked in order to apply Denoising.
This is the idea behind what's known as Sparse Representations.
In the context of Image Denoising you may have a look at Image Denoising with Shrinkage and Redundant Representations.
In the general context you may have a look at the course Sparse Representations in Signal and Image Processing: Fundamentals.


Adaptive equalizer for BPSK signals with frequency offset I'm trying to equalize a BPSK signal which has a residual frequency offset on it. This frequency offset is dealt with at a later stage in the DSP chain, along with phase estimation.

I have implemented a FIR filter which updates its filter taps using gradient descent CMA:

The equalized signal (y) is a discrete-time convolution of the filter taps (W) and the input signal (X):

$$
y = W^HX
$$

The update error is:
$$
\epsilon = |y|^2-1
$$

Finally I compute the update as:
$$
W_{i+1} = W_i - \mu \epsilon y^*X
$$
Where $\mu$ is the step size = 0.001.

My problem is that due to the frequency offset, my input signal X is complex, which is resulting in complex filter taps in W. This in turn introduces a huge phase error into my signal. One possible solution I've looked at so far is putting the equalizer behind the frequency offset compensation stage. However, since the signal still contains some phase error, it's not entirely constrained to the real axis in the constellation diagram and I have the same problem. My workaround so far has been to update the filter taps as:
$$
W_{i+1} = W_i - real(\mu \epsilon y^*X)
$$
My reasoning being that after frequency offset compensation and phase estimation, there should only be noise in the quadrature part of the signal, and I can safely discard it.

I'm unhappy with this solution for 2 reasons:

I'm not entirely sure it's theoretically sound, and at the very least it's rather unorthodox.
I would prefer to equalize the signal upstream in the DSP chain.

I would appreciate any help or insights!

It sounds like your biggest issue is dealing with a phase offset and otherwise you have no other concerns with your equalization approach?  In that case, the issue is resolved by simply having a phase correction prior to making your decisions on the real signal. If done correctly your carrier recovery implementation should remove any residual long term phase offset. This post may help give you further insight: Phase synchronization in BPSK.  I duplicated the final block diagram below where in your case "I data demod" would be your BPSK output. I want to point out using this diagram, that any residual phase error would accumulate in the PI Loop Filter which would then adjust the NCO to rotate the signal to have zero error (on average). Note that the Loop as shown if completely blind will lock in quadrature or 180° phase positions - if that is your remaining issue then another question on carrier recovery implementations would be more appropriate specific to any issue you have there, but a simple solution is to derive the phase error itself from $\phi \propto Q \text{sign}(I)$ since the objective in this case is to drive any signal on Q down to the noise. The assumption is the signal has been AGC'd prior to this point so such an approach to phase detection would result in a consistent loop gain.  Here is another related post with regards to what gain constants you should use (what loop BW should be used): Loop bandwidth for symbol timing recovery

Also notice the similarity of the above diagram to a classical Costas Loop with the following figure from http://michaelgellis.tripod.com/mixerscom.html. This is because the Decision Directed CR Loop is a Costas Loop. The referenced link also diagrams the classical Costas Loop for BPSK specifically, which also shows how the digital implementation such as shown in the diagram above can be modified for the BPSK only case (which is modified the phase error metric as I described earlier).  This solution will still have a 0/180° phase ambiguity which is typically resolved with minimum knowledge of the message (header or some other recognized feature in the data message or encoding in the data message).


Is there a physical reason why phase should be uniformly distributed? I was thinking a little about how signal phase might be distributed. Usually in such discussions I see people make the assumption that phase is uniformly distributed between $0$ and $2 \pi$. This only makes sense to me if its a statement about our ignorance of the system in question. That is, if we don't know anything about the phase, then it makes sense to assume all phase values are equally as likely -- uniformly distributed.
Is this the only reason though? Is there a physical reason, or process, as to why this can be the case?

A zero mean complex white noise process has a uniform phase distribution.  The uniform phase distribution would occur for any complex noise process whose I and Q (real and imaginary) components are zero-mean, independent and identically distributed.
The thermal noise floor common to radio receivers for receiving signals without other specific interference is typically well approximated by a zero-mean complex white noise Gaussian distribution; meaning a white noise Gaussian distribution on I (the real component) and Q (the imaginary component) and these two are independent from each other. This makes sense if you consider that the real and imaginary values for any given sample have no relation to each other and are also therefore independent. Within this restriction of having an independent and identically distributed random process for I and Q, the actual phase of any sample is equally likely as any other phase and is therefore uniformly distributed (regardless if the system is white noise or band-limited).
Complex quantization noise (using two real A/D converters in practice; one for the I channel and the other for the Q channel) will also be well approximated by a white noise process where in this case the phase is uniformly distributed as above but unlike the thermal noise example given above where the signed amplitude was approximated as Gaussian distributed, the amplitude for quantization noise is also uniformly distributed over a quantization level.


Deconvolving known impulse response from a sampled noisy signal I am interested in measuring a signal with significant energy content up to 1-2 kHz. However, I am only able to sample the signal after it has passed through a first-order lowpass filter of known bandwidth (which significantly distorts the signal; in this case the bandwidth is 500 Hz). I sample the lowpass-filtered version at 5 kHz, and my measurements are noisy.

What are some options I can use to try to recover the original signal after taking noisy measurements of a lowpass-filtered version of it?

If it helps, I know the shape of the signal I am looking for; I am just trying to measure how it is scaled in amplitude and in duration.

This very, very much sounds like the classical problem for which you can derive that the matched filter is the optimal receive filter to maximize SNR.

You say you know the shape of the signal; I interpret this like your system looks like this

TX impulse $p$ -> Pulse shaping filter $g$ -> 500 Hz low pass $h$ -> +Noise $z$ -> Sampling

That implies that the RX signal (given we sample sufficiently fast) is the discrete

$$r[n] = p*g*h + z$$
with $*$ being the convolution operator.

Since you know the order and bandwidth of $h$, I'd argue you actually know what $h$ is precisely. Thus, you can "pre-compute" the combination

$$m=g*h$$

The receive filter that maximizes the SNR is the matched filter – sadly, the matched filter to an IIR is not realizable (when does the time-inverse of something without end start?), you can't directly work with your actual system model.

However, if you can find a FIR filter that approximate your $h$ well enough, you should be able to  bound the error, and still get a reasonably good amplitude/power estimate by calculating the time-inverse conjugate (${}^\dagger$ operator)
$${\hat m}^\dagger =\left(g*\hat h\right)^\dagger$$
 based on the approximative filter $\hat h$. You'd then have
 $$
\begin{align}
\hat r &= p*g*h*\hat m^\dagger &+ z*\hat m^\dagger\\
&= p*m*\hat m^\dagger  &+ z*\hat m^\dagger
\end{align}$$

Since the noise process $z$ should not be correlated to the channel/filter, the second term should become small for sufficient observation.


determining FEC or parity bits from percent error I am very new to signal processing and my background is in physics.  I would like to know if it is possible to determine the number of parity bits needs to theoretically get 100% transmission from a transmission channel with known percent error.  Over a given transmission 25-30% of the transmission is incorrect.

Is it possible to know with say a reed-Solomon technique how many paity bits I would need to add to each block to move in the neighborhoos of 100%?

I am processing my signal using MATLAB right now...

Is it possible to know with say a reed-Solomon technique how many paity bits I would need to add to each block to move in the neighborhoos of 100%?

Normally yes, if you know that you will never have more than $e$ bit errors in a block of length $N$, then you could design a RS code that would fix them all.  The problem in your case is that 30% bit errors is a LOT of errors.  Assuming the errors are randomly distributed, no RS code will be able to fix that many errors.  The number of parity symbols you would need would exceed $N$.

The analysis goes as follows- RS codes work on symbols, not bits.  All of the RS codes that I have seen have symbol widths of 6 to 9 bits, with the most common symbol width being 8 (CD's use an 8-bit wide RS code).  At 30% bit errors, well over half of the symbols will have errors.  You need two parity symbols for each error that you want to correct, so when over half of the symbols have errors it is impossible to correct them.

You may think "why not make the symbols narrower, reducing the chance that they will contain a bit error?"  The problem is that the block length, $N$, is $2^m-1$, where $m$ is the RS symbol width.  If you reduce the symbol width to 1, then your block length is 1 and you can't correct any data.  Thus, no RS code can correct your errors.

I would say that your best bet is a low-rate Turbo code, though turbo codes cannot guarantee that they will fix your errors.


What is the difference between 'de-noising' and what we commonly know as filtering methods to remove noise? Pretty simple question - I am trying to figure out what exactly is different between 'de-noising' a signal, and simply filtering it (as we commonly know) to remove noise. Is this a case of lexical overlap or is there something fundamentally different? Why is it called 'de-noising'?

Edit: Perhaps crucially, when we talk of filtering a signal to maximize its SNR, we usually mean AWGN in the colloquial context. So is the 'noise' being referred to in de-noising also AWGN, and if so, is de-noising simply a different way of removing it, or is it a different type of noise (non-gaussian, colored, etc) to begin with?

De-noising is about the goal, and filtering is about the technique you employ.

You can obviously de-noise via filtering. For example, if you know that your system cannot transmit frequencies above a certain threshold, you can apply a low-pass filter. However, you can de-noise by other techniques as well, such as by averaging multiple recordings of a signal. Similarly, you can apply filters with other goals than (just) noise reduction.

Note that "noise" may not necessarily be additive random noise - depending on your context, non-pertinent signals may be lumped with noise as well, such as the sound of vuvuzelas.


Addition of two signals in polar coordinates There is quadrature detector which receive radio signal and return its amplitude and instantaneous phase after performing hilbert transformation and transferring the signal to zero frequency. The main purpose of this device is save information about phase of received signal after digitizing (here is example of such receiver). Simplified structure of such detector shown below:

The following octave code try emulating the detector (I'm not sure that is correct):

% Prepare to work.
close all;
clear all;
clc;

% Main constant.
frequency_local_oscillator = 30e6;
frequency_intermediate     = frequency_local_oscillator;
frequency_sampling         = frequency_intermediate * 4;
time_begin                 = 0;
time_end                   = 10e-4;
time                       = linspace (time_begin, time_end, frequency_sampling * (time_end - time_begin) );
size                       = length (time);

% Pulse.
magnitude_pulse            = 1.0;
phase_pulse                = 0;
carrier_pulse              = magnitude_pulse * cos (2 * pi * frequency_intermediate .* time + phase_pulse);
timeshift_pulse            = time_begin + (time_end - time_begin) / 4;
duration_pulse             = (time_end - time_begin) / 2;
envelope_pulse = zeros (1, size);
envelope_pulse = envelope_pulse + heaviside (time - timeshift_pulse) .* heaviside (duration_pulse - time + timeshift_pulse);
signal_pulse   = carrier_pulse .* envelope_pulse;
if 0 % Make signal more realistic.
% Noise.
magnitude_noise            = 0.01;
carrier_noise              = magnitude_noise * randn (1, size);
signal                     = signal_pulse + carrier_noise;
else
signal                     = signal_pulse;
end

% Moving to zero frequency and and quadrature detection.
signal_i =  signal .* cos (2 * pi * frequency_local_oscillator * time);
signal_q = -signal .* sin (2 * pi * frequency_local_oscillator * time);

% LPF.
[b a] = ellip (4, 0.1, 60, 0.1);
signal_i_filtered = filter (b, a, signal_i);
signal_q_filtered = filter (b, a, signal_q);

% Decimation (performed at two stages with common ratio: 15 х 16 = 240).
signal_i_decimated = decimate (signal_i_filtered,  15);
signal_q_decimated = decimate (signal_q_filtered,  15);
signal_i_decimated = decimate (signal_i_decimated, 16);
signal_q_decimated = decimate (signal_q_decimated, 16);
% Time scale after decimation.
time_decimated = linspace (time_begin, time_end, length (signal_i_decimated) );

% Magnitude and instantaneous phase.
signal_magnitude = sqrt  (signal_i_decimated .^ 2 + signal_q_decimated .^ 2);
signal_phase     = atan2 (signal_q_decimated, signal_i_decimated);

% Display signal.
figure;
grid on;
subplot (2, 1, 1);
plot (time, signal, 'b');
hold on;
plot (time, envelope_pulse, 'r');
title ("Input signal.");
xlabel ("Time [s].");
ylabel ("Magnitude [V].");
subplot (2, 2, 3);
plot (signal_magnitude, 'b');
title ("Magnitude of the signal.");
xlabel ("Sample number.");
ylabel ("Magnitude [V].");
subplot (2, 2, 4);
plot (signal_phase, 'g');
title ("Phase of the signal.");
set (gca, 'YTick',     -pi: pi / 2: pi)
set (gca, 'YTickLabel',{'-pi', '-pi / 2', '0', 'pi / 2','pi'})
xlabel ("Sample number.");
ylabel ("Phase [rad].");

Here is result of the code which looks in my opinion more or less real:

Question:
I have records of two different signals from output of the detector (desired signal and noise). And now I'd like to adding those signals (combine as if they was received simultaneously). Is it possible to do it and how can I do this?
The main purpose of this task is performing measurements of the the once recorded signal with different calibrated noises.

My attempts:
About possibility I think it's possible if all stages of the detector are linear. If it is true this allows swap processing stages.

In my opinion add together magnitudes and phases of signals is wrong. I think that this task can be solved with the following steps:

revert conversion at I and Q;
get sum of components of two signals;
convert I and Q to magnitude and phase again.

Am I right?

Mixing is not the correct terminology for what you want to do.  Mixing is multiplication and you want additive noise.  Convert back to I/Q, add complex noise, then convert back to mag/phase!


How I can perform SIMO channel convolution in MATLAB In OFDM system, I need to transmit a signal $x$ in SIMO channel $H$ of  single Tx antenna and $4$ Rx antennas. The initial equation for that convolution is $r = Hx + n$ where $n$ is the noise.

As known, in OFDM, after serial/Parallel conversion with dimension $(Number of subcarriers, Number of symbols)$, we perform the $IFFT$ and then add Cyclic prefix. Than, P/S conversion in order to do convolution between signal and channel. In case if we are using SISO system, we perform in $MATLAB$ the command $conv(H,x)$, regardless the dimensions in length of channel and transmitted signal,  but the problem I'm facing now, what's about if we are using SIMO channel? the convolution operation should be replaced by product since we have matrices, as mentioned,  $r = Hx$

Below is the code, the issue is in the last line, how can I get the received signal $r$ in case if $H$ and $x$ don't have similar dimensions to do the matrices product operation?

clear all;
clc;

    N = 128;       %Number of subcarriers
    P = 32;        % CP length
    Q = N + P;    % subcarriers + CP
    AA = 1000;     %number of symbole
    M = 15;        % IR order
    p = 4;         % # of Rx antennas

    sn=sign(randn(N,AA))+1i*sign(randn(N,AA));     %generate a signal ()
    xn=ifft(sn);   %    FFT
    cp(1:P,1:AA)=xn(N-P+1:N,1:AA);                  % Get CP
    x=[cp;xn];                                      %Add cp into signal
    H = (randn(p,M+1)+1i*randn(p,M+1));            % randum channle (4x16)
    r = H*x;                                       %% The received signal should be conv(H,x)

Your question is specific to how to do the convolution with a matrix multiplication in the case of a SIMO system.

Below is a graphic from Compensating Loudspeaker frequency response in an audio signal showing how the Toeplitz matrix is used to do SISO convolution with matrix multiplication:

In your case "r" in the graphic is "x" in your code, and the Toeplitz Matrix A is formed by placing x in the column and shifting this along the diagonal over the number of columns given by your channel length (in your case 12). So if your transmit signal was 160 samples long and your channel 16 long, A would be a 175 x 16 matrix. As shown in the figure above, if you only had a single channel, the channel would be multiplied as a column vector resulting in the matrix product of 175x16 with 16x1 resulting in a solution as a column vector of 175 x 1. (Remember to do the matrix multiplication, the inside dimensions must match and the result will be the outside dimensions).

For multiple channels as your case, simply multiply by the 16x4 channel matrix instead. This will result in 175x16 mulitplied with 16x4 resulting in a matrix result of 175 x 4, representing a column for each received signal.


How can Recognize handwriting number Which algorithms should use to recognize handwriting numbers from an image?

Depends on where you're trying to implement/the resources available to you.

I've used Support Vector Machine learning for that task; you have to feed in individual characters to train your algo.

take a look at python's scikit-learn and if you ned to implement on FPGA check out: MyHDL


Steps for proving Tweedie's Formula Assuming that $f_{X|D}(x|d)=e^{\alpha d x - \psi(d)}f_0(x)$ for some constant $\alpha$ such that it is a valid pdf for every value of $d$, i want to establish that:
     $$\mathbb{E}[D|X=x]=\frac{1}{\alpha}\frac{\partial}{\partial x}\lambda(x)\quad \text{, where}\quad  \lambda(x)=\log\frac{f_X(x)}{f_0(x)}$$

*There's no missing information. It's from a collection of Information Theory problems.

We can compute the marginal PDF of $X$ as follows:
\begin{align}
 f_X(x)=\int_\delta f_{X|D}(x|\delta)f_D(\delta)d\delta = f_0(x) \int_\delta e^{\alpha \delta x - \psi(\delta)} f_D(\delta)d\delta
\end{align}

We can express the previous equation in terms of $\lambda(x)$ by dividing both sides by $f_0(x)$ and applying the natural logarithm.

\begin{align}
 \lambda(x) := \log \frac{f_X(x)}{f_0(x)} = \log \left[ \int_\delta e^{\alpha \delta x - \psi(\delta)} f_D(\delta)d\delta\right]
\end{align}

Now, if we derive w.r.t the parameter $x$, we obtain the following equality:

\begin{align}
 \frac{\partial}{\partial x} \lambda (x) = \frac{\int_\delta \alpha \delta e^{\alpha \delta x - \psi(\delta)} f_D(\delta)d\delta}{\int_\delta e^{\alpha \delta x - \psi(\delta)} f_D(\delta)d\delta}= \alpha \frac{\int_\delta \delta e^{\alpha \delta x - \psi(\delta)} f_D(\delta)d\delta}{\int_\delta e^{\alpha \delta x - \psi(\delta)} f_D(\delta)d\delta}
\end{align}

With this equality, we are now able to express the MMSE estimator of $D$ in terms of $\lambda(x)$:

\begin{align}
 \mathbb{E}[D|X=x]&=\int_{\delta}\delta f_{D|X}(\delta|x)d\delta=\int_{\delta}\delta\frac{f_{X|D}(x|\delta)f_D (\delta)}{f_X(x)}d\delta=\frac{\int_\delta \delta f_0(x)e^{\alpha \delta x - \psi(\delta)}f_D(\delta) d\delta}{f_X(x)} = \\
 &= \frac{\int_\delta \delta f_0(x)e^{\alpha \delta x - \psi(\delta)}f_D(\delta) d\delta}{\int_\delta f_0(x) e^{\alpha \delta x - \psi(\delta)} f_D(\delta)d\delta}=\frac{1}{\alpha}\frac{\partial}{\partial x} \lambda (x)
\end{align}


Block diagram transfer function of a line What is the transfer function of a line, can i ignore it?

After simplification i got the diagram:

I don't know what to do with the empty feedback

The transfer function of a line is 1.

And hence you can reduce your blocks into a single block with transfer function
$$\frac{C(s)}{R(s)} = \frac{G(s)}{1+G(s)}$$

where, $$G(s) = \frac{K}{s(s+2)(s+3)+s}$$


position of useful filtered data inside DFT based filtered output I'm trying to filter a 400 samples signal with various bandstop FIR filters (constant group delay), one at a time to see which one gives me the desired result. Every filter was build using the Kaiser Window method in Matlab. Each filter has a different number of coefficients. The smaller one has 782 and the biggest one 2460. I'm performing the filtering via fast convolution using the following Matlab code:

temp = ifft( fft(filterCoeffs,ffsSize)' .* fft( signalSamples, ffsSize) );

FilteredSignal = temp(offset+1:offset+length(signalSamples));

I'm trying to identify where inside the FilteredSignal vector of ffsSize are the 400 filtered samples. So far depending on the value of ffsSize and the amount of filter coefficients is the offset that I have to impose to recover the original 400 samples filtered. Even though I already know the right offset (where the useful filtered signal starts) I want to know the theoretical explanation of why there. I also want to know why is the group delay 1/2 of the amount of filter coefficients?

By the way if I use the Matlab function:

filter(filterCoeffs,1,signalSamples)

It returns 400 samples but not the right ones. Any idea why? Is like it doesn't identify the right offset

I have 2 cases:

No time aliasing (circular convolution equivalent to linear convolution)

ffsSize = length(filterCoeffs)+length(signalSamples) - 1;
aliasingSamples = 0
offset = groupDelay(Filter)+ length(signalSamples);

Time aliasing (circular convolution not equivalent to linear convolution)

ffsSize = length(filterCoeffs);
aliasingSamples = ffsSize - length(signalSamples);

if(groupDelay(Filter) < aliasingSamples) %this only happens when I use the 782 coefficient filter
    offset = groupDelay(Filter);
else
    offset = aliasingSamples;
end

Lets put forward the intuition behind the concept of the group delay before further discussing how to find the delay of FIR filters.

Consider an input signal $x[n]$ of length $L_x$ which is nonzero between $n=0$ and $n=L_x-1$. And let a simplistic filter impulse response to be $h[n] = \delta[n-d]$. The output is immediately shown to be $y[n] = x[n-d]$ which is input shifted by $d$ samples to the right (for $d$ positive).

Let's look at the frequency response of the filter; $H(e^{j\omega}) = e^{-j \omega d}$ which is nothing but a phase term $\phi(\omega) = - \omega d$.

Consider a slightly complex filter $h[n] = 0.1 \delta[n] + 0.8\delta[n-1] + 0.1 \delta[n-2]$. Its frequency response can be written as $H(e^{j\omega}) = e^{-j 1 \omega} B(\omega)$ where $B(\omega) = 0.2 \cos(\omega n) + 0.8$ is a real valued function and has no phase-term. The output of this filter can be written in terms of delayed inputs due to each distinct impulse such as :
$$y[n] = 0.1 x[n] + 0.8 x[n-1] + 0.1 x[n-2]$$

It can be proposed that $y[n]$ is dominated by the middle term $x[n-1]$ of a shift $1$ sample. (a more detailed analysis would say that low frequencies of $x[n]$ are dominated by the middle term and high frequencies are tending to cancel etc. which makes it a lowpass filter eventually)

In this latter example the phase function of the filter was $\phi(\omega) = -1\omega$ whereas the dominating component $x[n-1]$ at the output also had a delay of one sample to the left.

This observation can be used as an intuitive basis for the claim that the delay (or the group delay) associated with an FIR filter is given by the negative of the derivative of its phase response wrt frequency $\omega$; i.e.,
  $$\boxed{ \tau = - \frac{d \phi(\omega)}{d\omega} }$$

Lets restrict the discussion to LTI, FIR symmetric filters which have linear phase responses. An LTI filter is described by its impulse response $h[n]$ which enables us to compute its response to an arbitrary input $x[n]$ via the discrete convolution sum:

$$ y[n] = x[n] \star h[n] = \sum_{k=-\infty}^{\infty} x[k] h[n-k] $$

where the dummy index $k$ would actually run over those valid range of nonzero multiplicands.

The filter also has a frequency response which is the DTFT of its impulse response:
$$ H(e^{j\omega})=  \mathcal{F} \{ h[n] \} =\sum_{n=-\infty}^{\infty} h[n]e^{-j\omega} = |H(e^{j\omega})| e^{j \phi_h(\omega)}$$

Now we define the group delay (in samples) of filter as:
$$\tau = - \frac{d\phi_h(\omega)}{d\omega}$$

In general group delay will be a function of frequency $\omega$, however linear phase FIR filters have the property that their phase function is a linear function of frequency such as $$\phi_h(\omega) = K - d \omega$$

From which we deduce that the group delay will be $\tau = d$ which is independent of frequency; i.e. all components of an applied signal $x[n]$ will shift the same amount when filtered by a linear phase filter.

Given such a linear phase filter whose delay is $\tau = d$ and assume that you apply a very narrowband signal $x[n]$ such as:
$$ x[n] = \cos(\omega_0 n) \frac{ \sin(\omega_1 n)}{\pi n} $$
where $\omega_1 \ll 1$  and $0 < \omega_0 < \pi$.  Then the associated output for such a signal will be of the form
$$ y[n]=A(\omega_0) x[n-d] = A(\omega_0) \cos(\omega_0 (n-d)) \frac{ \sin(\omega_1 (n-d))}{\pi (n-d)}$$

Finally, for example,  it can be shown that a causal and symmetric FIR filter of length $M = 2L + 1$ from $n=0$ to $n=2L$ such that $h[n] = h[2L-n]$ will have a linear-phase term of the form $\phi_h(\omega) = K - \omega L$. Where $K$ can be zero. Hence its group delay will be $L$ samples. Lets show it here:

$$
\begin{align}
H(e^{j\omega}) &= \sum_{n=0}^{2L} h[n] e^{-j \omega n} \\
&= \sum_{n=0}^{L-1} h[n] e^{-j \omega n} + \sum_{n=L+1}^{2L} h[n] e^{-j \omega n} + h[L] e^{-j \omega L} \\
&= \sum_{n=0}^{L-1} h[n] e^{-j \omega n} + \sum_{n=L+1}^{2L} h[2L-n] e^{-j \omega n} + h[L] e^{-j \omega L} \\
\end{align}
$$
Substitute $2L-n = m$ in the second sum:
$$
\begin{align}
H(e^{j\omega}) &= \sum_{n=0}^{L-1} h[n] e^{-j \omega n} + \sum_{m=L-1}^{0} h[m] e^{-j \omega (2L-m)} + h[L] e^{-j \omega L} \\
&= \sum_{n=0}^{L-1} h[n] e^{-j \omega n} + e^{-j\omega 2L} \sum_{m=L-1}^{0} h[m] e^{j \omega m} + h[L] e^{-j \omega L} \\
&= \sum_{n=0}^{L-1} h[n] e^{-j \omega n} + e^{-j\omega 2L} \sum_{m=0}^{L-1} h[m] e^{j \omega m} + h[L] e^{-j \omega L} \\
\end{align}
$$
Recognize the similarity of the first sum and second sum such that denoting the first sum as $A(\omega)$ then the second sum is $A(\omega)^*$ (assuming real $h[n]$), Hence

$$
H(e^{j\omega}) = A(\omega) +  e^{-j\omega 2L} A(\omega)^* + h[L] e^{-j \omega L} \\
$$

Taking all of them into the $e^{-j\omega L}$ parenthesis yields:
$$
H(e^{j\omega}) = e^{-j \omega L} \left( A(\omega)e^{j \omega L} +  A(\omega)^* e^{-j \omega L} + h[L] \right) \\
$$

Since $A(\omega)e^{j \omega L} +  A(\omega)^* e^{-j \omega L} = 2 \mathcal{Re} \{ A(\omega)e^{j \omega L} \}$

We reduce the frequency response of the FIR filter to:
$$H(e^{j\omega}) =  e^{-j \omega L} ( 2 \mathcal{Re} \{ A(\omega)e^{j \omega L} \} + h[L] ) $$

Where the thing in the parenthesis is a real valued function and has zero phase term, therefore we call it $B(\omega)$ and deduce that the DTFT $H(e^{j\omega})$  of the causal, symmetric, real valued impulse response $h[n]$ is

$$ \boxed{ H(e^{j\omega}) = B(\omega) e^{-j\omega L} }$$

from which we find that the associated phase response of the filter is:
$$\phi_h(\omega) = -\omega L $$

and the associated group delay is therefore
$$\boxed{ \tau = - \frac{d \phi_h(\omega)} {d\omega} = L}$$

samples...

% This script demonstrates the use of time-aliasing in an advantageous way
% in obtaining the "central" portion of the circular convolution implemented
% via the DFT approach...
%
% Let x[n] of length Lx and h[n] of length Lh = 2*L+1, are convolved linearly
%
%                     y[n] = x[n] * h[n]
%
% so that "full" length of y[n] is:  Ly = Lx + Lh -1
%
% However in some cases we are only interested in those samples of length Lx
% that relate to input x[n] which are at the "central" portion of y[n], whereas
% the initial and final L samples carry transient information. Hence we take
%
%                    yc[n] = y[L:L+Lx-1]
%
% When a linear convolution is implemented with DFT, the length of DFT must
% satisfy the following; if we are to avoid time-aliasing: M >= Lz=Lx+Lh-1
% However this method computes ALL the full samples of y[n] in addition to the
% central ones. Hence if we wish to compute only the central portion of
% linear convolution y[n], then we can allow some aliasing in the DFT method
% The amount of which is indicated in the length of DFT as:
%
%                    M = Lz - L
%
% Then this will produce a result ya[n] whose last Lx samples carry the
% required information such that
%
%                    yc[n] = y[L:L+Lx-1] = ya[Lz-L-Lx:Lz-L-1]
%
% in matlab notation the indices will be added by 1 as:
%
%                    yc(:) = y(L+1:L+Lx) = ya(Lz-L-Lx+1:Lx-L);
%
% ==============================================================================
clc; clear all; close all;

% S0 - Define the parameters:
% ---------------------------
Lx = 32;
L = 15;
Lh = 2*L+1;
Lz = Lx+Lh-1;

% S1 - Generate the signals:
% --------------------------
x = randn(1,Lx);
h = randn(1,Lh);

% S2 - implement time-domain convolution first:
% ---------------------------------------------
yt = conv(x,h);
y1 = yt(L+1:L+Lx);           % get the central portion of the convolution

% S3 - Implement frequency domain DFT based aliased-circular conv:
% ----------------------------------------------------------------
M = Lz-L;       % DFT length (with aliasing)

yt = real( ifft( fft(x,M).*fft(h,M) , M) );     % compute aliased result
y2 = yt(L+1:end);

% S4 - compare the results
% ------------------------
figure,stem(y1,'b')
hold on
stem(y2,'g')

figure,stem(y1-y2);title('the error');


How to verify and plot expression for sampled signal is correct? Given the signal $x_a(t)=4cos(2\pi\times40t+\pi/3)+10cos(2\pi\times160t-\pi/6)$, I am to sample it at $f_s=200Hz$ and find the expression for $x[n]$
My process:
$$x[n]=x_a(nT)=x_a({n\over f_s})=4cos(0.4n\pi+\pi/3)+10cos(1.6n\pi-\pi/6)$$
and due to aliasing:
$$x[n]=4cos(0.4n\pi+\pi/3)+10cos(0.4n\pi+\pi/6)$$
How do I verify that the expression for the discrete-time signal is correct?
What I've tried so far:

A1 = 4; A2 = 10; f1 = 40; f2 = 160; th1 = pi/3; th2 = -pi/6; fs = 200;
dt = 0.0001;
t = 0:dt:0.2;
x_a = A1*cos(2*pi*f1.*t + th1) + A2*cos(2*pi*f2.*t + th2);
plot(t, x_a);
hold on;
n = 0:0.2*fs;
x_n = A1*cos(0.4*pi.*n + pi/3) + A2*cos(0.4*pi.*n + pi/6);
stem(n/fs,x_n,':', 'LineWidth', 1.5);
hold off;

And I get the following plot:

So it seems like the discrete signal's values match those of the continuous signal.
Is plotting both signals and checking if the sampled values match the right way of verifying?
I'm still wrapping my head around sampling so any additional information that can help me understand the process and the outcome better is more than welcome.
Also, it seems like there's no way to properly reconstruct the signal with the given sampling frequency so I'm assuming the point of this exercise is just to test my knowledge of sampling and this is not a practical way of sampling, correct?
Thank you!

The best way is to let the program to sample the function for you.
You can define a function. And you don't even need a separate program, you can use an anonymous function, @(t) expr for MATLAB/Octave (That seems to be the program you are using).
  x_a = @(t) A1*cos(2*pi*f1.*t + th1) + A2*cos(2*pi*f2.*t + th2);

Then you could plot with
  plot(t, x_a(t), '-', n/fs, x_a(n/fs), 'o');

since you are ploting x, x_a(x) there is nothing that could go wrong in the sampling, the time values you pass will be replaced automatically in your expression.
If you want the expression you can use symbolic computation, however when I need to do so I use wxmaxima or sympy


Detecting direction of sound using several microphones First of all, I've seen a similar thread, however it's a bit different to what I'm trying to achieve. I am constructing a robot which will follow the person who calls it. My idea is to use 3 or 4 microphones - i.e. in the following arrangement in order to determine the from which direction the robot was called:

Where S is source, A, B and C are microphones. The idea is to calculate phase correlation of signals recorded from pairs AB, AC, BC and based on that construct a vector that will point at the source using a kind of triangulation. The system does not even have to work in real time because it will be voice activated - signals from all the microphones will be recorded simultaneously, voice will be sampled from only one microphone and if it fits the voice signature, phase correlation will be computed from the last fraction of second in order to compute the direction. I am aware that this might not work too well i.e. when the robot is called from another room or when there are multiple reflections.

This is just an idea I had, but I have never attempted anything like this and I have several questions before I construct the actual hardware that will do the job:

Is this a typical way of doing this? (i.e. used in phones for noise cancellation?) What are other possible approaches?
Can phase correlation be calculated between 3 sources simultaneously somehow? (i.e. in order to speed up the computation)
Is 22khz sample rate and 12bit depth sufficient for this system? I am especially concerned about the bit depth.
Should the microphones be placed in separate tubes in order to improve separation?

To extend Müller's answer,

  Should the microphones be placed in separate tubes in order to improve separation?

No, you are trying to identify the direction of the source, adding tubes will only make the sound bounce inside the tube which is definitely not wanted. The best course of action would be to make them face straight up, this way they will all receive similar sound and the only thing that is unique about them are their physical placements which will directly affect the phase. A 6 kHz sine wave has a wavelength of $\frac{\text{speed of sound}}{\text{sound frequency}}=\frac{343\text{ m/s}}{6\text{ kHz}}=5.71\text{ mm}$. So if you want to uniquely identify the phases of sine waves up to 6 kHz, which are the typical frequencies for human talking, then you should space the microphones at most 5.71 mm apart. Here is one item that has a diameter that is less than 5.71 mm. Don't forget to add a low pass filter with a cut-off frequency at around 6-10 kHz.

Edit

I felt that this #2 question looked fun so I decided to try to solve it on my own.

  Can phase correlation be calculated between 3 sources simultaneously somehow? (i.e. in order to speed up the computation)

If you know your linear algebra, then you can imagine that you have placed the microphones in a triangle where each microphone is 4 mm away from each other making each interior angles $60°$.

So let's assume they are in this configuration:

       C
      / \
     /   \
    /     \
   /       \
  /         \
 A - - - - - B

I will...

use the nomenclature $\overline{AB}$ which is a vector pointing from $A$ to $B$
call $A$ my origin
write all numbers in mm
use 3D math but end up with a 2D direction
set the vertical position of the microphones to their actual wave form. So these equations are based on a sound wave that looks something like this.
Calculate the cross product of these microphones based on their position and waveform, then ignore the height information from this cross product and use arctan to come up with the actual direction of the source.
call $a$ the output of the microphone at position $A$, call $b$ the output of the microphone at position $B$, call $c$ the output of the microphone at position $C$

So the following things are true:

$A=(0,0,a)$
$B=(4,0,b)$
$C=(2,\sqrt{4^2-2^2}=2\sqrt{3},c)$

This gives us:

$\overline{AB} = (4,0,a-b)$
$\overline{AC} = (2,2\sqrt{3},a-c)$

And the cross product is simply $\overline{AB}×\overline{AC}$

$$
\begin{align}
\overline{AB}×\overline{AC}&=
\begin{pmatrix}
4\\
0\\
a-b\\
  \end{pmatrix}
×
  \begin{pmatrix}
2\\
2\sqrt{3}\\
a-c\\
  \end{pmatrix}\\\\
&=\begin{pmatrix}
0\cdot(a-c)-(a-b)\cdot2\sqrt{3}\\
(a-b)\cdot2-4\cdot(a-c)\\
4\cdot2\sqrt{3}-0\cdot2\\
  \end{pmatrix}\\\\
&=\begin{pmatrix}
2\sqrt{3}(b-a)\\
-2a-2b-4c\\
8\sqrt{3}\\
  \end{pmatrix}
\end{align}
$$

The Z information, $8\sqrt{3}$ is just junk, zero interest to us. As the input signals are changing, the cross vector will swing back and forth towards the source. So half of the time it will point straight to the source (ignoring reflections and other parasitics). And the other half of the time it will point 180 degrees away from the source.

What I'm talking about is the $\arctan(\frac{-2a-2b-4c}{2\sqrt{3}(b-a)})$ which can be simplified to $\arctan(\frac{a+b+2c}{\sqrt{3}(a-b)})$, and then turn the radians into degrees.

So what you end up with is the following equation:

$$\arctan\Biggl(\frac{a+b+2c}{\sqrt{3}(a-b)}\Biggr)\frac{180}{\pi}$$

But half the time the information is literally 100% wrong, so how.. should one.... make it right 100% of the time?

Well if $a$ is leading $b$, then the source can't be closer to B.

In other words, just make something simple like this:

source_direction=atan2(a+b+2c,\sqrt{3}*(a-b))*180/pi;
if(a>b){
   if(b>c){//a>b>c
     possible_center_direction=240; //A is closest, then B, last C
   }else if(a>c){//a>c>b
     possible_center_direction=180; //A is closest, then C last B
   }else{//c>a>b
     possible_center_direction=120; //C is closest, then A last B
   }
}else{
   if(c>b){//c>b>a
     possible_center_direction=60; //C is closest, then B, last A
   }else if(a>c){//b>a>c
     possible_center_direction=300; //B is closest, then A, last C
   }else{//b>c>a
     possible_center_direction=0; //B is closest, then C, last A
   }
}

//if the source is out of bounds, then rotate it by 180 degrees.
if((possible_center_direction+60)<source_direction){
  if(source_direction<(possible_center_direction-60)){
    source_direction=(source_direction+180)%360;
  }
}

And perhaps you only want to react if the sound source is coming from a specific vertical angle, if people talk above the microphones => 0 phase change => do nothing. People talk horizontally next to it => some phase change => react.

$$
\begin{align}
|P| &= \sqrt{P_x^2+P_y^2}\\
&= \sqrt{3(a-b)^2+(a+b+2c)^2}\\
\end{align}
$$

So you might want to set that threshold to something low, like 0.1 or 0.01. I'm not entirely sure, depends on the volume and frequency and parasitics, test it yourself.

Another reason for when to use the absolute value equation is for zero crossings, there might be a slight moment for when the direction will point in the wrong direction. Though it will only be for 1% of the time, if even that. So you might want to attach a first order LP filter to the direction.

true_true_direction = true_true_direction*0.9+source_direction*0.1;

And if you want to react to a specific volume, then just sum the 3 microphones together and compare that to some trigger value. The mean value of the microphones would be their sum divided by 3, but you don't need to divide by 3 if you increase the trigger value by a factor 3.

I'm having issues with marking the code as C/C#/C++ or JS or any other, so sadly the code will be black on white, against my wishes. Oh well, good luck on your venture. Sounds fun.

Also there is a 50/50 chance that the direction will be 180 away from the source 99% of the time. I'm a master at making such mistakes. A correction for this though would be to just invert the if statements for when 180 degrees should be added.


Parameters of Butterworth filter that approximates the amplitude response of a Gaussian filter? For a given narrowband Gaussian filter with a specific central frequency and filter width, I need corners of a bandpass Butterworth filter whose amplitude response is close enough to the Gaussian filter. How should I choose the order and corner frequencies of the Butterworth filter ? The reason is I need to calculate statistics of the initial part of a time-series filtered with a Gaussian filter. However, the Gaussian filter is zero-phase and always shifts energy back and contaminates the initial part of the time-series. With the Butterworth filter, I can have a one-pass filter so that the filtered time-series is causal and there is no contamination of the initial part of the time-series.

The solution is to implement a FIR filter with the coefficients using the Gaussian response desired. The center dominant tap will set the delay of the filter, which will of course not be zero phase, and the result will be the causal one-pass filter desired.

Also related, see my post at the link below on how you can implement this filter as a cascade of 2 tap unity gain filters (it doesn't get much simpler than that!).

Gaussian FIR filter with no multipliers?


How to find smoothed estimates of the derivative and second derivative of a signal? I have a signal sampled at $\Delta t$: $f_i(t_i=i\Delta t)$ where $i = 0,\ldots,n-1$.
I want to find the first and second derivative of the signal: $f'(t)$ and $f''(t)$.

My first thought was to estimate the derivatives by central differences:

\begin{align}
f'(t_{i})&=\frac{f(t_{i+1})-f(t_{i-1})}{2\Delta t}\\
f''(t_{i})&=\frac{f(t_{i+1})-2f(t_{i})+f(t_{i-1})}{(\Delta t)^2}
\end{align}

However the signal may have a lot of high frequency noise that may cause quick fluctuations in $f'$ and $f''$.

What would be the best way to find "smoothed" estimates of $f'$ and $f''$?

It probably depends more on your data. Just know, since differentiation is a linear operation, if you choose any linear filter to smooth f' and f'', it is equivalent to smoothing f using that same filter, then taking its derivatives.

Can you post some pictures or more information about the signal you want to differentiate? Probably what you're looking for is some sort of lowpass filter to smooth the signal. A couple really simple options include a single-pole recursive filter like $y(n) = a \cdot x(n) + (1-a) \cdot y(n-1)$, or a Hann filter, which is just convolving the signal with a Hann window. The Hann filter option is nice because it's linear-phase. If you know the frequency range you care about, you can just design a suitable lowpass filter in the frequency domain.


What is meant by "stochastic sampling"? What exactly is meant by "stochastic sampling" and is it profoundly different from the regular Nyquist-Shannon sampling theorem? Is it related to sampling a stochastic process?

Stochastic sampling doesn't have anything to do with sampling stochastic waveforms. It simply means that instead of sampling at regular intervals, the waveform is sampled randomly.

Recall that in a sampling scheme per the Nyquist-Shannon sampling theorem, a continuous signal $x(t)$ on $\mathbb{R}$ is sampled as $x[n]=x(nT),\ n\in\mathbb{Z}$, where $T$ is the sampling interval and $f_s=1/T$ is the sampling frequency. If the maximum frequency in the signal is $f_{max}$, then $f_s$ must be such that $f_s\geq 2f_{max}$ so as to avoid aliasing. For ease of comparison with stochastic sampling later on in the answer, let me redefine the sampling in a slightly different form than usual as

$$
\begin{align}
s(t)&=\sum_{n=0}^{f_s\tau -1}\delta(t-nT)\\
x[n]&=x(t)\cdot s(t)
\end{align}
$$
where $\delta(t)$ is the Dirac delta function and $x(t)$ is only sampled on the interval $[0,\tau]$.

If you actually think about it, regular sampling is pretty limiting in practice. Aliasing crops up in several places, and probably a well known and visible effect is the Moiré patterns which can be reproduced at home by taking a photo of regular patterns displayed on a television (examples below).

However, this is always a problem with cameras, but never with your eyes if you were to see the pattern directly! The reason is because the photoreceptors in your retina are not laid out in a regular pattern unlike the CCD in a camera. The idea behind (not necessarily the idea that led to its development) stochastic sampling is very similar to the non-regular layout of photoreceptors in the eye. It is an anti-aliasing technique which works by breaking up the regularity in the sampling.

In stochastic sampling, every point in the signal has a non-zero probability of being sampled (unlike regular sampling where certain sections will never be sampled). A simple uniform stochastic sampling scheme can be implemented over the same interval $[0,\tau]$ as

$$
\begin{align}
s(t)&=\sum_{n=0}^{f_s \tau -1}\delta(t-t_n),\quad t_n\sim \mathcal{U}(0,\tau)\\
x[n]&=x(t)\cdot s(t)
\end{align}
$$

where $\mathcal{U}(0,\tau)$ is the uniform distribution on the interval $[0,\tau]$.

By sampling stochastically, there is no "Nyquist frequency" to talk about, so aliasing will no longer be a problem as before. However, this comes at a price. What you gain in anti-aliasing, you lose by noise in the system. The stochastic sampling introduces high-frequency noise, although for several applications (especially in imaging), aliasing is a much stronger nuisance than noise (e.g., you can see the Moiré patterns easily in the above images, but to a lesser extent the speckle noise).

As far as I know, stochastic sampling schemes are almost always used in spatial sampling (in image processing, computer graphics, array processing, etc.) and sampling in the time domain is still predominantly regular (I'm not sure if people even bother with stochastic sampling in the time domain). There are several different stochastic sampling schemes such as Poisson sampling, jittered sampling, etc., which you can look up if you're interested. For a general, low key introduction to the topic, see

  M. A. Z. Dippé and E. H. Wold, "Antialiasing Through Stochastic Sampling", SIGGRAPH, Vol. 19, No. 5, pp. 69-78, 1985.


Filtering discrete data Is is possible to use a digital filter (eg. Butterworth) for smoothing sensor data that is sampled at a non-constant rate? The time between two consecutive samples is between 12 to 40 ms.

The most naive and easiest method would be to interpolate the data points. If you know what is the maximum frequency present in your signal then you can specify the sampling period. Let's say you pick the lowest possible: 12ms. Now you can perform some interpolation of your data. Normally I would suggest you to perform the sinc interpolation, but you can start with linear/spline/cubic.

In MATLAB/Octave you can use the interp1 function to interpolate the data points. If your data is called x, known time stamps are t, and new time stamps are t_new, then you can do the following to calculate new dataset X:

t_new = 0:0.012:max(t); %spaced by 12ms up to maximum value of vector t
X = interp1(t, x, t_new, 'spline');

This is very simple approach but should work - although you might get some aliasing. You can also check out this article for more intelligent methods: Nonuniform Sampling and Reconstruction in Shift-Invariant Spaces.

Now when you have your dataset, you can perform any filtering you like. Can I ask you what your sensor is so it's providing such data?


Why can Fourier series be used on a non-repetitive function? I was just wondering why Fourier series can be used on the function in the linked image. This is since I thought the function had to repeat itself to use Fourier series on it. Or is it saying a period is from 0-5s and so it will repeat every 0-5s? like the 2nd graph i made?

I can find the answer easily just not sure why fourier series can be used?

Fourier series can be calculated for any signal in close section (from $T_0$ to $T_1$). In fact, I was taught that this is how Fourier series is defined.

But Fourier series is made up of ($T_0 - T_1$) cyclic functions, therefore, if you go beyond the original section you will just get repetitions as you displayed.


FFT amplitude or magnitude Can I use the word amplitude instead of magnitude when I describe FFT bins? I dont see any similar word in my language.

Even i had confusion differentiating between these two terms at the beginning ,
have a look at this explanation from one of the Award winning DSP books.

Amplitudes,

Magnitudes,


why the number of filter coefficients in FIR filter has to be an odd number? I calculated the order of the FIR filter to be 31(N). so, the number of coefficients has to be 32 (N+1). so, I have to increase it to 33 to make it an odd number.

Why the number of filter coefficients is required to be an odd number?

There is no such requirement. However, for linear-phase (symmetric) FIR filters, the delay can be expressed as:

$$
D = \frac{order}{2} = \frac{N_{taps} - 1}{2} \text{ samples}
$$

If the number of taps is odd, then the delay of the filter is an integer number of samples, which is desirable for some applications. If you have an even number of taps, then you end up with a half-sample delay.


Expectation maximization of moving average with binary source input I am trying to do blind system identification of a univariate linear FIR model:

I am unsure if the approach is correct or not and any help to further proceed with the maximization will be great. Thank you.

I think your E-step is correct(only one term missing in the last expression: $-N\ln\sigma_w$ ).
To obtain the M-step you have to differentiate with respect to all your parameters. You don't have to include $u(n)$  in $\theta$, since it's defined by $p$. So you compute $\frac{\partial Q}{\partial\theta}=0$. Solving for $\theta$, and this is your new $\theta$ i.e. $\theta_{m+1}$. E.g. for $h$ we gets $\frac{\partial Q}{\partial h}=-(-2\sum y_nz_n^s+2h\sum P_n^s)/(2\sigma_w^2)=0$.
Solving: $h=h_{m+1}=\sum y_nz_n^s(\sum P_n^s)^{-1}$.

-For $p$:
$Q'(p)=\sum z_n^s/p -(N-\sum z_n^s)/(1-p)=0,p=\sum z_n^s/N$.

This is typical Bernoulli/binomial. So it seems right.

-For $\sigma_w$:
$Q'(\sigma_w)=-(\sum y_n^2-2h^T\sum y_nz_n^s+h^Th\sum P_n^s)\cdot (-2/(2\sigma_w^3))-N/\sigma_w=0,
\sigma_w^2=(\sum y_n^2-2h^T\sum y_nz_n^s+h^Th\sum P_n^s)/N$

and updating the parameters in a loop until convergence.


How to conclude LTI, causality and BIBO stability of a system represented by a differential equation? I have started to learn about systems represented by differential equations in Oppenheim's Signals & Systems, and I got really confused about it. I am trying to understand how I can show that a system is LTI (or not), causal (or not) and BIBO stable (or not) for any input.

I know the following (please correct me if I'm wrong):

If the differential equation has zero auxiliary conditions then the system is linear. Therefore, showing linearity is easy.
A system is defined by the response to the Dirac input.
Causal, LTI $\iff$ initial rest (mathematically meaning if $x(t)=0$ for $t<0$ then $y(t)=0$ for $t<0$).

I am not sure if I got it right but if I'm willing to use #3 in order to show that the system is causal and LTI, then it has to have initial rest for any input. How can something like that be shown if as I say it has to be for any $x(t)$? I thought that maybe using #2 will be helpful (solving for the differential equation for delta function and then somehow to show), can it be?

In conclusion, I would thankful if you explain to me how show those features when I got a system that represented by differential equation. For example, in case that we got that simple ODE:

$$\frac{d^2y}{dt^2}-3\frac{dy}{dt}+2y(t)=x(t) \ \ \ ,y(0)=\frac{dy}{dy}=0$$

I take your equation:
$\frac{d^2y}{dt^2}-3\frac{dy}{dt}+2y(t)=x(t)$

Laplace transform will be:
$ s^2Y(S) - 3sY(S) +2Y(s) = X(s) $

Now I can find transfer function: $H(s) \triangleq \frac{Y(s)}{X(s)}$

$H(s)=\frac{1}{s^2-3s+2}$

Poles are the solution of denominator = 0 and zeros are the solution of numerator = 0. In this case there are not zeros because numerator is always equal to 1, denominator = 0 is a second order equation that can be easily solved.

$ Poles[H(s)] = \{s=1 ; s=2\} \\ Zeros[H(s)] = \{\} $

There are two real poles.

To determine if the system is BIBO stable you can see if all poles have real part $ < 0 $. In this is example there is at least one pole that has real part  $ \geq 0 $ (both in reality) so this system is UNSTALBE.

To be causal the number of zeros must be $\leq$ of the number of poles in this case we have $ 0 \leq 2 $ that it's true and means that your system is CAUSAL.

Time invariant system is a system that doesn't depends on absolute time but only on relative time from time origin. In a transfer function generally there are not information to determine if it's a time invariant system or not beacuse it should be an external specification of your problem. Btw in DSP you consider systems time invariant almost ever.

Finally this is only a quick overview of zero-pole analysis you should study the argument as well and the prove of the concepts I  used.


How to estimate the period of a pulse train? I am trying to find the period of a pulse train like signal.

I have reviewed the following questions but I did not find a clear-cut answer to the problem, 12892, 16502, 6260.

I do not know the period nor the pulse width. I work with windows that will contain between 2 and 16 pulses. Pulses may vary in length but will usually be shorter than half the period. The signal may be noisy, although I would be happy with a solution for the case of the (relatively) clean signal. I can work out a minimal pulse width if it helps.

Here are plots of the type of signal I am working with.

I am primarily interested in the period of the signal. The analysis needs not to be online.

For the first example, I'm looking to find a result of about 136.5.

My first attempt was to look at the signal in the frequency domain. But a pulse train in time domain is also a pulse train in frequency domain. I got this plot for the clean signal.

Next I tried autocorrelation. I got the following:

It is much better. I need to get the value automatically though, and I'm not too sure how to find that second local maxima without introducing too much ad-hoc thresholds that may fail when there are more or less peaks in the window. Additionally, I'm looking for a floating point result.
Also, the autocorrelation for the noisy signal was much less clear cut, even after smoothing the input.

There is probably another approach by detecting the rising edges in time domain?

Question:

Is there a standard way to estimate the period of a pulse train?
If autocorrelation is the right approach, how to automatically find the floating point coordinate of the peak from the autocorrelation result without too many assumptions?

If all you need to estimate is the period of your pulse train then you can threshold the original waveform and work on the clean one by detecting the rising and falling edges.

An FFT based approach would return an average estimate of its "frame of observation" anyway whose accuracy would depend on the length of the frame. (that is, an average estimate of the buffer that is passed at its input)

If you are concerned about the noise, use a short moving average filter (for instance a length 3 or length 5) to smooth the curve prior to estimating and applying the threshold. If absolute timings are essential for your application don't forget to take into account the group delay of the filter.


how to extract telephone audio only from an mp3 I have mp3s of radio talk shows in which the host interviews people on the phone. The host uses a microphone, and the interviewee uses either a landline phone or a cell phone.

I want to save just the telephone audio to a different file. Is that feasible?

There might be much discrepency in the signal's characteristics of the interviewee that depends partialy on the phones they're using, their mobile operator, etc...
It is much easier to extract the radio's speaker signal. Then, everything that does not come form the radio's host must be from one of the interviewees.
What you need to do is to analyse the signal of the radio's speaker that might contains more prominent low-pitched components for example (depends obviously on the station you recorded from). Also use a window of analysis sufficiently large to take into account that every stream must last at least 500/1000ms. That might introduce a little stability in your extraction.
If you want to split the different interviewees, you can pose the simple condition that interviewee/host must alternate.


MP3 Decoding Discrepancy Is it possible for two different decoders to output two different sample results for the same MP3? I'm decoding an mp3 on Matlab and it produces two columns (presumably one for each channel). When I extract the sample data from the same file using AVAssetReader (iOS), I get different results. When I tried the same code with a WAV file, the two results matched, so I'm guessing the code is fine (although the Matlab result had a few extra blank frames).

If I, non-scientifically, choose random values from my AVAssetReader MP3 result, and search for them in the Matlab results, they tend to exist, but not always. I don't know if this is just coincidence, different rounding off errors or if they genuinely match. The problem is, when they match, the surrounding numbers differ from one another. For example, 0.0021 occurs in both results, but the surrounding numbers on Matlab are -0.0649,0.0147 etc.. (which don't appear near the 0.0021 in the AVAssetReader results).

What could cause this discrepancy, does it cause a serious lack of integrity in my app and how would I fix it? Thank you.

It is possible that two decoders will have different results and be compliant with standard. The standard "ISO/IEC 11172-4 (part 4) compliance testing" specify two classes of decoders:
"accurate" and with "limited accuracy" *).

The accuracy is measured by comparing decoded stream with reference stream (reference input and output is specified in standard). If the RMS of difference is below $2^{-15}/\sqrt{12}$ the decoder is classified as accurate and  if below $2^{-11}/\sqrt{12}$ it is classified as "limited accuracy" one.

So if you are going to use different decoders or even the same decoder on different processor architectures you have to be prepared for this kind of discrepancy in output. Do not know your application or intend but maybe you can mimic the standard and work with RMS of error too?

*) I do not have text of the standard in English in front of me so actual name used in standard may be different.


What is z equal to in z-transform? In some places, it is said that z is equal to:
$$z = e^s \quad where \quad s = \sigma + j \Omega $$
But in some places, it is said that z is equal to:
$$z = e^{sT_s} \quad  $$
where Ts is a sampling interval.
Why are there two different definitions?

When discussing the z-transform, it is important to note that there are two different frequencies that are innately involved. In many text books, the two frequencies are denoted by $\omega$ (discrete-time) and $\Omega$ (continuous-time). We'll use this convention.
The z-transform uses the "discrete frequency" $\omega$ which is derived from the "analog frequency" $\Omega$ via sampling at a rate $T_s$ as
$$\omega = {\Omega}T_s$$
If we're on the unit circle, which your question implies as the starting condition, then the $\sigma$ term in $s = \sigma + j\Omega$ goes to zero and so
$$s = \sigma + j\Omega \Rightarrow s = j\Omega$$
We can now relate the z-transform as the discrete-time equivalent of the Laplace transform and so
$$z = e^{sT_s} = e^{j{\Omega}T_s} = e^{j\omega}$$
Which is exactly the substitution made between the Fourier transform of a discrete sequence $x[n]$ and it's z-transform:
$$X(e^{j\omega}) = \sum_{-\infty}^{\infty}x[n]e^{-j\omega n}$$
$$X(z) = \sum_{-\infty}^{\infty}x[n]z^{-n}$$
I have never seen $z$ be defined by $z = e^{s}$, which implies the use of the continuous-time frequency and is not how the z-transform is defined.
So the confusion might come down to abuse of notation between how $\omega$ and $\Omega$ are defined in the sources you are referencing.


Is it better to use sampling rates which are a power of 2? 32-64-128-256-512Hz? instead of using for example 30- 60-100-250-500Hz? I am recently working with a sensor that I its sampling frequency can only be configured to either 32hz- 64hz- 128hz-256hz... and I was wondering whether using such sampling frequencies have any benefits. For example is it better to use 128 Hz than use 100 Hz? or 64Hz rather than 60 Hz?

There's a lot of things that influence your choice of sampling rate.

First of all, of course, the Nyquist theorem, which says you need to sample at more than twice of the signal's single-sided bandwidth $B$ (for real signals)¹.

Then, you often choose a sampling rate that is useful for your system. For example, a lot of 3G/4G/5G systems use sampling rates of 30.72 MHz (or multiples thereof). That's not an elegant multiple of anything (as far as I know) – it just happens to be what the symbol rates in the standard suggests, so it's very useful if you're already sampling at rate that is rational relative to what you'll be using in the end.

That's a generally important rule: Pick a sampling rate that makes sense for what you'll be doing to the signal afterwards.

Then, especially in audio and lower-frequency sampling, you choose sampling rates that are multiples of other rates, so that it's easy to keep multiple subsystems in synchronization. For example, your CD audio sampling rate is 44100 Hz – first of all, that's an integer multiple of 100 Hz, and that's good, because it means it's easy to keep not only the second counter on your CD player correct by simply counting the samples and dividing by 44100, but also things like studio hardware, which tends to have higher "human-readable" time resolution. But, why the strange rate 44100? Well, that rate happens to be $44100=2\cdot2\cdot3\cdot3\cdot5\cdot5\cdot7\cdot7$. Back in the 1980s, when that rate came up, it was harder than it is now to derive good clocks from each other, unless they had an integer factor between them. So, a clock with "all" the small integer factors was chosen. A curious choice!

Another reason to pick a specific sampling rate is noise alignment, or intentional non-alignment. While it's not a sampling rate itself, the NTSC frame rate of 29.97 frames a second is still something that you find in video systems today, so you find a lot of video capturing ADCs that run at a high multiple of 29.97 Hz. And that frequency was chosen so that the zeros of the color information shaping filter didn't interfere with the brightness info.

In the very same motivation, you'll find noise-sensitive systems where the ADC rate is as far as possible from being rationally related to the other dominant digital clocks in the system – so that the noise from these clocks distributes more evenly to the digitized signal.

In other systems, you painstakingly try to align all these clocks – so that clock noise ends up being narrow interferers in the digital signal, and thus can be digitally filtered out using notch filters.

Of course, you also need to be able to get any rate you'd want to use – you can't just by a 876.543 kHz oscillator, even if you wanted to, these don't exist commercially. Sure, you can synthesize one with either a lot of jitter or a lot of effort, but in the end, complexity and cost restraints are relevant for almost any system design.

As you can see, a lot of conflicting reasons!

¹ that's not necessarily twice the maximum frequency in the signal! A lot of systems, especially superheterodyne systems, i.e. systems that have a strictly bandslimited signal on an intermediate frequency (IF) use subsampling, i.e. the sample rate is $> 2\cdot B$, but $\ll 2\cdot f_\text{max}$. The aliasing that occurs doesn't distort the information content, but aliases the signal down from the IF to the Nyquist band. It's a very useful, very common technique.


Filter design for an unusual EEG experiment I'm designing a quite unusual (in terms of signal analysis) EEG experiment, which forces me to design my own filter. I have never done this before, so kindly please check my work and suggest corrections.
The analysis consists of cutting a narrow frequency band (2Hz) from 6-second EEG recordings(4*) (sampling frequency 256 Hz). I expect that after such filtering, I will receive a "somehow oscillated" signal (not in a strict mathematical sense, just with notable local maxima (peaks) and local minima (troughs))(5*). Now I'm interested just at the very end of this 6-sec window, the last dozen or so data points. I want to determine (by the simple algorithm provided by myself) if the signal ends "near the peak" or "near the trough" or in either of the above. So, all 6-sec windows are compressed to one of the three categories. You can think of this as a very rough estimate of the phase at the end of the window.
The tricky part is that I need to analyze overlapping frequency bands ( 4-6 Hz, 4.25 - 6.25 Hz, 4.5 - 6.5 Hz, etc.), so it will be great to have an extremely steep roll-off even at the price of amplitude disturbance (most of this information will be lost anyway) - that's why I came up with an elliptic IIR filter. Zero-phase shifting is also crucial so, forward-backward (filtfilt) filtering is applied. However, the devil is in the details (filter order, maximum ripple, minimum attenuation), and I must honestly admit that I determined these parameters by trial and error.

from scipy import signal
from mne.viz import plot_filter

sfreq = 256.
f_s = 6.
f_p = 8.
flim = (f_s - 1.0, f_p + 1.0)  # limits for plotting

nyq = sfreq / 2.  # the Nyquist frequency is half our sample rate
freq = [0, f_s, f_s, f_p, f_p, nyq]
gain = [0, 0, 1, 1, 0, 0]

ftype = "ellip"
order = 8
sos = signal.iirfilter(order, [f_s / nyq, f_p/nyq],
                       btype='bandpass', rp=5, rs=35, ftype=ftype, output='sos')
plot_filter(dict(sos=sos), sfreq, freq, gain,
            f'{ftype} order={order}', flim=flim, compensate=True)

Filters Visualisations:

So my questions are as follow:

Is my filter "reasonable"?
Can I do better?
Can you recommend to me some learning materials for a filter design.
(*) This is obviously cut from continuous EEG signal - can it be just cut (boxcar window), or should I apply something different (eg. Hamming window)?
(*) Is my expectation justified?

Is my filter "reasonable"?

I would say no. The time domain ringing time of your filter far exceeds the length of your signal, there isn't enough signal for the output to stabilize.

Can I do better?

Difficult to say. Filter design involves complicated tradeoffs between many different factors and the "best" really depends on the specific needs & sensitivities of your application. You have partially explained what you want to do but not in enough details to give a recommendation. I'm skeptical about "looking at the phase" approach. Signal phase has only meaning if you have a well defined time reference for $t=0$, which most signals don't.

Can you recommend to me some learning materials for a filter design.

I personally like that one https://ccrma.stanford.edu/~jos/filters/ but I'm also mostly and audio guy.

(*) This is obviously cut from continuous EEG signal - can it be just cut (boxcar window),
or should I apply something different (eg. Hamming window)?

Depends a bit on your application details. Ideally your filter "ringing" is significantly shorter than your signal. Then you can just take a "large" chunk of signal, filter it and cut off the edges to remove the filter transients.

(*) Is my expectation justified?

Can't say without looking at your signal first. A sufficiently narrow band signal will look like an amplitude and/or phase modulated sine wave. However, if your original signal is very noise-like, you'd have to make it VERY narrow to look like that.


Ambiguity in the term 'dimension'? We used to classify signals as 1D and 2D etc ie one dimensional and two dimensional. For example a periodic square wave signal is 1D and an image is a 2D signal etc (reference - Signals and systems by Simon Haykin and Barry Van Veen, 2nd edition , page 2).

But the same periodic square wave signal can be decomposed using fourier series to infinite sinusoids with different frequencies. In the linear algebra terms these infinite orthogonal sinusoids forms the basis and the the dimension of a periodic square wave is infinite.

So actually which is the actual dimension or what is dimension?

Dimension in mathematics is defined as the number of independent components of a structure (with some generalizations, like fractal dimensions). Now which number of dimensions you assign to e.g a discrete image (as in picture) depends on what aspect of the image you are describing. As a vector space where each pixel is a component, you get as many dimensions as pixels. The same for a discrete signal, the number of samples is your dimension of the vector space.

However, sometimes you're rather interested in the manifold the object is a function on. The real line is a 1-dimensional real manifold, and a real function on it would be a signal. So you can call the signal "one dimensional", explicitly referring to this property. Sometimes your function may map to the complex numbers or to R^2 (still from the real line!) instead, and you might want to call that a "two dimensional signal", and that's also correct. However, a two dimensional signal may also just be a function on a two dimensional manifold, like a plane.

So there's a lot ambiguity in verbal descriptions, which is why to be clear you should use exact mathematical notation. Opposite to common belief it was not just invented to torture students.


How to Simulate AWGN (Additive White Gaussian Noise) in Communication Systems for Specific Bandwidth I am trying to generate a AWGN waveform to add it to the signal of my simulated communication system. The operating bandwidth of the communication system is about B=3GHz and suppose that T=300K (my signal is comprised of very short pulses, I am simulating an IR-UWB system).

What I have done so far is:

N = k*T*B
sigma = sqrt(N./2)
noise = sigma.*randn(1,s_length)

where s_length is the number of samples for my useful signal and noise the awgn waveform.

Is this procedure correct? I have browsed through different questions but was unable to clarify it.

To create Band Limited AWGN all you need is randn in MATLAB.

The question lies only in how to set its Standard Deviation.
To illustrate that, let's assume our AWGN generator has PSD which equals to $ {N}_{0} $.
Namely we have AWGN with zero mean and variance equals to $ \delta(0) {N}_{0} $.

Assuming we have limited bandwidth channel, hence an ideal LPF is applied.
Assume its cutoff frequency is $ {F}_{LPF} = \frac{{F}_{S}}{2} $.

Hence the Variance of the band limited AWGN is (Multiplication in frequency) the integral over its PSD multiplied by the Norm of the LPF (The factor 2 is for integration over the range $ -{F}_{LPF}:{F}_{LPF} $):

$$ {Var}_{BandLimitedAWGN} = 2 {F}_{LPF} {N}_{0} = {F}_{S} {N}_{0} $$

Now, generate in MATLAB, using randn noise with the corresponding STD (By the data of your simulation).

Good Luck!

Some Remarks

The Variance of the noise is independent of the signal (At least in the classic model).
The Variance of the noise is a function only of the analog channel and the Analog to Digital converter.
The classic model assumes that if the signal is sampled at $ {F}_{s} $ an ideal LPF with a cut off frequency of $ {F}_{LPF} = \frac{{F}_{S}}{2} $.
As I said, this is the frequency which sets the variance of the Band Limited noise.
The input signal (Which is the transmitted signal + noise) may have any bandwidth it might have, after the LPF its bandwidth is limited.
In order to minimize the energy of the noise in the system the LPF band width and the sampling rate should be as low as possible (Namely, the bandwidth of the signal in interest).
Though if the next step is "Matched Filter" the SNR will be maximized for any finite energy white noise (Or colored if the "Colorization" is known and the Matched Filter is accordingly updated).


Circular Convolution Formula Deduction from DFT I am trying to get to the circular convolution formula from the product of the DFT of two input signals. I came across this webpage, which contains the following.

Could someone help me to solve that equation?

I have something like this

$$
x_3(n) = \frac{1}{N} \sum^{N-1}_{k=0} \exp{\left [ \frac{j 2 \pi k n}{N} \right ]}
\left [
\sum^{N-1}_{m=0} x_1(m) \exp{\left [ \frac{-j 2 \pi k m}{N} \right ]}
\right ]
\left [
\sum^{N-1}_{m=0} x_2(m) \exp{\left [ \frac{-j 2 \pi k m}{N} \right ]}
\right ]
$$

I am not sure what is next. The following doesn't seem to be right.

$$
\begin{aligned}
x_3(n) &= \frac{1}{N} \sum^{N-1}_{k=0}
\left [
\sum^{N-1}_{m=0} x_1(m) \exp{\left [ \frac{-j 2 \pi k m}{N} \right ] \exp{\left [ \frac{j 2 \pi k n}{N} \right ]}}
\right ]
\left [
\sum^{N-1}_{m=0} x_2(m) \exp{\left [ \frac{-j 2 \pi k m}{N} \right ]}
\right ]\\
\\
&= \frac{1}{N} \sum^{N-1}_{k=0}
\left [
\sum^{N-1}_{m=0} x_1(m) \exp{\left [ \frac{-j 2 \pi k (m-n)}{N} \right ]}
\right ]
\left [
\sum^{N-1}_{m=0} x_2(m) \exp{\left [ \frac{-j 2 \pi k m}{N} \right ]}
\right ]\\
\end{aligned}
$$

Let $x$ and $y$ be signals of $N$ samples each, numbered as $x(0),\ldots,x(N-1)$. Then their DFTs are $X$ and $Y$, which also have $N$ entries each:
\begin{eqnarray}
X(k) &=& \sum_{n=0}^{N-1}x(n)e^{-2\pi i kn/N},\\
Y(k) &=& \sum_{m=0}^{N-1}y(m)e^{-2\pi i k m/N},
\end{eqnarray}
where the indices run from $0$ to $N-1$.

The $k^{\textrm{th}}$ entry of the entry-by-entry product of $X$ and $Y$ is
\begin{equation}
\begin{split}
X(k)Y(k) ~=& \left(\sum_{n=0}^{N-1}x(n)e^{-2\pi i kn/N}\right)\left(\sum_{m=0}^{N-1}y(m)e^{-2\pi i k m/N}\right)\\
~=& \sum_{n=0}^{N-1}\sum_{m=0}^{N-1}x(n)y(m)e^{-2\pi i k(n+m)/N}
\end{split}
\end{equation}

Now we consider the $\ell^{\textrm{th}}$ entry of the IDFT of this entry-by-entry product:
\begin{equation}
\begin{split}
\mathsf{IDFT}(XY)(\ell) ~=& \frac{1}{N}\sum_{k=0}^{N-1}X(k)Y(k)e^{2\pi i\ell k/N }\\
~=& \frac{1}{N}\sum_{k=0}^{N-1}\left[\sum_{n=0}^{N-1}\sum_{m=0}^{N-1}x(n)y(m)e^{-2\pi i k(n+m)/N}\right]e^{2\pi i\ell k/N }\\
~=& \frac{1}{N}\sum_{n=0}^{N-1}\sum_{m=0}^{N-1}x(n)y(m)\sum_{k=0}^{N-1}e^{-2\pi i k(n+m-\ell)/N}
\end{split}
\end{equation}
The sum over $k$ is equal to 0 unless
\begin{equation}
n+m-\ell=0~\textrm{mod}~N~~~(m=\ell-n~\textrm{mod}~N),
\end{equation}
in which case it is equal to $N$:
\begin{equation}
\sum_{k=0}^{N-1}e^{-2\pi i k(n+m-\ell)/N} = N\delta_{m,\ell-n~\textrm{mod}~N},
\end{equation}
where $\delta_{a,b}$ is the Kronecker delta. One intuitive way to see this is to see that

$\bullet$ if the exponent is not 0, we are adding all $N$ of the $N^{\textrm{th}}$ roots of unity (view them in the complex plane), which will cancel one another out when added; As noted by the OP, my struck-through claim is true for all $n$, $m$, and $\ell$ only if $N$ is prime. The correct proof is to note that if $e^{-2\pi i(n+m-\ell)/N}\neq 1$, then
\begin{equation}
\begin{split}
\sum_{k=0}^{N-1}e^{-2\pi i k(n+m-\ell)/N} &=~ \sum_{k=0}^{N-1}\left(e^{-2\pi i (n+m-\ell)/N}\right)^k\\
&=~\frac{ 1- \left(e^{-2\pi i (n+m-\ell)/N}\right)^N}{1 - e^{-2\pi i (n+m-\ell)/N}}\\
&=~\frac{ 1- e^{N\times(-2\pi i (n+m-\ell)/N)}}{1 - e^{-2\pi i (n+m-\ell)/N}}\\
&=~\frac{ 1- e^{-2\pi i (n+m-\ell)}}{1 - e^{-2\pi i (n+m-\ell)/N}}\\
&=~\frac{ 1- 1}{1 - e^{-2\pi i (n+m-\ell)/N}} ~~=~~ 0.
\end{split}
\end{equation}

$\bullet$ if the exponent is 0, then we are adding 1 $N$ times.

(Another not as intuitive away can be found in this video.)

So far, we have
\begin{equation}
\begin{split}
\mathsf{IDFT}(XY)(\ell) ~=& \frac{1}{N}\sum_{n=0}^{N-1}\sum_{m=0}^{N-1}x(n)y(m)N\delta_{m,\ell-n~\textrm{mod}~N}\\
~=& \sum_{n=0}^{N-1}\sum_{m=0}^{N-1}x(n)y(m)\delta_{m,\ell-n~\textrm{mod}~N}
\end{split}
\end{equation}
When we perform the sum over $m$, the only nonzero term is the one for which $m=\ell-n~\textrm{mod}~N$, so
\begin{equation}
\begin{split}
\mathsf{IDFT}(XY)(\ell) ~=& \sum_{n=0}^{N-1}x(n)y(\ell -n~\textrm{mod}~N).
\end{split}
\end{equation}
The expression on the right-hand side is the $\ell^{\textrm{th}}$ entry of the circular convolution of $x$ and $y$.
\begin{equation}
\mathsf{IDFT}(XY)(\ell) = (x\circledast y)(\ell),
\end{equation}
so $\mathsf{IDFT}(XY) = x\circledast y$.


Perfect reconstruction: what properties must the analysis LPF satisfy I have been reading on conjugate quadrature filters (CQF). The linked article explains that, if we have a lowpassfilter $h_0$, we can create a set of filters for a perfect reconstruction filterbank. The analysis highpass filter $h_1$ and synthesis low pass/highpass filters $f_0$ and $f_1$ can be generated as follows:

  \begin{align}
h_1(n)&=-\left(-1\right)^n h_0\left(L-1-n\right)\\
f_0(n)&=h_0\left(L-1-n\right)\\
f_1(n)&=-\left(-1\right)^n h_0(n)=-h_1\left(L-1-n\right).
\end{align}

This works excellent for a set of Smith-Barnwell coefficients. Yet it does not work for every lowpass filter that could be designed. For instance a Kaiser-Bessel windowed sinc-function with a cutoff at $\pi/2$ (as described in the article as well) does not lead to perfect reconstruction.

Are there any particular requirements on $h_0$/$H_0$ before this will work ?

So the central question here is what the restrictions for the filters of this class of filterbanks (i.e. $N=2$-filter bank with the HPF generated from a lowpass/highpass transform of the LPF) is.

Let us for that first introduce a couple of helping definitions:

$$\begin{align}
\mathbf{h}(z) &= \begin{bmatrix}h_0(z)\\h_1(z)\end{bmatrix}&\text{analysis filterbank matrix}\\
\mathbf{r}(z) &= \begin{bmatrix}r_0(z)\\r_1(z)\end{bmatrix}&\text{reconstruction filterbank matrix}\\
\end{align}$$

With these in place: Let's consider the polyphase TypeI representation $A(z)$ (analysis) and TypeII representation $R(z)$ (reconstruction) of these filter banks (by the way, this is why we have these different types of polyphase deconstructions); I hope my notation is clear, $h_{1,0}$ is the zeroth polyphase component of the high pass filter.

$$\begin{align}
\mathbf{h}(z) &=
   \begin{bmatrix}
     h_{0,0}^{\mathrm I}(z^2) & h_{0,1}^{\mathrm I}(z^2)\\
     h_{1,0}^{\mathrm I}(z^2) & h_{1,1}^{\mathrm I}(z^2)
   \end{bmatrix}
   \begin{bmatrix}
     z^{-0}\\
     z^{-1}
   \end{bmatrix}\\
&= \mathbf A(z^2) \mathbf e_z(z)\\
\mathbf{r}^T(z) &=
   \begin{bmatrix}
     z^{-1} &
     z^{-0}
   \end{bmatrix}
   \begin{bmatrix}
     r_{0,0}^{\mathrm {II}}(z^2) & r_{0,1}^{\mathrm {II}}(z^2)\\
     r_{1,0}^{\mathrm {II}}(z^2) & r_{1,1}^{\mathrm {II}}(z^2)
   \end{bmatrix}\\
&= \mathbf{ \tilde e}_z(z) \mathbf R(z^2) \\
\end{align}$$

Now, looking at this, we simply define perfect reconstruction "intuitively" as the case where applying $\mathbf R$ to the output of $\mathbf A$ leads to the input of $\mathbf A$, of course delayed, and maybe scaled with a constant $c$:

$$\begin{align}
  h,\,r &\text{ form a perfect reconstruction fb}\\
  &\iff\\
  \mathbf R(z)\mathbf A(z) &= c z^{-k} \mathbf I_N & c \text{ const.},\, k\in\mathbb N
\end{align}$$

That gives us a clear method to move forward:

$$\begin{align}
  \mathbf R(z)\mathbf A(z) &=    \begin{bmatrix}
     h_{0,0}^{\mathrm {II}}(z) & h_{0,1}^{\mathrm {II}}(z)\\
     h_{1,0}^{\mathrm {II}}(z) & h_{1,1}^{\mathrm {II}}(z)
   \end{bmatrix}
   \begin{bmatrix}
     r_{0,0}^{\mathrm I}(z) & r_{0,1}^{\mathrm I}(z)\\
     r_{1,0}^{\mathrm I}(z) & r_{1,1}^{\mathrm I}(z)
   \end{bmatrix}\\
&= \begin{bmatrix}
h_{0,0}^{\mathrm {II}}(z)r_{0,0}^{\mathrm I}(z) +  h_{0,1}^{\mathrm {II}}(z)r_{1,0}^{\mathrm I}(z) &
h_{0,0}^{\mathrm {II}}(z)r_{0,1}^{\mathrm I}(z) +
h_{0,1}^{\mathrm {II}}(z)r_{1,1}^{\mathrm I}(z)\\
h_{1,0}^{\mathrm {II}}(z)r_{0,0}^{\mathrm I}(z) +
h_{1,1}^{\mathrm {II}}(z)r_{1,0}^{\mathrm I}(z) &
h_{1,0}^{\mathrm {II}}(z)r_{0,1}^{\mathrm I}(z) +
h_{1,1}^{\mathrm {II}}(z)r_{1,1}^{\mathrm I}(z)
\end{bmatrix}\\
&\overset!= \begin{bmatrix}c z^{-k}&0\\0&c z^{-k}\end{bmatrix}
\end{align}$$

EDIT now, let's do our jobs of inserting our conditions:

$$r_0  = h_0 (L-1-n)$$

notice that the substitution of $n$ by $L-1-n$ implies that the typeI polyphase decomposition of $r_0$ is the same as the typeII decomposition of $h_0$!

Hence, the diagonal entries of above matrix get a lot easier to calculate:

$$\begin{align}
  \mathbf R(z)\mathbf A(z) &=\begin{bmatrix}
h_{0,0}^{\mathrm {II}}(z)r_{0,0}^{\mathrm I}(z) +  h_{0,1}^{\mathrm {II}}(z)r_{1,0}^{\mathrm I}(z) &
h_{0,0}^{\mathrm {II}}(z)r_{0,1}^{\mathrm I}(z) +
h_{0,1}^{\mathrm {II}}(z)r_{1,1}^{\mathrm I}(z)\\
h_{1,0}^{\mathrm {II}}(z)r_{0,0}^{\mathrm I}(z) +
h_{1,1}^{\mathrm {II}}(z)r_{1,0}^{\mathrm I}(z) &
h_{1,0}^{\mathrm {II}}(z)r_{0,1}^{\mathrm I}(z) +
h_{1,1}^{\mathrm {II}}(z)r_{1,1}^{\mathrm I}(z)
\end{bmatrix}\\
&=\begin{bmatrix}
h_{0,0}^{\mathrm {II}}(z)h_{0,0}^{\mathrm {II}}(z) +
h_{0,1}^{\mathrm {II}}(z)-h_{1,0}^{\mathrm{II}}(z) &
h_{0,0}^{\mathrm {II}}(z)h_{0,1}^{\mathrm {II}} +
h_{0,1}^{\mathrm {II}}(z)-h_{1,1}^{\mathrm{II}}(z)\\
h_{1,0}^{\mathrm {II}}(z)h_{0,0}^{\mathrm{II}}(z) +
h_{1,1}^{\mathrm {II}}(z)-h_{1,0}^{\mathrm{II}}(z) &
h_{1,0}^{\mathrm {II}}(z)h_{0,1}^{\mathrm{II}}(z) +
h_{1,1}^{\mathrm {II}}(z)-h_{1,1}^{\mathrm{II}}(z)
\end{bmatrix}\\
&=\begin{bmatrix}
h_{0,0}^{\mathrm {II}}(z)h_{0,0}^{\mathrm {II}}(z) -
h_{0,1}^{\mathrm {II}}(z)h_{1,0}^{\mathrm{II}}(z) &
h_{0,0}^{\mathrm {II}}(z)h_{0,1}^{\mathrm {II}} -
h_{0,1}^{\mathrm {II}}(z)h_{1,1}^{\mathrm{II}}(z)\\
h_{1,0}^{\mathrm {II}}(z)h_{0,0}^{\mathrm{II}}(z) -
h_{1,1}^{\mathrm {II}}(z)h_{1,0}^{\mathrm{II}}(z) &
h_{1,0}^{\mathrm {II}}(z)h_{0,1}^{\mathrm{II}}(z)
h_{1,1}^{\mathrm {II}}(z)h_{1,1}^{\mathrm{II}}(z)
\end{bmatrix}\\
\end{align}$$

To achieve the zero entries in the identity matrix, it is necessary that

$$\begin{align}
h_{0,1}^{\mathrm {II}} h_{0,0}^{\mathrm {II}}(z)&= h_{0,1}^{\mathrm {II}}(z)h_{1,1}^{\mathrm{II}}(z)\\
h_{1,0}^{\mathrm {II}}(z)h_{0,0}^{\mathrm{II}}(z) &= h_{1,0}^{\mathrm{II}}(z)h_{1,1}^{\mathrm {II}}(z)
\end{align}$$

This breaks down to filters that need to have zeros at every $N$th position, since the 2-polyphased decomposition of a HPF/LPF transformed HPF is either the same as the original's decomposition, or the original's decomposition with a factor of $-1$.


Fast / Efficient Way to Decompose a Separable integer 2D Filter Coefficients without the SVD I would like to be able to quickly determine whether a given 2D kernel of integer coefficients is separable into two 1D kernels with integer coefficients. E.g.

 2   3   2
 4   6   4
 2   3   2

is separable into

 2   3   2

and

 1
 2
 1

The actual test for separability seems to be fairly straightforward using integer arithmetic, but the decomposition into 1D filters with integer coefficients is proving to be a more difficult problem. The difficulty seems to lie in the fact that ratios between rows or columns may be non-integer (rational fractions), e.g. in the above example we have ratios of 2, 1/2, 3/2 and 2/3.

I don't really want to use a heavy duty approach like SVD because (a) it's relatively computationally expensive for my needs and (b) it still doesn't necessarily help to determine integer coefficients.

Any ideas ?

FURTHER INFORMATION

Coefficients may be positive, negative or zero, and there may be pathological cases where the sum of either or both 1D vectors is zero, e.g.

-1   2  -1
 0   0   0
 1  -2   1

is separable into

 1  -2   1

and

-1
 0
 1

I have taken @Phonon's answer and modified it somewhat so that it uses the GCD approach on just the top row and left column, rather than on row/column sums. This seems to handle pathological cases a little better. It can still fail if the top row or left column are all zeroes, but these cases can be checked for prior to applying this method.

function [X, Y, valid] = separate(M)    % separate 2D kernel M into X and Y vectors
  X = M(1, :);                          % init X = top row of M
  Y = M(:, 1);                          % init Y = left column of M
  nx = numel(X);                        % nx = no of columns in M
  ny = numel(Y);                        % ny = no of rows in M
  gx = X(1);                            % gx = GCD of top row
  for i = 2:nx
    gx = gcd(gx, X(i));
  end
  gy = Y(1);                            % gy = GCD of left column
  for i = 2:ny
    gy = gcd(gy, Y(i));
  end
  X = X / gx;                           % scale X by GCD of X
  Y = Y / gy;                           % scale Y by GCD of Y
  scale = M(1, 1) / (X(1) * Y(1));      % calculate scale factor
  X = X * scale;                        % apply scale factor to X
  valid = all(all((M == Y * X)));       % result valid if we get back our original M
end

Many thanks to @Phonon and @Jason R for the original ideas for this.


Tools for 3D shape analysis I have a 3D shape in a 3D binary image. Therefore, I have a list of all of the x,y,z points. If I am to analyze this shape for various identification, such as "sphericity", volume, surface area, etc., what are some of the choices do I have here?

To approximate the volume ($V_{counted}$) you need to count all the voxels.

To roughly approximate the surface area ($A_{approx}$) count all the voxels that have an "empty" voxel as a neighbor. (Have a look at Marching cubes or Marching tetrahedrons for more inspiration and a more detailed discussion on how to determine a surface element).

In order to determine the sphericity I would calculate the minimal bounding box$^1$, take the largest axis and regard it as the diameter: $V_{calc} = \frac{4}{3}\pi r^3$

Since your features ($f_{vol}$, $f_{area}$, $f_{sphere}$, $\dots$) should be in the interval $[0, 1]$ (or 0% - 100%) you simply need to divide the smaller value by the larger value:

$f_{vol}(V_{counted}, V_{ref}) = \frac{\min(V_{counted}, V_{ref})}{\max(V_{counted}, V_{ref})}$

$f_{area}(A_{approx}, A_{ref})$ and $f_{sphere}(V_{counted}, V_{calc})$ accordingly.

$f_{compactness} \tilde{=} f_{sphere}$, since a sphere is the most compact volume.

I hope you can understand, what I am trying to say, and hopefully the answer you are looking for is somewhere in here.

$^1$ if your data is not aligned, you may have to align it first (PCA).


Difference between freqz vs freqresp in scipy? I'm always confused that scipy has

scipy.signal.freqz, but also
scipy.signal.lti.freqresp / scipy.signal.dlti.freqresp

According to the documentation they all compute the frequency response. Surprisingly, the convention for the b coefficients are exactly opposite (as pointed out in this answer as well), and they use a different default of sample points. Other than that, my interpretation of the documentation is that scipy.signal.freqz should be yield the same as scipy.signal.dlti.freqresp. However I can't get the same results trying it out on a small example:

import scipy
import matplotlib.pyplot as plt

filter_coef = [0.2, 0.8]
w1, h1 = scipy.signal.freqz(filter_coef)

dlti = scipy.signal.dlti(filter_coef[::-1], [1]) # accounting for flipped convention
w2, h2 = dlti.freqresp(n=len(w1))

fig, axes = plt.subplots(2, 1)
axes[0].plot(w1, np.abs(h1), label="freqz")
axes[0].plot(w2, np.abs(h2), label="dlti.freqresp")
axes[1].plot(w1, np.angle(h1), label="freqz")
axes[1].plot(w2, np.angle(h2), label="dlti.freqresp")
plt.legend()
plt.show()

Produces:

It looks like they match in terms of magnitude, but for some reason the phase is flipped.

Any ideas what is the difference here? And more generally, when should I use which version of freqz/freqresp?

The equation of the filter is given as

$$H(z) = 0.2 + 0.8z^{-1}$$

This can be equivalently described with positive powers of z as:

$$H(z) = \frac{0.2z+0.8}{z}$$

The two equations above are mathematically equivalent. The function freqz specifies that the numerator and denominator polynomials are to be entered in decreasing (negative) powers of z, while freqresp specifies that they are to be entered in positive powers of z. So the OP should use the first form for freqz and the second form for freqresp. In contrast the OP has flipped the order which is the not the same filter.

In general, in order for the filter to be causal, if expressed in positive powers of z it is implied that it also must be divided by poles at the origin, in this case multiplied by $z^{-1}$.

Doing the following will result in the same answer as the result with freqz:

dlti = scipy.signal.dlti(filter_coef, [1,0]) # accounting for positive power convention

Flipping the coefficients, as the OP has done, is the "reverse filter". The reverse filter has the same magnitude response but has zeros at the reciprocal locations on the z-plane, given by $H(z^{-1})$, giving the result that the OP achieved.


Different ways of decomposing signals into sinusoids For purposes of generating a spectrogram, the signal can be divided into overlapping chunks, each chunk multiplied by a window function and transformed to frequency domain using Discrete Fourier Transform (DFT).

But isn't DFT just one of infinitely many ways to decompose a finite-length discrete-time signal into a sum of sinusoids? What other methods are there to do such a decomposition? How would this choice of method affect the appearance of the spectrogram?

Setting aside the issue of a Spectrogram or a sequence of windowed Fourier transforms.

A single windowed Fourier Transform can be shown to correspond to a Holder p norm minimization where p=2.

One can also use any p norm. so yes there are many possible Fourier types of compositions but the p=2 Holder norm can be bound in terms of the p=1 and p=infinity norms, and vice verse. This implies that these alternate spectra will not be radically different from the conventional DFT ( for discrete time).

A lot of sparse techniques depend on the zero norm.

The conventional Fourier Transform is a direct calculation while the others require iterative optimization.

In some cases like a sparse expansion, it is worth it but in general, the conventional Fourier transform is satisfactory. Also through Parsevals theorem, the conventional transform has a straightforward physical interpretation in terms of power.

Given the Wiener Khinchin theorem for random processes, the other norms would have seem to have less justification.

In a nutshell, yes you can have these other possible expansions but not without a lot of justification that the conventional Fourier series has.

There are cases where other expansions have been shown to be useful but not as a general approach


What is meant by finite action signal In section 1.3 from Signals, Systems and Inference, Alan V. Oppenheim, George C. Verghese. Finite-action signals, which are also called absolutely summable signals, are defined by the condition
\begin{equation}
 \int_{-\infty}^{\infty}|x(t)|dt<\infty
\end{equation}
whereas for discrete time signal, its as
\begin{equation}
 \sum_{k = -\infty}^{\infty}|x[k]|<\infty
\end{equation}
The integration and sum on the left are called the action of the signal. Therefore also known as finite-action signals.
So I would like to know what does finite action signal means? Or is there intuitive explanation about it.

It's something like "finite energy signals"

$$ \int_{-\infty}^{\infty}|x(t)|^2 dt<\infty $$

$$ \sum_{n=-\infty}^{\infty}|x[n]|^2 <\infty $$

except that the "finiteness" is a tougher constraint with

$$ \int_{-\infty}^{\infty}|x(t)| dt<\infty $$

$$ \sum_{n=-\infty}^{\infty}|x[n]| <\infty $$

for example, the signal

$$ x(t) = \frac{1}{\sqrt{1+t^2}} $$

is finite energy, but is not finite action.