Coherent averaging of polyphase components of a signal I want to subsample a multi-band signal to introduce controlled aliasing. My frequency components are uniformly distributed across the spectrum (comb structure) and occupy a relatively narrow bandwidth each. The spectrum can be said to be sparse.

I achieve decimation through picking every $n$-th sample from the original signal. This of course does not take advantage of the remaining samples, which are simply rejected. By adding a phase (sample) shift in my decimation routine, I can get $n$ under-sampled signals in total, which are known as polyphase components. I can also incoherently average them (periodogram averaging) in the spectral domain, since the spectrum is of my interest rather then the time-domain representation.

Is it possible to coherently (in phase) average such undersampled polyphase components of my signal? What phase shift should be applied after down-sampling with a factor $n$? In other words, can we simultaneously undersample (decimate) and oversample (average) to take advantage of the remaining samples for quantization noise suppression?

A typical scenario would be sampling a 1 GHz wide comb signal with a spacing of 10 MHz between the teeth, which we can decimate 10 times to have a spacing of 1 MHz and a bandwidth of 100 MHz. If we sample the signal with an 8-bit analog-to-digital converter, can we use the remaining 9 polyphase components to enhance the vertical resolution of the signal after subsampling?

Thank you,
Lukasz S.

Yes. Simply sum all the polyphase outputs and the sum result will have higher resolution.

Consider that each polyphase output is a delayed version of the same signal, so that if you commutated through all the outputs, you would get a higher sampled version of your same signal and the quantization noise of this signal would be approximately white across this higher digital frequency band. If you passed this signal through a low pass filter prior to taking every nth sample in a decimation process you would indeed benefit by an increased dynamic range to the extent the quantization noise components in each sample are independent. (they should be; the delayed replicas at each output would have independent quantization noise and coherent signal components). So the decimation filter could be a simple moving average in which case it is easy to see that your best solution is to sum average the outputs of your polyphase filters and use that average as the final result rather than just selecting any one of them.  The moving average introduces passband droop- so with more complication you could do more elaborate filtering. I would be inclined to use a CIC filtering approach as the moving average followed by a simple 3 tap inverse sinc shaper if I was concerned about the resulting passband droop.


What is signal-energy? What is energy in terms of signal processing? You can read here, that energy is

I just need a simple everyday-example for that energy. And yes, bringing an example like

  $x(t)$ representing the potential of an electrical signal propagating
  across a transmission line

is still too complicated for me.

Pick up an apple and raise it one meter. You have just spent (approximately) one joule of energy. You had to spend energy because the action (raising the apple) goes against a force (Earth's gravity).

The same is true for electrons. An electrical signal is nothing but electrons moving in an electric field. Causing these electrons to move takes energy, exactly like moving a train or an apple.

In the formula you posted, $x(t)$ is a voltage causing electrons to move in a resistance of $1\,\Omega$, and $E_s$ is the energy required to make those electrons move around in the resistor, causing a current to develop.


What is meant by "correlation" when referring to spectral coherence I've been reading that coherence measures the correlation between two waves as a function of frequency. I also read that difference in phase does not mean less coherence at a given frequency, and that amplitude is also not a factor. Some explanations say that its is a measure of the shift in phase. As I understand it, a shift in phase requires a shift in frequency, but coherence is the measure of correlation at a given frequency. So I guess what I want to know is this:

If I've taken the Fourier transform of two waves, and am looking to see how they are correlated in the frequency domain, what am I actually measuring?

EDIT: I'm a biologist, and my knowledge of digital signal processing is pretty basic. I'm trying to pick this stuff up on the fly. So intuitive explanations, if possible, would be greatly appreciated.

The coherence function, as used in signal processing, measures the normalized correlation between to power spectra:

$$
C_{xy}(f) = \frac{|G_{xy}(f)|^2}{G_{xx}(f) G_{yy}(f)}
$$

where $G_{xx}$ is the power spectral density (PSD) of the signal $x(t)$, $G_{yy}$ is the PSD of $y(t)$, and $G_{xy}$ is the cross-spectral density (CSD) of $x(t)$ and $y(t)$.

The CSD is computed as the Fourier transform of the cross correlation of $x$ and $y$:

$$
R_{xy}(\tau) = \operatorname{E} \left[ x(t) y(t + \tau) \right]
$$
$$
G_{xy}(f) = \mathcal{F} \left[ R_{xy} (\tau) \right]
$$

You can think of the CSD as measuring "hidden periodicities" in the cross correlation, just like the PSD measures periodicities in the autocorrelation. (This is where the name for the "periodogram" PSD estimator comes from.)

Incidentally, if the cross-correlation has harmonic content at some frequency $f$, you can think of the signals being correlated at that frequency.


Impulse response of a 3x3 PSF - how to find analytical expression for fourier transform of a 3x3 matrix? I have a filter $\mu[n_1, n_2]$ with taps:

$$ (1/8) (1/4) (1/8)$$
$$ (1/4) (1/2) (1/4)$$
$$ (1/8) (1/4) (1/8)$$

How do I find an analytical expression for $\hat\mu(w_1, w_2) $?

Since it looks so much like a triangle, I feel like it looks something of the form $sinc^2(w)$, but the exact answer I need help with.

You can use the 2D DFT formula

$$H(\omega_1,\omega_2) = \sum_{n_1} \sum_{n_2} h[n_1,n_2] e^{-j(\omega_1 n_1 + \omega_2 n_2)} $$

and simply the trigonometric algebra to get a closed form analytic expression for the 2D-DTFT. However, as @LaurentDuval has already mentioned, your 3x3 kernel is separable and one set of 1D filters is this

$$f[n_1]  = [\frac{1}{4} , \frac{1}{2}, \frac{1}{4}]^T $$
$$g[n_2]  = [\frac{1}{2} , 1 , \frac{1}{2}] $$

Then from the Fourier transform we know that

$$ h[n_1,n_2] = f[n_1]g[n_2] \implies H(\omega_1,\omega_2) = F(\omega_1) G(\omega_2) $$

Analytic expressions for 1D-DTFT's can be obtained from
$$H(\omega) = \sum_{n}h[n] e^{-j \omega n} $$

Assuming 3x3 kernel has zero phase, then we see that :
$$F(\omega_1) = \frac{1}{4} e^{j \omega_1 } + \frac{1}{2} e^{j 0 } + \frac{1}{4} e^{-j \omega_1 } = 0.5 \cos(\omega_1) + 0.5 $$

similarly
$$G(\omega_2) = \frac{1}{2} e^{j \omega_2 } + 1 e^{j 0 } + \frac{1}{2} e^{-j \omega_2 } = \cos(\omega_2) + 1 $$

Hence:
$$H(\omega_1,\omega_2) = 0.5 \left( \cos(\omega_1) + 1 \right)(\cos(\omega_2) + 1) $$


Nyquist Criterion in higher dimensions I was just wondering if there were any "obvious" resources, papers, or textbooks I should look into regarding the nyquist-criterion in higher dimensions. For example, this would be relevant if there is some curve that is i.e. not 1-to-1 w.r.t. time: $f: t \in \mathbb{R}^1 \rightarrow f(t) \in \mathbb{R}^1 $, but instead for some other parameterization, we have a line $l: u \in \mathbb{R}^1 \rightarrow l(u) = \big(\tilde t(u),\tilde f(u)\big) \in \mathbb{R}^2$, it is 1-to-1, but not onto, as $(\tilde t,\tilde f)$ is on a 1-dimensional manifold (a line) in a 2-D plane.

Then I'd like to know what sampling scheme I need to be able to recreate the curve through reconstruction. Using this approach, you could use nyquist's theorem for general spatial surfaces and manifolds as well.

This might be a math question, so I can migrate this. But I figured some experts on this sub-stack might know.

EDIT: several small edits for clarity/correctness.

EDIT 2: Fixed parameterization!!!

EDIT 3: I want to emphasize that although I am parameterizing by $t$, it may be more useful to think of $t$ as a spatial dimension, since time can only move forward.

i do not have a reference for you, but you are sampling among every one of the dimensions and the Nyquist Criterion must be satisfied for each one.

so fix $v$ to an arbitrary value $v_0$, then check if $\tilde{f}(u,v_0)$ is bandlimited some bandlimit, we'll call $B_u$.  find the max value of $B_u$ for all possible $v_0$, double that and that is the minimum sample rate along the $u$-axis.

then swap the roles of $u$ and $v$, fix $u$ to $u_0$, find the bandlimit of $\tilde{f}(u_0,v)$ for every possible value of $u_0$.  find the maximum (over all $u_0$) bandlimit $B_v$, double that and that is the minimum-sample rate along the $v$-axis.

if your curve can move at an arbitrary angle w.r.t. the $u$-axis and $v$-axis, then pick the maximum of $B_u$ and $B_v$ and sample both axes at double of that rate (or slightly higher).  then, along the path length, you must sample at a minimum of that rate, per unit of arc length, to satisfy Nyquist in two dimensions.


QPSK sync sequence detection problem I have IQ data collected

This is the beginning of the packet that includes preamble and sync sequence. Red is real part, blue is imaginary.
Also constellation diagram for this data

I don't have access to the radio that transmit this data.
I know that this is QPSK modulated data. Preamble and sync sequence should be BPSK modulated, but I cannot demodulate this data, they don't look like a "good" BPSK. I also know that the data have been processed by RRC filter before transmission.
What is this phenomenon, and what I can do with this data?
Thank you,

In this case you have very clear symbol transitions and the appearance of RRC shaped symbols with a small frequency offset. What you could do is pass the waveform as is through a second RRC filter (matched filter) assuming you know the RRC parameters and then estimate the decision locations in the preamble (since simpler there with just two symbols) as the midpoint between transitions where a 1010 pattern exists and observe successive complex decisions for rotation; from inspection it appears that rotation will indeed be less than 90 degrees from symbol to symbol in which case you should be able to determine the frequency offset and correct that with an opposite complex rotation (stop the symbol decisions from rotating).  Once that is done all successive symbols at the decision location should be clustered around their expected positions.
In practice a timing recovery loop (Gardner loop for example) can determine the symbol locations and a carrier recovery loop can remove the carrier offset; along with equalization when significant multipath interference exists- but this waveform looks to be of very good SNR with little distortion that doing it manually as described above will be a good learning exercise as to what those algorithms need to do in a tracking and acquisition loop.


Distinguishing Spoken Speech from Played-back Recording of Speech Do any general-purpose techniques exist to distinguish live spoken speech from played-back speech samples? For example, if I were to deploy a speaker verification system to authenticate users by voice, I would not want the system to be able to be spoofed by playback of voice recordings.

I think I can imagine at least one technique which might work for certain known combinations of microphones / speakers: If a microphone of sufficient quality captures a novel signal outside the speech band that corresponds (in time) to a signal within the speech band, this would seem to indicate a high likelihood of the sample not being a live speaker (because live speakers don't emit signals outside the speech band). Is this reasoning sound? (No pun intended).

Do any general-purpose techniques exist to distinguish live spoken
  speech from played-back speech samples?

No, at least not how I understand your definition.

Consider the case of  "live speech" : some one talks into the microphone, you digitize the signal and analyze it. You can also store the digitized speech and analyze it a day later. Now the exact same samples constitute "play back". Since the samples are identical, you can't distinguish them from "live speech".


Discrete time Final Value Theorem applied to feedback system I wish to calculate the Final Value of systems in which a high pass filter of the output feeds back into the input.
A simple example would be:

where  is a 1st order high pass filter with transfer function:

I was expecting the y in the above example to have an infinite final value to a step in x, because  keeps feeding
However, the workings below give a different answer:

Re-writing hp1(z) in terms of its inputs only:

Add  to both sides of the system's equation:

Write the system's transfer function:

Re-write the infinite sum in the denominator:

Apply the Final Value Theorem to the response of this system to a step in x:

Taking the limit:

The above suggests that the system  has a well defined terminal value to a step in x. However I don't think that can be the case.
Where am I going wrong?
Help much appreciated

I would derive the total transfer function directly in the transform domain. Your input-output equation can be written as
$$Y(z)\big(1-z^{-1}\big)=\alpha G(z)z^{-1}Y(z)+\beta z^{-1}X(z)\tag{1}$$
where $G(z)$ is the transfer function of the high-pass filter. From $(1)$ you directly obtain the transfer function
$$H(z)=\frac{Y(z)}{X(z)}=\frac{\beta z^{-1}}{1-z^{-1}\big(1+\alpha G(z)\big)}\tag{2}$$
which clearly has a pole at $z=1$ because $G(1)=0$.

The mistake in your calculation is in step 2. You increased the lower bound of the summation from $1$ to $2$ without changing the power of $\phi$ from $i$ to $i-1$. In the transfer function this results in a $\phi^2$ in the last part of the denominator instead of a $\phi$. The denominator should be:
$$D(z)=1 - (1+\Psi\phi)z^{-1}-\Psi(\phi-1)\frac{\phi z^{-2}}{1-\phi z^{-1}}\tag{3}$$
For $z=1$ this evaluates to
$$D(1)=1 - (1+\Psi\phi)-\Psi(\phi-1)\frac{\phi }{1-\phi }=1 - (1+\Psi\phi)+\Psi\phi=0\tag{4}$$


Why is the median filter called 'edge-preserving'? I want to know about median filtering. Recently, I read that it is best for preserving edges. But how does it preserve edges? I need complete knowledge of median filters. Why it is non-linear filter and why it is best, compared to linear filter?

Non-linearity

A linear filter is mathematically described by the convolution sum (for discrete signals) and the convolution integral for continuous signals. The median cannot be found using a linear function except in the trivial case where you have a discrete filter of size 1, which is why the median filter is non-linear.

Edge Preserving Properties.

  Median filtering is one kind of smoothing technique, as is linear
  Gaussian filtering. All smoothing techniques are effective at removing
  noise in smooth patches or smooth regions of a signal, but adversely
  affect edges. Often though, at the same time as reducing the noise in
  a signal, it is important to preserve the edges. Edges are of critical
  importance to the visual appearance of images, for example. For small
  to moderate levels of (Gaussian) noise, the median filter is
  demonstrably better than Gaussian blur at removing noise whilst
  preserving edges for a given, fixed window size. However, its
  performance is not that much better than Gaussian blur for high levels
  of noise, whereas, for speckle noise and salt and pepper noise
  (impulsive noise), it is particularly effective. Because of this,
  median filtering is very widely used in digital image processing.

Source: Wikipedia


Doubt : Derive DFT from CFT In the below screenshot I was unable to understand how the equation 4.3 is derived from the above equation.. How $x(nT)$ has become $x(n)$ ? Please clear my doubt and suggest some link to understand DFT. I am sorry for the poor image quality.

It looks like the author implicitly makes the substitution $x(nT) \rightarrow x(n)$, but even if they were making this explicit, there are more clear ways to get to Eq. 4.3.

One way to do this is to denote the continuous (analog) and discrete time-domain signals as $x_a$ and $x$, respectively. Your text, using this convention, would then look like this (starting from Eq. 4.2):

$$
x_a(t) \xrightarrow{sampling} x_a(nT)
$$

and

$$
X(e^{j\omega})\ =\ \sum_{n = -\infty}^{\infty} x_a(nT) e^{-j\Omega nT}.
$$

(Note the right-hand-side exponent: it is $-j\Omega nT$, not $-j\Omega nt$ as it appears in your text; this looks like an error to me.) Then if you take into account the following equality (assuming uniform sampling):

$$
x(n)\ =\  x_a(nT), \quad -\infty < n < \infty,
$$

you end up with:

$$
X(e^{j\omega})\ =\ \sum_{n = -\infty}^{\infty} x(n) e^{-j(\Omega T)n},
$$

or (noting that $\omega = \Omega T$):

$$
X(e^{j\omega})\ =\ \sum_{n = -\infty}^{\infty} x(n) e^{-j\omega n}.
$$

By the way, this is a derivation from CFT to DTFT operation, not CFT to DFT one (as you say in the title of your question).


Double the "speed", double the frequency? Usually, when some sample is played-back at twice the rate it was recorded, the pitch of the recording sounds an octave higher.

For example: $g(t) = \sin(2 \pi f t)$ and $g(2t) = \sin(2 \pi f 2 t) \implies \sin(4 \pi ft)$. Where $f$ is a given frequency.

There are some strange "sample" cases however that when played back at twice their recording rate their pitch appears to...not have changed as expected. And some even stranger "sample" cases that when played back at twice their recording rate their pitch appears to be...lower!

How do these strange "sample" cases might "look like" (in terms of a mathematical description) and what are, broadly, the conditions that can lead to such cases?

EDIT (hint 1):

Variants of the technique of pitch-shifting / time stretching have been mentioned so far by more than one persons in relation to this question.

Pitch shifting is not the core of the question and the technique of pitch-shifting does not change the nature of the question.

As it was pointed out in one of the comments to this post and as a form of a hint: Yes, pitch shifting does alter the pitch of the "sample" but a) what sort of sample would sound almost identical to its unshifted version if it was passed through a $2 \times$ pitch shifter and b) what sort of sample would sound at a lower pitch if it was passed through a $2 \times$ pitch shifter. (the $2 \times$ here is just an example).

EDIT (hint 2):

Another two things that have appeared in responses so far are aliasing and some aspect of human hearing. In other words, "Maybe you have pitch shifted in such a way that the frequency jumps back on to itself" and also "Maybe you are taking advantage of some psychoacoustic phenomenon in which the frequency has changed but the perceived frequency does not"

Aliasing and the limits of human hearing are relevant to this question but are not absolute requirements for the "paradox" that is mentioned here to occur.

White noise: sounds exactly the same if you play it back at double the speed
Stack up diminished fifth starting at a low enough fundamental (say 20 Hz). Play back at double speed: sounds exactly the same.
Stack octaves with equal amplitude starting with a low enough fundamental. Play back at twice the speed: sounds exactly the same
Detune the octaves slightly (say multipler of 2.05 starting at 20 Hz). Play back at double the speed: perceived pitch goes DOWN.

You get plenty more when fiddling with the harmonic structure.
Code example for the last one:

%% make some disharmonic sounds
fs = 44100; % sample rate
q = 2.05; % harmonic multiplier
f0 = 20;  % fundamental
n = 10; % number of harmonics
duration = 2; % length of file at fs
nx = duration*fs;
t = (0:nx-1)'/fs; % time vector
x0 = sin(2*pi*f0*t); % fundamental
f1 = f0;
for i = 2:n  % add harmonics
  f1 = f1*q;
  x0 = x0+sin(2*pi*f1*t);
end

x0 = 0.5.*x0/max(abs(x0)); % normalize
% and play
sound(x0,fs); pause(duration);
sound(x0,2*fs);


Does $y[n] = x[n] \star (u[n]-u[n-2])$ have memory or is it memoryless? $y[n] = x[n] \star (u[n]-u[n-2])$, by its definition is has to be a system with memory since it is depended from a fraction of time in the past, but if we calculate the difference $u[n]-u[n-2]$, it results $u[n]$ for $0 \le n < 2$, which is memoryless. And that sign $\star$ is convolution.

The system $$y[n] = x[n] \star (u[n]-u[n-2])$$ where $u[n]$ is the unit step function, has memory.

Indeed the system is equivalent to $$y[n] = x[n] \star ( \delta[n] + \delta[n-1] ) \implies y[n] = x[n] + x[n-1]$$ and as it's clear from the given I/O relationship, the current value of the output $y[n]$, depends on the values input $x[n]$ at other times namely $x[n-1]$ (one step past value); hence the system has memory.


High FIR order in time critical classification task I have a question regarding filtering (bandpass in my case) as a preprocessing step before a classification task.
The frequency band I am interested in lies between 0.5 - 22 Hz which results in a 6000th order bandpass filter due to the steep slope needed for the highpass.
The phase response is linear and I use the filtfilt matlab function to get a zero phase output.
After filtering I try to classify on sample level, to determine when a certain state of two is present in my data.
My question now is whether the high filter order distorts the time sequence, since the output samples of the filter are a superposition of samples $n$ to $n-5999$?

I would handle this as a chain of FIR filters with decimation to get much lower tap count FIR filters that don't need to be nearly as steep in the early stages.
Maybe using Blackman-Harris windowed filters for the first two stages (for superior alias fold-in rejection).
Then a Kaiser windowed filter for the 3 stage to get a flatter top and a sharper roll-off.
I guarantee you'll be using much less than 6000 taps total for the same filtering performance.
I've worked out fairly decent scheme in some Octave/MatLab code:
LPF -> decimate by 2 -> LPF -> decimate by 5 -> BPF
This gets you to a final 100 samples/second sample rate.
The three filters, that I came up with, have tap lengths of 9, 51, and 205 respectively.
The two low-pass anti-alias filters keep any alias fold-in into your band of interest at -40 dB or lower.
The final band-pass filter is -35 dB at DC, -6.3 dB at 0.5 Hz, -0.1 dB at 22.5 Hz, and -40 dB at 24.5 Hz.

Here's the Octave code.
Main script:
% First anti-alias filter, designed for decimation by 2 after the filter
Fs1 = 1000;  % [-500,500]
a = blackman_harris_lpf(1.0, Fs1, 100.0, 300.0, 55);
length(a)

% Second anti-alias filter, designed for decimation by 5 after the filter
Fs2 = Fs1/2;  % [-250, 250]
b = blackman_harris_lpf(1.0, Fs2, 50, 20.0, 45);
length(b)

% Band pass filter for the frequency band of interest
Fs3 = Fs2/5;  % [-50, 50]
c = kaiser_lpf2(1.0, Fs3, (23.5-0.5)/2, 2, 9);
radian_phase_inc = pi/(Fs3/2) * (24.5-0.5)/2;
c = 2*real(c .* exp(1j * radian_phase_inc .* [0:(length(c)-1)]'));
length(c)

% Get response curves for each filter for plotting
[ha,wa] = freqz(a, [1], 2*5*1000);
[hb,wb] = freqz(b, [1], 2*1000);
[hc,wc] = freqz(c, [1], 2*1000/5);

% Create plotting x-axis (aka frequency-axis) value arrays to
% fold freq axis to show aliasing into final passband.
% Scale to original Fs1
waturn = [wa(1:(2*5*1000/2))' wa((2*5*1000/2):-1:1)'] * 500/pi; % Fs = 1000, fold for decimation by 2
wbturn = [wb(1:(2*1000/5))'   wb((2*1000/5):-1:1)' ...
          wb(1:(2*1000/5))'   wb((2*1000/5):-1:1)' ...
          wb(1:(2*1000/5))'                         ] * 250/pi; % Fs = 500, fold for decimation by 5
wcturn = wc * 50/pi; % Fs = 100, no turns

clf;
plot( ...
     waturn, 20*log10(abs(ha)), ';1st filter, Blackman-harris LPF, 1000 sps folded to 500 sps;', ...
     wbturn, 20*log10(abs(hb)), ';2nd filter, Blackamn-harris LPF,  500 sps folded to 100 sps;', ...
     wcturn, 20*log10(abs(hc)), ';3rd filter, Kaiser BPF, 100 sps;'  ...
);
title('Filters Alias Folding Performance');
ylabel('Gain (dB)');
xlabel('Frequency (Hz)');
grid on;
hold off;

blackman_harris_lpf.m:
function b = blackman_harris_lpf(gain, Fs, Fc, W, atten_dB)
% B = blackman_harris_lpf(GAIN, FS, FC, W, ATTEN_DB)
%
% Design a Blackman-harris windowed low pass filter
%
%    GAIN Passband linear gain. E.g. 1.0
%
%    FS   Sampling frequency in Hz
%
%    FC   Cutoff frequency (the edge of the ideal LPF response before
%         before truncation and windowing of the sinc() pulse) in Hz.
%
%    W    Transition width in Hz
%
%    ATTEN_DB  stop band attenuation in dB

M = floor(atten_dB * Fs/(22.0*W)); % rule of thumb from fred harris
if mod(M,2) == 0  % ensure M is odd
   M = M + 1;
end

% Windowed, truncated sinc()
w = blackmanharris(M);
h = sinc([-(M-1)/2:1:(M-1)/2]*Fc/(Fs/2)).' .* w;

% Normalize passband response to 0 dB and apply gain
b = gain * h / sum(h);

kaiser_lpf2.m:
function b = kaiser_lpf2(gain, Fs, Fc, W, Beta)
% B = kaiser_lpf2(GAIN, FS, FT, W, BETA)
%
% Design a Kaiser windowed low pass filter
%
%    GAIN Passband linear gain. E.g. 1.0
%
%    FS   Sampling frequency in Hz
%
%    Fc   Cutoff frequency (the edge of the ideal LPF response before
%         before truncation and windowing of the sinc() pulse) in Hz.
%
%    W    Transition width in Hz
%
%    BETA Kaiser window parameter dictating the width of the main lobe
%         and the height of the stopband lobes.  Beta == 0 is equivalent
%         to a rectangular window.

M = floor((Beta/0.1102 + 8.7) * Fs/(22.0*W)); % rule of thumb from fred harris
if mod(M,2) == 0  % ensure M is odd
   M = M + 1;
end

% Windowed, truncated sinc()
w = kaiser(M, Beta);
h = sinc([-(M-1)/2:1:(M-1)/2]*Fc/(Fs/2)).' .* w;

% Normalize passband response to 0 dB and apply gain
b = gain * h / sum(h);


What are the most common algorithms for adaptive thresholding? Adaptive thresholding has been discussed in a few questions earlier:

Adaptive Thresholding for liver segmentation using Matlab

What are the best algorithms for document image thresholding in this example?

Of course, there are many algorithms for Adaptive thresholding. I want to know which ones you have found most effective and useful.

Which Adaptive algorithms you have used the most and for which application; how do you come to choose this algorithm?

I do not think mine will be a complete answer, but I'll offer what I know and since this is a community edited site, I hope somebody will give a complimentary answer soon :)

Adaptive thresholding methods are those that do not use the same threshold throughout the whole image.

But, for some simpler usages, it is sometimes enough to just pick a threshold with a method smarter than the most simple iterative method. Otsu's method is a popular thresholding method that assumes the image contains two classes of pixels - foreground and background, and has a bi-modal histogram. It then attempts to minimize their combined spread (intra-class variance).

The simplest algorithms that can be considered truly adaptive thresholding methods would be the ones that split the image into a grid of cells and then apply a simple thresholding method (e.g. iterative or Otsu's method) on each cell treating it as a separate image (and presuming a bi-modal histogram). If a sub-image can not be thresholded good the threshold from one of the neighboring cells can be used.

Alternative approach to finding the local threshold is to statistically examine the intensity values of the local neighborhood of each pixel. The threshold is different for each pixel and calculated from it's local neighborhood (a median, average, and other choices are possible). There is an implementation of this kind of methods included in OpenCV library in the cv::adaptiveThresholding function.

I found another similar method called Bradley Local Thresholding. It also examines the neighborhood of each pixel, setting the brightnes to black if the pixels brightness is t percent lower than the average brightness of surrounding pixels. The corresponding paper can be found here.

This stackoverflow answer mentiones a local (adaptive) thresholding method called Niblack but I have not heard of it before.

Lastly, there is a method I have used in one of my previous smaller projects, called Image Thresholding by Variational Minimax Optimization. It is an iterative method, based on optimizing an energy function that is a nonlinear combination of two components. One component aims to calculate the threshold based on the position of strongest intensity changes in the image. The other component aims to smooth the threshold at the (object)border areas. It has proven fairly good on images of analog instruments (various shading and reflection from glass/plastic present), but required a careful choice of the number of iterations.

Late edit:
Inspired by the comment to this answer. There is one more way I know of to work around uneven lighting conditions. I will write here about bright objects on a dark background, but the same reasoning can be applied if the situation is reverse. Threshold the white top-hat transform of the image with a constant threshold instead of the original image. A white top hat of an image is nothing but a difference between the image $f$ and it's opening $\gamma(f)$. As further explanation let me offer a quote from P. Soille: Morphological Image Analysis:

  An opening of the original image with a large square SE removes all relevant image structures but preserves the illumination function. The white top-hat of the original image or subtraction of the illumination function from the original image outputs an image with a homogeneous illumination.


Reducing one frequency in song How would I take a song input and output the same song without certain frequency ranges?

Based on my research so far, the song should be broken down into chucks, FFT it, reduce the target frequency ranges, iFFT it, and stitch the chunks back together. However, I am unsure if this is the right approach to take, and if so, how I would convert from the audio to FFT input (what seems to be a vector matrix), how to reduce the target frequency ranges, and how to convert back from the FFT output to audio, and how to restitch.

So far, I've understood the basic premise of fft (through basic 3blue1brown yt videos and the like) and that it is available through scipy/numpi, and figured out how to convert from youtube to 0.25 second chunks in a wav format. I've tried googling the problem, but while many people say it's trivial, I cannot find an explanation with psudocode or actual code.

For background, my grandpa loves music. However, recently he cannot listen to it, as he has become hypersensitive to certain frequency ranges within songs. I'm a high school student who has some background coding, and am just getting into algorithmic work and thus have very little experience using these algorithms. Please excuse me if these are basic questions; any pointers would be helpful.

Zeroing out bins / attenuating them in discrete Fourier domain is universally a bad idea, due to the undesirable time-domain effects of that.

Instead, use a audio processing program to apply an adjustable equalizer to the song of choice. Let's walk you through Free software:

Get audacity;
load your song in that, open "Effects"->"Equalizer"
start with one  of the presets to see how things work, and then
sit down with your grandpa and figure out what frequency ranges to dampen.
Save that preset!
You can later apply it to all your music collection using batch processing.

You should also, in that Equalizer dialogue, "manage / export" your freshly designed filter frequency response. With a little work, you could, for example, build a digital equalizer (in music, your granddad probably won't mind the latency that brings, which would be a bad thing in live music) using e.g. a Raspberry Pi and a USB soundcard, which you could put in between you granddad's CD player and his amplifier, or so. (of course, it could also act as media player/server/center/whatevs.)


Percentage difference between two signals Apologies in advance for what I imagine is a simple question, however I can't work out how to do this.

I have two signals of the same type but with different sensor types, the figure below will give you an idea of some of the data I've got:

I very simply want to calculate the average percentage difference between the two signals. Methods I've attempted so far have calculated outrageous results, in the order of ~200% difference. My attempt is below:

  for i = 1:length(data)
      per(i) = abs((var1(i)-var2(i))/((var1(i)+var2(i))/2)) * 100;
  end

I believe the problem in this method is that I'm only calculating the difference between the two immediate points in question (i) and not on a grand scale of all the data.

Any help would be greatly appreciated.

At the least you should normalize it by the length of the sequence. But what you're trying to do is a pretty standard thing. I'm not sure what percentage difference really means. What you should probably be using is something like the root mean squared error:

$$RMSE=\sqrt{\frac{1}{N}\sum_{n=0}^{N-1}(x(n)-y(n))^2}$$

where $x(n)$ and $y(n)$ are your two sequences.


Removing Gaussian Noise from a Signal to Get Minimum Value I have a signal that has a minimum value that I'm trying to read.
The issue I'm having is that the signal is spread out by gaussian noise.
I have the signal at a lot of timesteps (and expect the function to be roughly exponentially distributed and being truncated at the minimum value), but can't work out how to remove that gaussian noise.
Any help would be appreciated!
(Also this may need transplanting to another stackexchange, but I feel like this was something which should have been solved multiple times previously by this community)

If you have no prior data on the signal of interest there is nothing to do actually.
The more prior you have the better you can do.
For instance, if the only information you have is the Bandwidth of your signal, the best you can do is apply an Band Pass / Low Pass Filter.
If you know a sparse representation of your signal it will be great, as you can remove all the other coefficients (Generalization of the LPF / BPF).
In case you know something about how smooth it is, you can apply Total Variation like smoothing.
Remark
When you write Gaussian Noise I assume you mean Additive White Gaussian Noise (AWGN).


One representative time series signal from multiple time series data I came across this basic problem in DSP which is not my area owf expertise but I consider this problem to be a common problem which others might be addressing.

I have multiple time series data and I need to find a time series data that is a representative of all those data. One method that I am applying is averaging the signals but I am sure there will be some other sophisticated methods to do the same.

What are the other methods that I can apply?

Thanks

Decision of representative is usually formulated as an optimization problem. The optimal representative minimizes sum of some dissimilarity measures from the observed signals. The average you indicated is optimal in the sense of minimization of Euclidean distance. If you want to minimize feet of perpendiculars, you can use the first principal component.

Now the problem is how to choose the dissimilarity measure. In signal processing, we often consider some probabilistic model and use the dissimilarity measure proportional to the negative log likelihood. This is often referred to as maximum likelihood estimation. For example, Euclidean distance corresponds to the negative log likelihood of mean parameter of Gaussian distribution. Even using the same Gaussian distribution, there can be various ways of modeling, and various dissimilarity measures are derived such as Mahalanobis distance and Itakura-Saito distance.


Raised-Cosine Filter - calculating magnitude response Here is the formula to compute magnitude of a raised-cosine filter:

I know that 0 <= Beta <= 1. I am not sure what are the values for f (frequency) and Ts. I assume that some normalization will take place. I also know that Ts is symbol period and it is reciprocal of symbol rate. What are these parameters and relationship between them? Can someone provide an example of typical values of these parameters as used in the raised-cosine formula? Here is the graph from Wikipedia. How the normalization of frequency and magnitude is done?

The frequency $f$ is no parameter but it's the independent variable. The frequency response $Z(f)$ is a function of frequency. The only real parameter is $\beta$, defining the excess bandwidth, but you already know that $0\le \beta\le 1$. You also know that $T_s$ is the symbol period, so $1/T_s$ is indeed the baud rate, which is part of the system's specifications.

From the formula you can see that the minimum bandwidth is half the baud rate (for $\beta=0$, i.e. zero excess bandwidth). For $\beta=1$ (i.e., $100\%$ excess bandwidth), the bandwidth is twice as large, i.e. it equals the baud rate $1/T_s$.


Why is the bandwidth doubled with AM? Let's say I have an arbitrary signal which has the following frequency spectrum and a bandwidth of BNF:

Now I use Amplitude Modulation to transport the signal. The frequency spectrum will change accordingly:

Now I know that the required bandwidth doubles. I don't understand where the doubled bandwidth applies as I can't transport negative frequencies.

What you've observed is just the difference between baseband modulation (your original plot) and carrier modulation (the second plot). Note that for the AM case, since the signal that modulates the carrier is real, the resulting carrier-modulated signal is symmetric. This means that you can get away with only transmitting one of the sidebands. This approach is called single sideband (SSB) modulation. With an SSB system, the carrier-modulated signal only has the bandwidth of approximately one of the sidebands of the original baseband signal.


How to apply Hamming Window and compute FFT? I'm borrowing an EEG device whose accompanying software also computes the FFT of the data it collects. The software provides the user with a file of EEG values and a file of FFT values computed from the EEG file, and I wanted to learn how to compute the FFT values on my own.

Here's the manufacturer's description of how they do it: "Each path contains 129 decimal values with a range of roughly -40.0 to 20.0. Each array represents FFT coefficients (expressed as Power Spectral Density) for each channel, for a frequency range from 0hz-110Hz divided into 129 bins. We use a Hamming window of 256 samples(at 220Hz), then for the next FFT we slide the window 22 samples over(1/10th of a second). This gives a 90% overlap from one window to the next. These values are emitted at 10Hz."

Going through MATLAB's documentation I ended up with this code:

X = data(:, 1);
h = hamming(256);
buf = buffer(X, 256, 234);
hWin = buf.*repmat(h(:), [1, size(buf, 2)]);
Y = fft(hWin);
P2 = (abs(Y))';
P1 = P2(:, 1:256/2 + 1);
P1(:, 2:end-1) = 2*P1(:, 2:end-1);
P1 = 10 * log10(P1);

But it doesn't produce the sames values that I see in the four FFT files given to me by the manufacturer. First, let me start by saying that their EEG file contains a 193011 X 4 matrix of data, and each column corresponds to the data collected by one of the four channels/sensors on the EEG device. And with each column of data, they were able to compute four separate files, which they named raw_fft0, raw_fft1, raw_fft2, and raw_fft3. Each raw_fft file contains a 8607 X 129 matrix whose values should be somewhere between -40 and -20. The actual observed values for their FFT files are somewhere between -30 to -20 to the low-to-mid 40s.

On the other hand, my FFT files are of the dimension 8774 X 129, and the range of values in my files are bigger than what the manufacturer says they should be. Would anyone mind pointing out what I've done wrong or pointing me to information that could clarify what I've done wrong?

Three possible issues:

The Discrete Fourier Transform definition that the MATLAB fft uses does not preserve the energy of the vector.  Multiply it by $1/\sqrt{N}$, where $N$ is the FFT size, if you want to keep the energy constant.  The DFT used by the other code might (or might not) have different scaling.
You are computing $10log_{10}$ of the magnitude of the signal, but to compute the power spectral density you should be computing $10log_{10}$ of the magnitude squared - or, alternatively, $20log_{10}$ of the magnitude.
I don't know how the other code got 8607 rows, but if you form windows such that each window contains data in the original vector, you would get 8762 rows by my calculation.  In your use of buffer, you get some windows that are stuffed with 0s at the end.


What is the time it takes to measure the SINR in LTE? I would like to know how much time it takes the receiver to evaluate the SINR?

Scenario:
Assume that a base station (eNodeB) would like to transmit to one user. Let us say that it sends a pilot signal of L bits to that user. Then the user calculates the SINR (assuming that there are a set of N other interfering base stations) and compares it to a certain threshold and then feeds back 1 bit to the base station to tell it if the SINR is above or below the threshold. How much time it takes for this process?

In general, how the standard LTE works to evaluate the SINR at the receiver?

So you're giving your answer yourself in the question:

  In general, how the standard LTE works to evaluate the SINR at the receiver?

The point is that measures like SNR and SINR are dependent on what you consider your signal of interest, what you consider interference and what you consider noise¹.

LTE actually specifies multiple different things that describe "received signal strength", whcih you can pu in relation to unwanted effects, and all of these could be labeled as SINR (but have different values for the same transmission!). The two "S" in ("SINR") measurements you probably care about are

RSRP, Reference Signal Received Power: Simply put, LTE uses distributed pilot symbols to figure out the channel state; just adding up the strengths of these will give you the magnitude of the impulse response over all OFDM subchannels.
RSRQ, Reference Signal Received Quality: The reference tone powers above don't say overly much about how well-suited the channel actually is – they are really just a sum. But in an OFDM system, where all but let's say 5% of subchannels are heavily attenuated, but these 5% are extremely well power-conserving, you just won't get the same performance in a system where all channels are attenuated, but things are more evenly distributed, even if the two systems have the same RSRP; hence, RSRQ takes other things into consideration, too.

You can (weightedly) subtract these from RSSI, which is also specificially defined for LTE, to get an estimate for the "IN" part of "SINR".

Frankly, RSRQ and RSRP are probably much better suited to describe a complex multicarrier system than the simplistic SINR – so if I could, I'd just try to argue based on these, and drop the SINR terminology.

¹ This is actually the answer I most often give at my job, "SNR depends on what you consider to be signal and noise. No general formula can be given.".


Why is phase range between $-\pi$ and $+\pi$ (instead of $0$ and $2\pi$)? I wonder why the phase range in DSP applications is between $-\pi$ and $+\pi$.
For example, why is that not between $0$ and $2\pi$?
This question popped when I was reading about synchronization in digital communications.
It said instead of $\arg( Im(x) / Re(x) )$ we can use just $Im(x)$ because they are of the same sign because phase (here $x$) is always between $-\pi$ and $+\pi$ !!!

It is just a convention, but it is useful in some cases. For example, the phase of the DFT of a real discrete-time signal is odd only if the angles are expressed in the range $[-\pi, \pi)$. Sometimes you just have to adapt to the convention used by your tools -- for example, MATLAB functions like angle and atan2 return angles in $[-\pi, \pi)$.
Note that the frequency of a discrete-time signal, measured in radians per sample, is in the range $(-\pi, \pi]$, but in this case it has a different interpretation.


Replacing "e" in Euler's formula with another number Does Euler's formula remain valid if we use any real number other than the constant $e$?   For example replacing $e$ with 5 would make the formula look like this:  $5^{it}$.

I tried this idea in Matlab and replaced $e$ with few other real numbers (eg 1.5, 10, 2.1) and each time the plot was still showing what seemed like cosine and sine waves. The frequency of the cos and sin was changing depending on the base.

Here's roughly my approach:

w = freq * 2 * pi;
t = 0:0.001:1000 ;

a = real( number ^ (i*wt) ) ; % cos in Euler's formula
b = imag( number ^ (i*wt) ) ; % sin in Euler's formula

Say you're interested in $$M^{j2\pi f_0 t}. \tag{1}$$ Note that $$M = e^{\log M},$$ so $(1)$ can be written as

\begin{align}
M^{j2\pi f_0 t} &= \left( e^{\log M} \right) ^ {j2\pi f_0 t} \\
&= e^{j2\pi (f_0\log M) t} \\
&= \cos(2\pi (f_0\log M) t) + j \sin(2\pi (f_0\log M) t),
\end{align}
which is a complex sinusoid with frequency $f_0 \log M$. That is why using $M$ instead of $e$ results in a change in frequency.


Calculate Q factor of a Low Shelf and High Shelf filter I'm developing an app that has an EQ using eq10q plugin. My eq has 10 bands. I'm using peaking filters for the middle bands and a low shelf filter for the first band and a high shelf filter for the last band. I calculated the Q factor for the middle bands by dividing the middle frequency by the bandwidth desired for that band but I don't know how to calculate the Q factor for the High Shelf filter and the Low Shelf Filter. How should I do it?

it looks like this eq10q uses cascaded parametric EQs, initially set up to be bell, but switchable to other types (like shelving).

i wouldn't doubt that they're using the audio EQ cookbook.  if you're using that, there is a parameter called "shelf slope" or $S$ that is often set to 1, because that gets you the steepest slope without dips or lips or bumps.  so it's a straight incline or decline to the shelf level.  if the control is $Q$, instead, there is a relationship between $S$ and $Q$:

$$ \frac{1}{Q} = \sqrt{ \left( A + \frac{1}{A} \right) \left( \frac{1}{S} - 1 \right) + 2} $$

where

$$ A = 10^{dB_\text{gain}/40} $$

and $ dB_\text{gain} $ is the shelf gain in dB.

i would just plug in $S=1$ and get a $Q=\frac{1}{\sqrt{2}}$.


Sub-band partitioning of filter-banks In ECE-700 Filterbank Notes the partitioning is $ \Delta = \frac{2\pi}{M}$ with center-frequencies at $ \omega_k = \frac{2\pi}{M}k$ which makes me wonder, considering that frequencies in the range from $\pi$ to $2\pi$ are just the negative versions of the range from $0$ to $\pi$, which divide into  low- and high-pass frequencies.

In the follow-up picture the partitioning then stretches from $0$ to $\pi$ - is there someone so nice as to explain to me what I'm not understanding here?

The notes made me wonder as well, and I cannot suspect the author made a mistake. The first graph with periodicity is generic. The second one corresponds to dyadic iteration (wavelets), often displayed in $[0,\pi]$, except when dealing with complex wavelets where  $[-\pi,\pi]$ is used. I agree, $[0,2\pi]$ could have been used instead, but I have never seen the $2^k$ partitioning mirrored at the end.

You can see from the rest of the paper that the author uses $[0,\pi]$.


How Come RIP Guarantees Unique Restoration of the Sparse Solution by $ {L}_{1} $ Minimization? We have a sensing matrix $\Phi$, satisfying the restricted isometry property (RIP), and a sparse signal $x$. We want to recover $\hat x$ from the measurement $y=\Phi x$ by using $l_1$-minimization.

I am wondering about, how the RIP can guarantee the uniqueness of the recovered sparse vector solution $\hat x$.

On slide 4, at 01:15 minute, Mark Davenport says in his presentation about compressed sensing:

  The RIP is what we need to ensure, that this hyperplane is oriented in a way, that guarantees us, that there is actually only one sparse vector in this hyperplane. So we are going to be guaranteed to get the right answer.

But why is that? I do not see that the RIP can guarantee this.

For example, consider this picture below. It's a slide from the linked presentation. The blue hyperplane also intersects at least one other axis. How can we be sure, that this other sparse vector isn't the actual solution, i.e. was our formler signal $x$?

Compressed sensing theory (see for example Candès & Wakin, 2008) states that, given you have enough measurements $y$, the shown $\ell_1$ minimisation problem recovers $x$ exactly with very high probability. This means that your correct solution must have the smallest possible $\ell_1$ norm.

This is what the figure shows; the "$\ell_1$ ball" intersects the plane given by $\Phi x' = y$ at the point $x'$ with the smallest $\ell_1$ norm. There are, as you mention, other sparse solutions - i.e. the ones where the plane intersects the other two axes. As we can imagine from the figure, these points however have larger $\ell_1$ norm and are therefore not the correct solution.


Advice on designing a digital filter that doesn't have phase-sensitive edge artifacts? I'm pretty well versed in statistics, but not really digital signal filtering. I have a data scenario where I expected to be able to pretty easily filter out some noise (human pulse) that's at a known frequency band, but I'm having a lot of trouble using the standard tools in the scipy.signal library and think I must be misunderstanding how to design digital filters. I have a notebook here that walks through my explorations thus far, but the gist is that the standard scipy filters seem to cause large distortions at the start and end of my signal, with the precise behaviour dependent on the phase of the noise signal I'm trying to subtract. Just in case the above binder link goes down, I'll include some of the key points below as well:
First generating some synthetic data that's similar to my real data:
#generate time vector
samples_per_sec = 10.0
total_time = 40.0
time = np.linspace(0, total_time, int(total_time*samples_per_sec))

#generate the pulse signal
pulse_hz = 1.0
pulse_phase = np.radians(0)
pulse = np.sin(time*(2*np.pi)*pulse_hz - pulse_phase)

#generate the BOLD signal (just something that goes up then down)
dist = stats.beta(2, 2)
bold = dist.pdf((time-10)/20) / 10.0 # division by 10 to make bold a small signal

#combine
pulse_plus_bold = pulse+bold
plt.plot(time, pulse_plus_bold);

Try a 1st order butterworth:
#1st order butterworth filter in ba mode
ba1 = signal.butter(
    output = 'ba'
    , N = 1 #needs to be low if using output='ba', else use output='sos' and sosfiltfilt
    , Wn = [0.5,1.5]
    , btype = 'bandstop'
    , fs = samples_per_sec
)
filtered_ba1_nopad = signal.filtfilt(
    b = ba1[0]
    , a = ba1[1]
    , x = pulse_plus_bold
    , padtype = None
)
plt.plot(time, filtered_ba1_nopad, 'b');
plt.plot(time, bold, 'r--');
plt.legend(['Filtered', 'Expected'], loc=(1.04,.5));

First-order butterworth with even padding:
filtered_ba1_pad_even = signal.filtfilt(
    b = ba1[0]
    , a = ba1[1]
    , x = pulse_plus_bold
    , method = 'pad'
    , padtype = 'even'
)
plt.plot(time, filtered_ba1_pad_even, 'b');
plt.plot(time, bold, 'r--');
plt.legend(['Filtered', 'Expected'], loc=(1.04,.5));

First-order butterworth with odd padding:
filtered_ba1_pad_odd = signal.filtfilt(
    b = ba1[0]
    , a = ba1[1]
    , x = pulse_plus_bold
    , method = 'pad'
    , padtype = 'odd'
)
plt.plot(time, filtered_ba1_pad_odd, 'b');
plt.plot(time, bold, 'r--');
plt.legend(['Filtered', 'Expected'], loc=(1.04,.5));

This latter looks really good! But after playing around I discovered that whether odd or even (or either) padding works better seems to be contingent on the phase of the signal being filtered-out. As an example, while the above obtained excellent filtering with odd-padding, here's the same scenario but with a phase-shift added to the pulse signal that yields edge artifacts in both odd and even:
phase = np.radians(45)
pulse_shifted = np.sin(time*(2*np.pi)*pulse_hz - phase)
pulse_shifted_plus_bold = pulse_shifted+bold

filtered_shifted_ba1_pad_odd = signal.filtfilt(
    b = ba1[0]
    , a = ba1[1]
    , x = pulse_shifted_plus_bold
    , method = 'pad'
    , padtype = 'odd'
)
filtered_shifted_ba1_pad_even = signal.filtfilt(
    b = ba1[0]
    , a = ba1[1]
    , x = pulse_shifted_plus_bold
    , method = 'pad'
    , padtype = 'even'
)

fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(5, 3))
axes[0].plot(time, filtered_shifted_ba1_pad_odd, 'b')
axes[0].plot(time, bold, 'r--')
axes[1].plot(time, filtered_shifted_ba1_pad_even, 'b')
axes[1].plot(time, bold, 'r--')
fig.tight_layout()
plt.title('Odd (left) and Even (right)')
plt.legend(['Filtered', 'Expected'], loc=(1.04,.5));

I've also tried the 'gust' padding method as well as higher order filters (using sos of course), and observe the same phase-dependent edge artifacts in everything I've tried. Any tips?

Your basic problem is that filtfilt (and most other linear filtering routines) take filters that are designed for infinitely long time expanses, and apply them to a chunk of data as if the data were extended infinitely in both directions with zeros.
So you have a legitimate bandpass filter, and it's "seeing" a legitimate jump in the signal at the starting point of your signal.
There's three basic approaches you can take; the first two are ad-hoc and easy, the third one is difficult if you're starting from first principles.  It's certainly been solved out there someplace, but a brief search here on "filter finite-length data" didn't find me joy.
Approach 1: window the input data
Take your input data, and multiply it by something that'll make it taper off at the ends.  I.e. a ramp from 0 to 1 over 10 samples at each end, or $\frac{1}{2}\left (1 - \cos \frac{\pi n}{N} \right)$ for N samples at each end (suitably reversed on the left end).  You'll have some artifacts (a rising sine wave isn't the same as a steady one, after all), but they'll be attenuated.
Here's python code implementing a cosine edge attenuation with the ability to customize what central % of the signal is kept as 1:
def attenuate_edges(signal,time,edge_attenuation_percent):
  start = int(np.floor(len(time)*edge_attenuation_percent))
  end = int(len(time)-start)
  ramp = (1-np.cos(np.pi*(np.arange(start)/start)))/2
  edge_attenuator = np.ones(len(time))
  edge_attenuator[0:start] = ramp
  edge_attenuator[end:len(time)] = np.flip(ramp)
  return(signal*edge_attenuator)

Approach 2: Trim the output data
Do what you're doing now, and lop off the nastiness at the ends.  This is probably the easiest, and if you can just collect a bit more data, doesn't lose you anything.
Approach 3: Do a proper estimate of the interfering signal, and subtract it out
This will be fun-fun if you love math and have the time.  Basically you'll use the fact that the value of your interfering signal at time $n$ correlates in a specific way with the values of your interfering signal at time $k$ for all values of $n$ and $k$ in your data set.  You'll probably end up with something that looks a lot like a Wiener or a Kalman filter, that takes the end-effects into account.  Your estimate will be worse at the ends, but this will show up as a bit of noise on the ends -- not as those honkin' big pulses.
If I couldn't figure out the search terms for this, it would take me a day to do and another day to verify, and supposedly an expert.  OTOH, Gauss or Laplace probably invented it in the 19th century, and may even have thought it important enough to write down, somewhere.  So I'm sure the method exists.


Upsampling and pulse duration How does upsampling factor relate to pulse duration?
If I want to do pulse shaping and I have rcosdesign(beta,span=8,sps=4) should I upsample by 8, 4 or 33(8*4+1) and why?
I've seen you'd normally upsample by 4 and I can't see why.
Like, that way you're not using the whole pulse to shape a symbol, shouldn't it be 33?

I'd appreciate very much someone explained this to me. Thanks!

If the number of samples per symbol (sps) is 4, then you should upsample by 4.  The input to the pulse shaping filter is 1 sample/symbol, and the output is 4 samples per symbol, so the upsampling rate is 4.

The span is the number of symbol periods in the filter's impulse response.  Intuitively, this means that a particular symbol influences the pulse-shaped output for 8 symbol periods, or that the pulse-shaped signal at any moment in time contains a linear combination of 8 different symbols.  The number of taps in the filter is span*sps+1, which is 33 in this case.  The delay of the filter is span/2.  For a fixed filter span in symbol periods, the number of taps in the filter goes up with the number of samples per symbol because the output sample rate increases.  However, the span of the filter in symbol periods and the upsampling factor are two independent variables.  Changing the span does not change the upsampling rate.  Changing the number of samples per symbol in the pulse shaping filter is equivalent to changing the upsampling rate.


Standard deviation is varying with change in amplitude of DC signal I am trying to make sense of the below code. I am taking a large sample so that I get good estimates for standard deviation.

N = 25;
h1 = ones(1, N)./N;
L = 1000000;
s = ones(1, L);
n = 0.1*randn(1, L);
x = s + n;
stdx = std(x);
y = filter(h1,1,x);
stdy = std(y);

This code will give me stdx = 0.1 and stdy = 0.02, which is expected since $\sigma_y=\sigma_x/{\sqrt{N}}$as we are essentially doing an average filter.

But giving more amplitude to the signal  s = 10*ones(1, L); and running the same code above is giving me a stdy = 0.0344. Shouldn't the standard deviation be same in both the cases since the signal s has no variance?

I believe your intuition is correct.

The filter you're using is a low pass filter with a DC gain of 1. Without seeing the output of the filtering operation, I suspect that you're seeing the effect of initial conditions in the filter: getting from 0 to 1 in your first instance takes fewer samples than getting from 0 to 10 in your second.

To test this, try taking the middle $L/2$ of your filtered signal and see what std(y(L/4:3*L/4) is then.


Why Does the Rudin Osher Fatemi (ROF) Method Use Variational Methods for Image Denoising When Denoising Problems Are Not Boundary Value Problems? I was reading the recent book Variational Methods in Image Processing by Vese and Guyader which is quite interesting. In the book the authors discuss different types of image processing problems, from denoising and deblurring problems to image segmentation and beyond. The common theme is that they use Variational Methods, read as Calculus of Variations methods to take the noisy image and remove the noise. My question is, why do these scholars use Variational Methods (Euler-Lagrange Equation) to solve a problem that is not a Boundary Value Problem. Let me lay some foundation to the question.
Now these image problems are naturally optimization problems where I have some original image $f$ and I want to find a new image $u$ that best replicates $f$ but is subject to some constraints. The Rudin-Osher-Fatemi method sets up the problem this way:
$$
F(u)=\int_{\Omega}|D u|+\lambda \int_{\Omega}(K u-f)^{2} d x
$$
Where the second term is the standard difference between the original image $f$ and the adjusted image $Ku$, and the first term is the Total Variation norm in the image, meaning the sum of the absolute value of the first difference between adjacent pixels.
The solution to this problem sets up the Euler-Lagrange equation for the above energy function, and then minimizes that energy using gradient descent or such algorithms.
BUT, in my mind Calculus of Variations was always taught as a way to solve a boundary value problem. I am given two points, and I want to find the least energy path between those two points. So I have a BVP and I can either solve the optimization problem directly by discretizing the space between the two points and using shooting methods + opimizer to find the minimal energy path. Or I would use the Euler-Lagrange equation to find the least energy path as the numerical solution ODE or PDE--and thus avoiding minimization.
In looking at image processing problems, I did not have a sense of a boundary value problem? The initial image is fine but the final image is the result of optimization and not known before hand. Hence, can anyone explain why is the Variational approach used, or Euler-Lagrange approach?
Are they just using the E-L to find a minimum without thinking of a boundary value, or is there some other intuition going on that I am just not catching?

This turned out to be easier than I thought. Image processing is a boundary value problem, and the boundary is the set of pixels along the edge of the image. The common notation for this boundary in the Vese Guyader book is $\partial(\Omega)$, so that is not the most obvious notation for an image. Of course $\partial$ is often the boundary operator in differential topology too. I think for the Rudin-Osher-Fatemi method, there is a Neumann boundary condition, meaning that the first derivative at the boundary pixels is not changing.
So the answer was very easy after all.


Transfer function of second order notch filter Can anyone mention the transfer function of second order notch filter to remove the line frequency of 50 Hz, in terms of frequency and sampling rate.

Just like for Low pass Butterworth filter as,

$$
H= \frac{1}{\sqrt{1+\left(\frac{\omega_n}{\omega_c}\right)^4}},
$$

where $\omega_n$ is the signal frequency and $\omega_c$ the cutoff frequency.

For digital notch filters, I like to use the following form for a notch filter at DC ( $ \omega $=0):
$$ H(z) = \frac{1+a}{2}\frac{(z-1)}{(z-a)} $$
where $a$ is a real positive number < 1. The closer $a$ is to 1, the tighter the notch (and the more digital precision needed to implement).
This is of the form with a zero = 1, and a pole = $a$, where $a$ is real. The multiplication by $\frac{1+a}{2}$ is just to normalize the magnitude back to 1. To move this to a frequency, rotate the pole and zero to the frequency desired. For a real filter we end up with complex conjugate pole zero pairs, resulting in a 2nd order filter:
Defining a digital frequency range of 0 to 2$\pi$, with the sampling frequency at $f_s=2\pi$ and the notch frequency is $\omega_n$, then if we rotated the pole and zero above to $\omega_n$ we would get:
$$ H(z) = \frac{1+a}{2}\frac{(z-e^{+j\omega_n})(z-e^{-j\omega_n})}{(z-ae^{+j\omega_n})(z-ae^{-j\omega_n})} $$
Multiplying this out results in:
$$ H(z) = \frac{1+a}{2}\frac{z^2-2z\cos\omega_n+1}{(z^2-2az\cos\omega_n+a^2)} $$
So for your case of 50Hz, if we assume a sampling frequency of 1KHz, $\omega_n$ would be:
$$ \omega_n =\frac{f_c}{f_s}= \frac{50}{1000}2\pi$$
The coefficient $a$ is chosen to balance precision needed and bandwidth (bandwidth is tighter as $a$ approaches 1), and $cos(\omega_n)$ is a value between +1 and -1 that sets the frequency of the notch (+1 corresponds to DC with $\omega_n=0$ and -1 corresponds to $F_s/2$ with $\omega_n=\pi$, and any values in between for all frequencies in the first Nyquist zone.)
One possible implementation (using transposed Direct Form II) for this transfer function is shown below.

For example, see below digital notch filter with $a = .99$ and $\omega_n$ = 0.707. (Frequency axis is normalized where 1 = $f_s/2$


What is the intuition of "averaging is a low pass filter"? As I know, the shape of a low pass filter in time and frequency are as follow:

But how averaging work like a normal function in time domain?
Is there any intuition in it?

Consider a sequence of numbers

x = [1, 100, 1, 100, 1, 100]

If we send this sequence through an averaging transformation

y(i) = 0.6*x(i) + 0.4*x(i-1)

we get

y = [1, 60.4, 40.6, 60.4, 40.6, 60.4]

Original signal x had a lot of sudden changes 1 -> 100, 100 -> 1 where as signal y has smaller changes 60 -> 40, 40 -> 60

When we say that a signal has high frequency components we mean that the values change rapidly with time. So x had rapid changes in amplitude, while y does not have that much of rapid changes in values. This is the intuition behind why averaging is equivalent to low-pass filtering (disallowing high frequencies).

--

Another way to look at it

x = [1, 100, 1, 100, 1, 100]

could be written as sum of two sequences, x1 and x2

Low(zero?) frequency component

x1 = [50, 50, 50, 50, 50, 50]

High frequency component

x2 = [-49, 50, -49, 50, -49, 50]

After averaging transformation, we get

y = [1, 60.4, 40.6, 60.4, 40.6, 60.4]

Which has low frequency component of

y1 = [50, 50, 50, 50, 50, 50]

and high frequency component

y2 = [-49, 10.4, -9.4, 10.4, -9.4, 10.4]

As you can see that y2 has much smaller absolute amplitudes than x2, meaning that we have diminished the high frequency components in x by doing averaging.


Using Matlab/Simulink to build Reader based on Software Define Radio I am a student. I want to learn about Software Defined Radio to build a Reader for Saw tag and simulate them on Matlab/Simulink.

These bocks below is what I want to build in Matlab/Simulink

I found I could simulate RF front-end block using AD9361 because they are supported by Analog Device.

And I want to ask, can I build the digital blocks of Tx/Rx side and connect it to Tx/Rx AD9361 model to run the loop-back test like this (connect Rx to Tx) ?

You probably can. But: the AD9361 is definitely a complex beast, and I honestly think you're setting yourself up for more trouble than your solving when using a model of a device with a couple thousand different settings.

Unless you're planning to build an SDR out of an actual AD9361, I'd recommend not using them in a simulation. It's probably among the most complex models available in Simulink!

You shouldn't be building your own SDR, but rather buy one; as far as I remember, these typically need relatively high bandwidth, and designing the digital side of a PCB that carries multiple Megasamples of complex data per second into your computing unit will not be a beginner's project (aside from being far more expensive to buy and manufacture in single digit quantities than just buying a commercially available SDR).

The question really is whether you need to have a simulation that incorporates the actual SDR hardware – wouldn't for a very long time a simple idealized (IQ) mixer be representative enough? Aren't you starting off at the wrong end by trying to model something that works quite close to ideal very intimately instead of trying to get a system understanding based on idealized models?

For example, your flow graph top row isn't correct, to be honest. You used "QAM modulation" for what should be a quadrature modulator (QAM modulation is a baseband symbol mapping method, not a "modulation" as in "shift of frequency").


Deconvolution Question on Article "Deriving Intrinsic Images from Image Sequences" by Yair Weiss there are n derivative filters: $f_i$, and denote $f_i^r$ as $f_i$'s reverse filter such that

$$f_i(x,y)=f_i^r(-x, -y)$$

$r_i, f_i$ given, to find $r$ from the equations:
$$f_i * r = r_i, (1 \leq i \leq n)$$

Professor Weiss recovers $r$ using the following method:
$$r = g * \left(\sum_{i=1}^n f_i^r*r_i\right)$$ and $g$ satisfies:
$$\delta = g * \left(\sum_{i=1}^n f_i^r*f_i\right)$$

Why can the method solve the equation? What books should I refer to?
I want to know how to find the $g$?

Professor Weiss uses these two derivative filters: $[0\ 1 -1]$ and $[0; 1; -1]$.
and in his code the function invDel() returns $g$, which I could not understand. Article link, implementation link.

Professor Weiss's MATLAB code for getting $g$:

function [invK]=invDel2(isize)    % isize is 2 * max(imgWidth, imgHeight)
    K=zeros(isize);
    K(isize/2,isize/2)=-4;
    K(isize/2+1,isize/2)=1;
    K(isize/2,isize/2+1)=1;
    K(isize/2-1,isize/2)=1;
    K(isize/2,isize/2-1)=1;

    Khat=fft2(K);
    I=find(Khat==0);
    Khat(I)=1;
    invKhat=1./Khat;
    invKhat(I)=0;
    invK=ifft2(invKhat);
    invK=-real(invK);
    invK=conv2(invK,[1 0 0;0 0 0;0 0 0],'same');% shift by one

I didn't have time to fully read the article yet I'll try to answer according to your question.

If you take the first equation and convolove it on the left by $ {f}_{i} $ and use the second identity you can see this specific $ g $ holds the equation.

Regarding how to get this $ g $, I think he uses "Deconvolution" in the Fourier domain.


Proof of fourier transformation of multiplication of two signals I've been trying to find a proof of the following, but still I m unable to proof it, can someone help me?
$$
ℱ[x(t)g(t)] = \frac{1}{2\pi}
 [X(\omega)*G(\omega)]
$$

I am very thankful to Dilip-sarwate and Gilles, who took their precious time to understand my problem and guide me.
So, Now I'm going to write the correct solution to my question. Which is as follows:
 $$
ℱ[x(t)g(t)] = \int_{t=-\infty}^{\infty}[x(t)g(t)]e^{-j\omega t}dt
$$
As we know :
$$
x(t) = \frac 1{2\pi}\int_{\alpha=-\infty}^{\infty}X(\alpha)e^{j\alpha t}d\alpha
$$
Therefore,
$$
ℱ[x(t)g(t)] = \int_{t=-\infty}^{\infty}\bigg[\frac 1{2\pi}\int_{\alpha=-\infty}^{\infty}X(\alpha)e^{j\alpha t}d\alpha\bigg]g(t) e^{-j\omega t}dt
$$
After rearranging, we have
$$
ℱ[x(t)g(t)] = \frac 1{2\pi}\int_{\alpha=-\infty}^{\infty}X(\alpha)\bigg[\int_{t=-\infty}^{\infty}g(t)e^{-j(\omega-\alpha)t}dt\bigg]d\alpha
$$
As we also know that :
$$
G(\omega) = \int_{t=-\infty}^{\infty}g(t)e^{-j\omega t}dt
$$
and replacing $\omega$ with $\omega-\alpha$, we get:
$$
G(\omega-\alpha) = \int_{t=-\infty}^{\infty}g(t)e^{-j(\omega-\alpha)t}dt
$$
Therefore,
$$ℱ[x(t)g(t)] = \frac 1{2\pi}\int_{\alpha=-\infty}^{\infty}X(\alpha)G(\omega-\alpha)d\alpha$$
Finally, we get
$$
ℱ[x(t)g(t)] = \frac 1{2\pi}[X(\omega)*G(\omega)]
$$

[$\because$ $x_1(t)*x_2(t) = \int_{\tau=-\infty}^{\infty}x_1(\tau)x_2(t-\tau)d\tau$]


Trying to down-sample audio I am trying to downsample from $44100Hz$ to $8000Hz$, and I am not sure what I am doing wrong here.  $44100 / 8000$ is $5.5125$, so I am getting hung up a little on the fractional remainder which is causing my implementation to play back slightly slower than it is supposed to.  It seems that I need to somehow use the fractional remainder to compensate for this, but I am not sure how...

def downsample
  s = low_pass_filtered_signal_by_half_target_sample_rate
  samples_to_discard = (sample_rate.to_f / target_sample_rate).floor //   (44100.0 / 8000.0).floor

  s1_n = 0
  for (n = 0; n < N; n++) {
     if (n % samples_to_discard == 0) {
       s1[s1_n] = averaged(n, samples_to_discard)
       s1_n += 1
    }
  }
end

def averaged(start_index, samples_to_discard)
  sum = 0
    for (n = start_index, n < start_index + samples_to_discard; n++) {
      sum += s[n] || 0
    }
  sum / samples_to_discard
end

First you need to interpolate between samples rather than just retaining and discarding samples (which introduces horrible jitter noise during non-integer downsampling).

Then, if the original audio data contained spectrum around or above 4000 Hz (half the sample rate), you will need to low-pass filter before or in conjunction with the interpolation to the lower sample rate, or else you will get aliasing.

A good combined method of low-pass filtering plus interpolation is to use a windowed Sinc function as the interpolation kernel, where the Sinc is scaled to be the transform of an ideal "brick-wall" or rectangular filter with the appropriate low-pass cut-off frequency.  The width of the window on the scaled Sinc function will control the transition width of the frequency cut-off.


Getting wrong periodogram I'm driving a white noise signal $e(n)$ through an AR random process:

$$ x(n)=e(n)+\sum\limits_{k=1}^N \alpha_k x(n-k) $$

I got five values of $\alpha_k$ and plotted its power spectrum in Matlab for $\omega$ from $0$ to $\pi$.

Then I generated the realizations of $x(n)$ with 10000 samples

Now I'm trying to get the estimate of the spectrum from the expectation data using the periodogram formula

$$ S_{per}(\omega)= \frac{1}{K} \left|\sum\limits_{k=0}^{K-1}x[k]e^{-j\omega k} \right|^2 $$

with $K=256$ the number of first values of $x(n)$. When I run my code what I get is a vector with the same value in every position, so I'm plotting a straight line. I think that just applying the formula I should get the estimate, so what am I doing wrong? :-S

This is the section of the code I'm using to find the estimates

var_e=0.01;                                  % noise variance
b_1=1.4261;                                  % coefficients
b_2=-.7634;
b_3=-.9002;
b_4=1.2548;
b_5=-.5707;
K=256;
wn=rand(1,10000)-0.5;                        % white noise
samp=1:1:1e4;
T=10000;
x=zeros(1,length(samp));
% here I generate the expectations
for n=6:T+5
    x(n)=b_1*x(n-1)+b_2*x(n-2)+b_3*x(n-3)+b_4*x(n-4)+b_5*x(n-5)+wn(n-5);
end
x=x(6:end);
figure
plot(samp,x)
    for i=1:length(w)                         % w is the frequency, has 1000
        for k=1:K
        per2(i)=sum(x(k)*exp(-j*k*i));        % sum of exponentials
        pe2(i)=per2(i)*conj(per2(i));         % squared
        end
    end

Thanks in advance for your help

I think that is because your exponentials are wrong.
There is no variable named w in your code, although you refer to it, so I believe it should be something like:

per2(i)=sum(x(k)*exp(-2j*pi*k*w(i))));

But consider using fft instead:

pe2 = abs(fft(x)).^2;

or

pe2 = abs(fft(x(1:K))).^2;

if you insist on using only $K$ samples for some reason (why generate 10000 sample then??).

Perhaps you should consider the Welch method, with length 256

pe2 = pwelch(x,256,128,256);

BTW - you can generate AR process using filter:

x = filter(1,a,wn);

where a is the vector of AR coefficients.


What is the best method for measuring of room reverberation time? I want to measure a reverberation time in the room. What is a best way to do that? What kind of signal should be used? I know about balloon and pink noise, but are there any better methods?

Well, there is plenty of ways for approaching this problem. Which one is best depends on room itself, your resources and pursued accuracy. I think that best way to start is to look into ISO 3382-2:2008 standard, you will find lot of information. You should remember that post processing of your signals is very important here.

Nevertheless you should consider following methods:

Hand clap, balloon, starter gun. Cheap and most intuitive. You are recording the "impulse" itself from which you are then retrieving the reverberation time and other parameters. Basically almost every other method is relaying on getting the impulse response of your room. After that you filter it in frequency bands of interest and find the slope. Order is not random here. Hand clap doesn't have wide frequency response in low frequencies, balloon is slightly better, and starter gun will give best results.
Intermittent noise. You play the band-pass filtered pink noise (usually in octaves) and register it's level in dB scale. You turn off playback and record the decay. From the slope in dB scale you find the time at which energy curve will cross -60 dB from steady state. This is your reverberation time for this given octave band. It is important to use omnidirectional source and microphone. It is also good if you repeat results for different source-microphone arrangement, i.e. for 2 positions of your source and 6 positions of microphone.
Sweep sine. For that you presumably need a piece of software. You can either write your own (here is how to do it), use a free one (such us REW) or even buy some professional tool (i.e. Dirac, EASERA, etc.). Similarly to above you play signal from speaker, record for few positions and then analyse using favourite application - as simple as that. This is my method of choice.
Maximum Length Sequence. This is white-noise-like signal that correlates nicely with itself and can be used for extraction of impulse response. Personally I do not recommend it as sweep sine seems to work better for me in every case. Although some people like to use it for measurement of speaker cabinets.


What is the Compression Algorithm used by Facebook? I would like to know what is the Compression Algorithm used by Facebook for compressing the images? I uploaded a image of size 242 Kb and I downloaded the same image from my Facebook album. The size of downloaded image is 46 Kb and the subjective quality is similar, I wasn't able to differentiate between the original (uploaded) and downloaded (from Facebook) image. The downloaded image was downsampled though, from 1226 x 1600 to 552 x 720. So Facebook compresses the image at a good ratio while we are uploading it to our Facebook account without affecting the subjective quality much.

I wonder what compression algorithm they are using for this purpose. What would be your best guess?

                           Original image: tamil_girl_painting (size 1226 x 1600)

                  Downloaded (from Facebook) image: tamil_girl_painting (size 552 x 720)

From the images you posted its quite clear that the image has been downsampled and re-compressed with lower quality jpeg settings.  If you look round the mouth you can clearly see JPEG-like artefacts.


Capacity of AWGN channel I am confused understanding  basic concepts of communication over AWGN channels. I know the capacity of a discrete time AWGN channel is:
$$C=\frac{1}{2}\log_2\left(1+\frac{S}{N}\right)$$
and it is achieved when the input signal has Gaussian distribution.

But, what does it mean that the input signal is Gaussian? Does it mean that the amplitude of each symbol of a codeword must be taken from a Gaussian ensemble?
What is the difference between using a special codebook (in this case Gaussian) and modulating the signal with  M-ary signaling, say MPSK?

Assuming a channel whose input at each time is a continuous random variable $X$ and its output is $Y=X+Z$, where $Z\sim\mathcal{N}(0,N)$ and $Z$ is independent of $X$, then $$C_{\text{CI-AWGN}}=\frac{1}{2}\log_2\left(1+\frac{P}{N}\right)$$ is the capacity of the continuous-input channel under the power constraint
$$\mathsf{E}X^2\le P$$
The mutual information $I(X;Y)$ is maximized (and is equal to $C_{\text{CI-AWGN}}$) when $X\sim\mathcal{N}(0,P)$.

This means that if $X$ is a continuous Gaussian random variable with the given variance, then the output has the highest possible mutual information with the input. That's it!

When the input variable $X$ is discretized (quantized), a new formulation is required. Indeed, things can easily become difficult. To see it a little bit, one can consider the simple case of a very coarse discritization of $X$ where it can only have two values. So assume that $X$ is selected from a binary alphabet, for instance let $X\in\{\pm1\}$ (or a scaled version to satisfy a power constraint). In terms of modulation, it is identical to BPSK.

It turns out that the capacity (even in this simple case) has no closed form. I report from "Modern Coding Theory" by Richardson and Urbanke:

$$\begin{align}C&_{\text{BI-AWGN}}=1+\\&\frac{1}{\ln(2)}
\left(\left(\frac{2}{N}-1\right)\mathcal{Q}\left(\frac{1}{\sqrt{N}}\right)-\sqrt{\frac{2}{\pi N}}e^{-\frac{1}{2N}}+\sum_{i=1}^{\infty}\frac{(-1)^i}{i(i+1)}e^{\frac{2i(i+1)}{N}}\mathcal{Q}\left(\frac{1+2i}{\sqrt{N}}\right)\right)\end{align}$$
A comparison between the two cases can be seen in the figure below:


How Are Images Considered Non Stationary Signal When They Are Invariant to Time? I have read Wavelets are better than Fourier in dealing with non-stationary signals such as images, but I don't understand how images are considered stationary??

Stationarity is a statitical concept defined through stochastic processes. However the following is a practical account of why most typical images are nonstationary :

Throughout the image window, consider a small rectangular test block, say 16x16, and observe the frequency spectrum estimation of that block. Now move this test block accross the image and re-compute the resulting spectrum estimation each time. Now, if the spectrum does not change while the test block is moving accross the image, then you can loosely say that the image under concern is stationary.

Since for most natural images, such a moving block will be crossing unrelated distinct objects during its travesral of the image, the frequency spectrum will also be changing, i.e., image being nonstatinary

Basic Fourier type transforms are not very efficient and proper for this kind of local variations in spectra whereas wavelet or packets or short-space Fourier methods will be more efficiently indicating those local features.


Can I use more than one Window on an FFT? I was wondering if it is possible to take a signal and run it through an FFT multiple times with a range of different windows, and then combine/average the results to get a more accurate Signal analysis than just one window. For example, could i run my signal through a Hamming window to predict the component frequencies and then run a copy of the signal through a Flat Top window to get better amplitude predictions. combine the two data sets with some math and achieve a more accurate overall picture.
Are there any major flaws in this method? Are there better methods?

That makes little sense; the FFT (which is just an implementation of the DFT) is a linear operation; summing/weightedly averaging the same signal windowed with different windows is mathematically identical to just windowing with the summed/weightedly averaged window to begin with. It's literally the same.
So, what you really want is to design / choose a window that's appropriate for your use case instead!


Convert 16-bit wav file to 12-bit raw audio file I want to convert a 16-bit wav file to 12 bit raw audio samples (to simulate the output of an ADC).

I'm trying to run

$ sox file.wav --bits 12 --encoding signed-integer --endian little 12_bit_signed_le/file.raw

And I get the error

sox WARN formats: raw can't encode Signed Integer PCM to 12-bit

Output File    : '12_bit_signed_le/file.raw'
Channels       : 1
Sample Rate    : 4410
Precision      : 16-bit
Duration       : 00:59:53.47 = 15847219 samples ~ 269511 CDDA sectors
Sample Encoding: 16-bit Signed Integer PCM
Endian Type    : little
Reverse Nibbles: no
Reverse Bits   : no
Comment        : 'Processed by SoX'

sox INFO sox: effects chain: input         4410Hz  1 channels
sox INFO sox: effects chain: output        4410Hz  1 channels

As you can see, the command resulted in 16-bit output. How do I force it to truncate to 12 bits?

Also, do I need to be worried about clipping (as this post suggests?)

Sox will not support bit widths for wav files other than multiples of bytes. If you really need such unusual binary format then I suggest to write a short Python program that packs values in arbitrary format. Just please keep in mind that the maximum amplitude will be 2047 (instead of 32767 as for 16 bits).


Is match filter a solution of some differential equation? Matched filter is very common and termed as optimal filter for detection purpose. The question is : can we formulate a differential equation and show that the solution is matched filter.

As the name implies, a matched filter is found to be an exact replica of the signal of interest to be detected; just a reversed version as a result of an optimization.

Hence your question reduces to whether the signal of interest is produced as a solution of some diferential equation.

Even though it's true that there are signals as solutions to differential equations of all sorts (linear, nonlinear, time varying, time invariant etc), whether an arbitrary signal can be imposed to be a solution to some solvable differential equation is not very clear to me.

Since the signal of interest is an arbitrary signal, you must show that every (arbitrary) signal can be shown to be a solution of a corresponing differential equation. I would say no.

Note that as a degenerate case $y(t) = x(t)$ can also be considered to be a (useless) differential equation (of order zero?) but I think it's not what you mean by a differential equation.


Removing a darkened caption with text on it in an image Image processing  is new for me, so forgive me if my question sounds vague.

I'll try to supplement it with what I have learned until now and also a couple of images.

See the caption in the image below?

What I'm essentially trying to do, it to remove it and restore the original image in (Python using OpenCV).

Now I have a couple of approaches in mind. First one I read about is a technique called Inpainting. Now I saw a tutorial on inpaiting here but this required me to create a separate mask where the non-zero pixels denote the stuff I want gone.

Now what I noticed is that the caption is not fully opaque. So i was wondering if there is any possible way to restore the original image by first removing the darkened part of the strip. (essentially something very the original image with only the whitened text on it) and then create a mask of the text and then use inpaiting.

Now I have a couple of questions.

What technique do I use to remove the darkened part (let the text be now, we can remove it in the second step using inpaiting)
Does this algorithm even make sense. Is there a better approach I should be looking at?

NOTE: In no way am I looking for any sort of code or specific implementation. I'm just looking for what techniques and procedures I can study up on so as to get the job done. The rest is on me

Does this algorithm even make sense.

Your caption box is actually pretty wide – the algorithm that you've linked to would need to replace the complete content by its neighboring pixels. That won't look right – it will look as if someone "stretched" the guitar there.

So if you ask me, no, this is not the right approach here.

  Now what I noticed is that the caption is not fully opaque.

I was confused for a while until I realized that with "caption" you mean the box, and I mean the text. So I'll refer to the box to avoid further confusion.

So, yes, that would be what I'd look at here. The image just seems to be darkened in the box area – what about you just take this part of the image, and reverse that darkening, by increasing the brightness? I'd assume that degree of darkening is static – so a bit of trying out in e.g. GIMP you should be able to figure this out; I tried this, and  got this:

as you can see, it's slightly blueish- although the brightness seems to match. That might either be because whatever draws the box isn't actually reducing the brightness only, but somehow also changes the hue of the image in the region, but it might very well also be an artifact of JPEG compression; notice that your picture has a terrible quality to begin with. Or this is a color space conversion effect.


Detecting a frequency swept sinusoid and its parameters? Given a (FFT-sized) frame of data, and detection of a spectral component statistically above the noise floor in the FFT of this window, what characteristics or signal analysis could be used to determine that this spectral component is more likely to be a linearly swept sinusoid, rather than one that is stationary across the frame?

And, assuming the dF/dt sweep across the data window is small (from a fraction of an FFT bin to a couple bins), how can one estimate the sweep parameter (but, beyond answer offered to this question, assuming this is an estimation of a detected signal in noise.)

One offered solution seems to be to segment the FFT frame into several shorter subframes, do shorter STFTs, and look for a linear best fit among the resulting set of subframe FFT peak magnitude frequency estimates (which all have poorer frequency resolution due to the shorter subframes).  Are there any other or better options for detection and estimation?

You could compute the instantaneous frequency of the signal in the frame. This can be done as outlined in this paper (e.g. Eq.(9)). You need to compute the analytic signal using a Hilbert transformer:

$$x_a[n]=x[n]+jx_h[n]\tag{1}$$

where $x_h[n]$ is the Hilbert transform of $x[n]$. The instantaneous frequency is given by the (discrete-time) derivative of the phase of (1). If you use a first order difference for approximating the discrete-time derivative you'll get the formula given in Eq. (9) of the paper cited above. Given an estimate of the instantaneous frequency, you could simply fit a line to estimate the sweep parameter.

This little Matlab/Octave script shows how this approach would work with a toy problem:

% create chirp signal
N = 256;
n=(0:N-1)';
fx=(0.1+0.01*n/N);
x=sin(2*pi*fx.*n);

% compute analytic signal
xh=hilbert(x);
xr=real(xh);xi=imag(xh);

% compute instanteneous frequency (could be done over fewer points)
tmp1 = xr(1:N-1).*xi(2:N) - xr(2:N).*xi(1:N-1);
tmp2 = xr(1:N-1).*xr(2:N) + xi(1:N-1).*xi(2:N);
f=atan2(tmp1,tmp2)/(2*pi);  % (biased) estimate of instanteneous frequency

% linear regression over center values of f
skip = floor(N/4);    % skip unreliable points at beginning and end of f
na = skip:(N-skip);
A=[ones(length(na),1),na(:)];
u=A\f(na);

f0=u(1); k=u(2)/2;  % remove bias
[f0,k]

% reconstructed linear frequency
fr = f0 + k*n;

plot(n,fx,(2:N),f-k*(2:N)',n,fr,':')
legend('exact','estimated','reconstructed')
axis([0,N,f0*.9,(f0+k*N)*1.1])

The plot shows the actual linear frequency of the signal (blue), the estimate from the analytic signal (green), and the reconstructed frequency obtained from fitting a line through the estimated frequency. Note that the estimate is unreliable at both ends, so the line is fitted to the center values of the frequency estimate.


Metrics for Noisy Images I have a doubt. We calculate metrics such as PSNR, SNR, and MSE to prove that the proposed algorithm is better than its counterparts. To do this we need Original image (Noise free), Noisy image (we add this noise) and then Filtered image. Now we compare original and filtered output to estimate these parameters.

My question is, how to estimate these parameters when we have only the noisy image? That is in case of an ultrasound image, speckle noise occurs in the original image itself. We don't have a noisy free images in case of ultrasound. If we do as usual (as stated above), we will get poor results (PSNR,SNR & MSE) even though the subjective quality of the image is improved considerably.

Well, one way to do is to add extra noise to this noisy (original) image and apply our filter and show that the added noise has been removed efficiently. I want you to suggest any other ways to prove that my proposed algorithm is performing better than its counterparts.

Example:

              Original Image                       Filtered Image

How would you expect the effect of a denoising model? Do you want it to exceed your ability to observe the noise pattern and remove them, or to replace you and remove the known pattern which you define as noise?

I think most denoising algorithms target the 2nd option. That's the reason why it makes more sense to quantify and compare different algorithms through simulations. The paper assumes the main source of noise comes from speckle, which can be further assumed as Gaussian noise. So they added normal distribution values onto the 'clean' image, and PSNR, SNR, and MSE all get involved with the variation of the image. Such variation is just due to the adding values. They are actually trying to find out which algorithm can reduce such variation to the closest level of 'clean' image. So the comparisons have to be implemented with the known 'clean' image.

Strictly speaking, the 'clean' image in the paper is not clean at all. It is experimental data from ultrasonic B mode scan, and includes the unknown noise. I remember you asked a question in SO very long ago, claiming that you failed to reproduce some of the comparison results from the paper. That's very possible, since as I mentioned above, you can only evaluate your algorithms via simulations, and the 'clean' image should have been generated as a computational phantom. Otherwise you don't even know which part of the image from experiment, other than your added value,  is treated as noise in your algorithm.

As far as I know, SNR can be used as the quantification of experimental data, but it is not calculated that way in your paper. The aim of your algorithm is to enhance the region of interest (organ, tissues in the ultrasound image), making it clearer to the diagnosis. As a result, you'd better choose such regions by yourself, and select other background regions as the  'noise' part in SNR calculation. If you have more than one regions not viewed as background, contrast to noise ratio (CNR) is also suited in your contrast quantification (this is more used in MRI image analysis). Besides, edge spread function (ESF) along with modulation transfer function are often used in optical imaging.

You may also find some books on image analysis useful in image quality evaluation.


generating ultrasonic sound in matlab I am working for my class project on making animal repellent , I need to know how I could go about taking a normal audio file and shift its frequency/Pitch very high to make it an ultrasonic sound.

To answer the "how to shift the frequency of an audio signal up" bit:

You could multiply the signal by a sine wave at a high frequency. This would shift and mirror the whole spectrum of the original signal into the high frequencies (multiplication by a sine in the time domain = convolution by a pair of symmetric Dirac in the frequency domain) - the mirror image could then be removed by high-pass filtering. Limitation: this would not preserve integer ratios between harmonics. For example, a harmonic mixture with partials at 100, 200 and 300 Hz would have partials at 1100, 1200 and 1300 Hz if shifted by 1kHz - and the integer ratio between harmonics is lost, making the resulting audio signal "metallic" or "dissonant" - at least to human hearing's perception criteria!.
You could simply play back the original data at a higher sample rate. Limitation: multiplying the sample rate will also reduce the duration of the original data.
Actual musical pitch-shifting is done with techniques such as the phase vocoder (Short-Term Fourier Transform).

However, you probably don't need to bother with any of these methods because:

There is no need to shift an existing audio signal - given that the content of the original signal is of little relevance and will be degraded in one way or the other by the transposition process, why not directly generating a tone at the required frequency to start with?
Keep in mind that whatever you do, the converters and transducers that come with a computer cannot reproduce ultrasounds (they have a cutoff at the upper limit of the human hearing range, about 20kHz).
The evidence that ultrasounds have a repellent effect on mosquitoes is very shallow.


What is the meaning of **Sample** in **Sample ACF** What is the meaning of Sample in Sample ACF (autocorrelation function)

It's the empirical ACF computed using the sample.


Bandpass butterworth filter frequencies in scipy I'm designing a bandpass filter in scipy following the cookbook. However, if I decrecrease the filtering frequencies too much I end up with garbage at high order filters. What am I doing wrong?

from scipy.signal import butter, lfilter

def butter_bandpass(lowcut, highcut, fs, order=5):
    nyq = 0.5 * fs
    low = lowcut / nyq
    high = highcut / nyq
    b, a = butter(order, [low, high], btype='band')
    return b, a

if __name__ == "__main__":
    import numpy as np
    import matplotlib.pyplot as plt
    from scipy.signal import freqz
    # Sample rate and desired cutoff frequencies (in Hz).
    fs = 25
    # Plot the frequency response for a few different orders.
    plt.figure(1)
    plt.clf()
    for order in [1, 3, 5, 6, 9]:
        b, a = butter_bandpass(0.5, 4, fs, order=order)
        w, h = freqz(b, a, worN=2000)#np.logspace(-4, 3, 2000))
        plt.semilogx((fs * 0.5 / np.pi) * w, abs(h), label="order = %d" % order)
    plt.xlabel('Frequency (Hz)')
    plt.ylabel('Gain')
    plt.grid(True)
    plt.legend(loc='best')

    plt.figure(2)
    plt.clf()
    for order in [1, 3, 5, 6, 9]:
        b, a = butter_bandpass(0.05, 0.4, fs, order=order)
        w, h = freqz(b, a, worN=2000)#np.logspace(-4, 3, 2000))
        plt.semilogx((fs * 0.5 / np.pi) * w, abs(h), label="order = %d" % order)
    plt.xlabel('Frequency (Hz)')
    plt.ylabel('Gain')
    plt.grid(True)
    plt.legend(loc='best')

    plt.show()

You're probably running into numerical precision issues on the higher-order filters with sharp cutoffs. It could be a limitation of the butter function in SciPy, or it could be due to the filter structure that it's using when evaluating the frequency response. A better way of implementing high-order IIR filters is using a second-order sections (SOS) structure.

A little Googling suggests that SciPy doesn't yet support converting a filter to its SOS representation. Furthermore, there are apparently some issues with how filter parameters are handled internally that cause bad results with high-order filters. That sounds like what you're describing. These are apparently in the process of being addressed.

With that said, I came across this example that shows how you might implement an SOS structure using SciPy. I haven't tested the code at all, so there's no warranty, but it may be worth a look.


Image warping - how to proceed I'm studying image warping and I'm unsure how to warp a starting rectangular image into another image where the box corner points have been changed as I want:

(image from http://www.inf.ethz.ch/personal/pomarc/courses/CompPhoto/cpv03.pdf)

As far as I've understood, I need to:

1) Decide the four points (e.g. by using direct mapping on the border cases of the transformation)

2) Create a lattice of integer points on it (how?)

3) Do reverse mapping on those integer points to get the values in the original image

4) ??

I can't say I understood a lot of that, can somebody help me out please?

ImageMagick does that sort of thing, and it's open source.  Look at the source code in magick/distort.c of the current ImageMagick distribution, which you can get from imagemagick.org.  This assumes you have already selected the four control points and passed them in as an option.

The particular distortion method you described is called "binlinear forward" distortion in ImageMagick.  See http://www.imagemagick.org/Usage/distorts/#bilinear_forward