Does a simple photograph contain more information than a complex painting? I hope this question is appropriate for this site.

I came across this passage in The Three Body Problem, a novel by Liu Cixin:

  The professor had put up two pictures: One was the famous Song Dynasty painting Along the River During the Qingming Festival, full of fine, rich details; the other was a photograph of the sky on a sunny day, the deep blue expanse broken only by a wisp of a cloud ...
  The photograph's information content - its entropy- exceeded the painting's by one or two orders of magnitude

Representative pictures:

Is this true? How does one explain this counterintuitive phenomenon?

It depends how you define the term "information" or "entropy".

The conventional definition of entropy of an image is to think the image as a two-dimensional matrix of pixels and
$$H = - \sum_k p_k \log_2(p_k),$$
where $p_k$ is the probability, which is calculated from histogram, associated with gray level $k$.

This kind of entropy is correct if we ignore the correlation between pixels. For example the two images have the same entropy by this definition.

It is not true if the correlation between pixels is considered. For example if the color first pixel in top-left has probability $p_k$, the next pixel surely has the same color and its color does not have the same probability $p_k$.

We human being, with you as an example, use this kind of correlation to perceive the images. This kind of correlation are called "details" and, therefore, we/you expect that rich-detail images should have more information/entropy than simple ones. This is the reason why you have found it counter-intuitive.

PS:

I tried calculating the entropy of the two images you have posted, but they are not different "by one or two orders of magnitude" !!!!

"Along the River During the Qingming Festival" entropy about 7

"The sky" entropy about 6

They must not be the same files of the professor.


Phase shift problem in Fast Fourier Transform I try to make graph/print for magnitudes and phase shifts for impulse response calculated by FFT.
For magnitude everything works perfect, but for phase shift I get some strange curve for higher frequencies. I can't figure out why. Could anyone help me?

I calculate phase shift by that:
atan2(fftOutput.imag(),fftOutput.real()) * 180.0/M_PI;

And for simply impulse, like: impulse[1024] = { 1, 0, 0, 0 ... 0 }

with no processing I expect straight line (phase shift for all freq bin should be zero).

But I get something like that (I drawn it in paint, cause I can't run my app at the moment, but it looks almost exactly the same):

Why is that?

Something is wrong with your FFT. This looks like your input signal is either time reversed or shifted (circular) by one sample to the left.


What exactly is the effect of a Hann Window on the FFT output? In the code I am trying to understand, there is an usage of Hann window multiple times (There's overlapping frames, the size of each is 512 and the overlapping factor is 75% so each time we capture  128 samples we build a new frame by shifting the previous one by 128 to the left). For each Frame we apply a Hann window before the FFT, so the first 128 samples were actually "framed" 4 times.

However I don't understand the purpose of this windowing. Can someone please tell me what its effect is on the FFT ?

See the following: http://www.dtic.mil/dtic/tr/fulltext/u2/a381695.pdf  especially p.631 Figs. 17+18.  As noted by the author, a 4 to 1 overlap using a Hann or Hamming window allows for proper reconstruction.

Also see: F. J. Harris, “On the use of Merged, Overlapped and Windowed FFT’s to Generate Synthetic Time Series Data with a Specified Power Spectrum,” Proc. IEEE 16th Asilomar Conf., 1982, pp. 316-321.


Difference between cv::subtract and operator - in cv::Mat OpenCV What's the difference between cv::subtract and operator - in OpenCV?

Mat original = imread...;
Mat smoothed, difference, result;

smoothed - original image processed with Convolution with Gaussian kernel;

1)

difference = original - smoothed;

2)

cv::subtract(smoothed, original, difference);


Confusion regarding overshoot? As shown underlined in attached photo
Overshoot is mentioned as 9%
I am confused in understanding it
It is 9% of error or size of jump(1)?

It is the error with respect to size of jump. It has nothing to do with error at discontinuity. As per the explanation given, suppose the perfect square is having values $0$ and $1$. Just after the discontinuity of a rise time, the perfect square would take value $1$, while the reconstructed wave would take the value $1.09$ at the peak of overshoot. Here the size of discontinuity is difference between perfect low value $0$ and perfect high value $1$ of the square wave. In general if it were $A$ and $B$, the overshoot percentage would have been
$$
\frac{(Ovs_{peak} - (B-A))}{(B-A)} \times 100
$$

The error of $0.5$ in the error graph is not related to overshoot. For perfect square wave of $0$ and $1$, at discontinuity it is either $0$ or $1$, but the reconstructed wave will be midway $0.5$ ($|0.5−0|$ or $|0.5−1|$). So max error will be at discontinuity and its absolute value will be $0.5$.


Frequency estimation of circularly shifted single tone signal I have a discrete signal $y[n] = <e^{j ~ 2 \pi f ~ n}>_J + ~w[n]$ with $n \in [0, N[$ and $w[n]$ AWGN, $<x[n]>_K$ denotes the signal $x[n]$ circularly shifted by $K$ samples. Let's define $Y_k = \mathrm{DFT}(y[n])$ the output of the $N$-point DFT of $y[n]$. By ignoring the AWGN and using well-known Fourier properties, we can deduce that $Y_k$ is a sinc function centered around $f$, multiplied by $e^{\tfrac{k J}{N} }$ due to the circular shift.

In my application, $J$ is unknown to me but I am interested in obtaining an estimation of $f$ (which is real-valued in $[0, N[$). Over the years, many estimators $\hat{f}$ for single tone signals $z[n] = e^{j ~ 2 \pi f ~ n} + ~w[n]$ have been designed, and they work fairly well.
These estimators rely either on the complex values $Z_k$ (i.e. the DFT of $z[n]$, see e.g. [R1], which I can't use because I don't know $J$), or use the amplitudes of interpolated DFT bins, e.g. $\hat{f}_{res}$ which computes the residual frequency offset inside two DFT bins from [R2]:

$$ \hat{f}_{res} = \frac{|\bar{Z_1}|^2 - |\bar{Z}_{-1}|^2}{u(|\bar{Z_1}|^2 + |\bar{Z}_{-1}|^2) + v|\bar{Z_0}|^2} $$

where $\bar{Z}_k$ is the zero-padded $2N$-point DFT ($N$ points are from $z[n]$, the following $N$ points are $0$). $\bar{Z}_0$ is the bin with greatest amplitude and $\bar{Z}_1$/$\bar{Z}_{-1}$ are its neighbors. It is interesting to note that in the case of a $N$-point DFT, $|Y_k| = |Z_k|$. However this property does not hold anymore for the interpolated bins in $\bar{Y}_k$ obtained through the zero-padding, and I am unable to figure why analytically.

Here's an illustration of the difference of amplitudes between $\bar{Z}_k$ (zero-padded DFT of true single tone signal) and $\bar{Y}_k$ (zero-padded DFT of the circularly shifted signal):

and the MATLAB code used to generate the figure:

N = 256;
f = 100.4/N;
J = 125;

% Compute z[n]
n = 0:(N-1);
z = exp(2*pi*1j*f*n);

% Compute zero-padded DFT Z[k]
z2 = zeros(1, N*2);
z2(1:N) = z;
Z2 = fft(z2);

% Compute zero-padded DFT Y[k]
y2 = zeros(1, N*2);
y2(1:N) = circshift(z, J);
Y2 = fft(y2);

subplot(2, 1, 1);
stem(abs(Z2));
hold on;
stem(abs(Y2));
xlim([1 2*N])
legend('|Z(k)|', '|Y(k)|', 'Location','northwest')

subplot(2, 1, 2);
stem(abs(Z2) - abs(Y2))
xlim([1 2*N])
legend('|Z(k)| - |Y(k)|', 'Location','northwest')

Does someone have any clue on how to model analytically the interpolated DFT bins in $\bar{Y}_k$, or the DTFT of $y[n]$ ? Thanks.

[R1] https://www.researchgate.net/publication/3321864_Fast_Accurate_Frequency_Estimators_DSP_Tips_Tricks

[R2] https://ieeexplore.ieee.org/document/5910415 (unfortunately behind a paywall)

The clue is that without zero-padding, the circulary shifted sequence is just a shifted version of the periodic continuation of the original sequence. That's why without zero-padding the magnitudes of the DFTs of both sequences are identical.

When you use zero-padding, the two sequences cannot be obtained from each other by pure shifting. So the DFTs must be different.

You can avoid the problem by zero-padding the circularly shifted sequence between the last and the first element of the original sequence:

% this assumes that J is in [0,N]
y2 = [z(N-J+1:N),zeros(1,N),z(1:N-J)];
Y2 = fft(y2);

This will guarantee that the magnitudes of the DFTs of the two sequences are the same, even after zero-padding.


How to calculate the data rate of IEEE 802.15.4 to be 250kbps? How to calculate the data rate of IEEE 802.15.4 to be 250 kbps

You have $16$-ary symbols, each symbol requires $4$ bits. With a chip rate of $1000$ K, the symbol rate is $1000/16=62.5$ K, hence the bit rate is $4\times 62.5=250$ Kbps. Did you mean this calculation?


Why does reversing the order of these two transfer functions give me different outputs? Consider these two systems:
\begin{align}
&u\ {\longrightarrow}\boxed{s}{\longrightarrow}{\boxed{\frac 1s}}{\longrightarrow}\ y\\
&u\ {\longrightarrow}\boxed{\frac 1s}{\longrightarrow}{\boxed{s}}{\longrightarrow}\ y
\end{align}

If I input a step ($u=1/s$), I get an output of $y=0$ from the first system, but $y=1$ from the second. (You can see this from considering the signal in between the two blocks - it's zero in the first case and a ramp in the second case).

Why does reversing the order of the blocks change the output?

If the input is a unit step, then the output of the first block in system 1 is not zero, but it is a Dirac delta impulse $\delta(t)$. Intuitively, the derivative is infinite at $t=0$ because of the step going from zero to one. Integrating the Dirac delta impulse will give you a step at the output of the first system.

The output of the first block in system 2 is a ramp, which will be transformed into a step by the differentiator. So both systems have the same output.

Note that in general you can only interchange the order of LTI systems without changing the output if all systems are stable. For unstable systems, as in your example, the order can make a difference, because the corresponding integrals might not converge, and consequently you can't change the order of integration.


Perform overlap add with available FFT size smaller than filter coefficients I would like to process 12 to 20 seconds of incoming audio at a sample rate of 44100. I must process this audio in real time in an STM embedded kit (perhaps also an Android Smartphone). I'm trying to detect and count the number of occurrences of a signal of roughly 6500 samples inside the incoming audio. The maximum FFT available is of 1024 samples.

I was thinking about applying overlap-add but the number of coefficients would be 6500 and that's larger than maximum FFT size of 1024. I tried to simulate this in Matlab using fftfilt but the function help says:

If you supply a value for n, fftfilt chooses an FFT length, nfft, of 2^nextpow2(n)and a data block length of nfft - length(b) + 1. If n is less than length(b), fftfilt sets n to length(b).

This makes me think that I'm forced to use an FFT of at least 6500 samples (which I can't) and then process 1 incoming audio sample at a time (super inefficient).

What can I do?

You can make a big FFT out of smaller FFTs

This code implements a 16384 point FFT with a 16 point FFT and  1024 point FFT.

You need only calculate the Twiddle matrix once.

clear all
M=16;
N=1024;
x=sin(linspace(1,M*N,M*N)*2*pi*60/(M*N));  % test signal
X=reshape(x,N,M).'; % form 2D matrix read data in as rows
Twiddle=zeros(size(X));    % make Twiddle matrix
for i=1:M
for k=1:N
Twiddle(i,k)=exp(-1j*2*pi*(i-1)*(k-1)/(N*M));
end
end
X=fft(X); % fft on each column
X=X.*Twiddle;% element by element product
X=fft(X.').' ; %fft on each row
y=reshape(X,N*M,1); % read out as columns
figure(1)
plot(abs(y),'linewidth',2)
title('Composite DFT')
figure(2)
plot(abs(fft(x)),'linewidth',2)
title('Direct DFT')

which is based on a section in

Rabiner, Lawrence R., and Bernard Gold. "Theory and application of digital signal processing." Englewood Cliffs, NJ, Prentice-Hall, Inc.,
1975. 777 p. (1975).

I don't have the book handy for a page number, but the section is the table of contents under something like 1D DFT as a 2D FFT.

update

section 6.8 page 371
A unified approach to the fft


Why are window functions used when calculating spectral estimates, really? Having just read the rather long article on Wikipedia about window functions, I still do not see the mathematical point in using window functions when calculating spectral estimates (using e.g. Welch's method).

For example,

  The Fourier transform of the function cos ωt is zero, except at frequency ±ω. However, many other functions and waveforms do not have convenient closed form transforms

what is this supposed to mean? I doubt any sampled and digitized signal has a convenient closed form transform anymore.

You are right that no real-world signal will have a "convenient closed form" transform. I also find the quoted sentence misleading, and in my opinion it does not motivate the use of windowing. From your question I believe that you do not actually ask about the pros and cons of different window functions, but if I understand you correctly you ask yourself why people cut signals into pieces anyway if they want to analyze their spectrum. The reason is that many real-world signals are simply too long and/or their length is not known. So you don't want to wait till the signal is 'complete', and you also don't want the complexity of processing such a long signal. So you do block-processing, which reduces computational complexity and gives you quickly a spectral estimate. When you cut the signal into blocks you have many parameters that will influence the result of the spectral estimate: the block length, the type of window, the overlap between blocks, etc. This is a complex topic but you'll find a lot of literature about it. If you have access to Oppenheim and Schafer's Discrete-time Signal Processing I would recommend the chapter on Fourier Analysis of Signals Using the Discrete Fourier Transform. It has a lot of valuable information about block processing and windowing.


Physically, what does the usage of two variables mean for convolution My intuition of convolution is that it is just a way to depict multiplication of two signals where each signal is made up of various frequencies and phases.
Since it isn't easy to find the value of $\cos({w1}+{\phi1})+\cos({w2}+{\phi2})$ multiplied by $\cos({w3}+{\phi3})$, we need a simpler method to represent our resulting signal, similar to how we would use phasors when adding signals of the same frequency.
And intuitively, convolutions are seen as sliding windows.
But I am confused about how that logic gets translated into our equation with $t$ and ${\tau}$
$$
y(t) = x(t) * h(t) = \int_{-\infty}^{\infty}x(\tau)h(t - \tau)d\tau
$$
${\tau}$ is very clearly some dummy variable but beyond that I feel a little confused.
The value $t$ is used to used as time to check out value at a certain instance. And I think using t as a time-shift allows us shift the impulse function to just see it as a moving periodic function function and we are purely changing the origin point of our impulse response.
Because that matches what you are taught as an undergraduate and you draw different regions of overlap.
Is two variables then just to depict:

How much overlap there is between signals since there will be points of the signal that "crossed" the LTI system or haven't reached our filter yet.

What the status of signals are during that overlap?

Would my intuition so far be correct? Which variable is representing this term and how does it relate to the pen and paper method for looking at regions of overlap?

My intuition of convolution is that it is just a way to depict multiplication of two signals where each signal is made up of various frequencies and phases.

Not really, no.  I can't think of a case where you'd actually receive two signals and convolve them (correlate them -- yes, but not convolve).  This may be why you're not understanding the rest.
Convolution is used to find the output of a linear system* when you know the system's impulse response.

And intuitively, convolutions are seen as sliding windows.  ... $\tau$ is very clearly some dummy variable but beyond that I feel a little confused.

$\tau$ is the dummy variable that implements the sliding window in the convolution.  For each time $t$, you multiply the signal by the time-reversed impulse response, and you integrate that product.  $\tau$ is just the variable over which you do the integration.
* I purposely left the "time-invariant" part of that out -- convolution applies to time-varying systems, too, it's just more complicated.


Designing a multi-band FIR filter on MATLAB I am trying to plot the frequency response of a dual recursive running sum (RRS) filter, with the $\mathcal Z$-transform of each given by:
$$\frac{1 - z^{-L_n}}{1 - z^{-1}}$$
This is the code I've written:
% (1 - z^-L1) / (1 - z^-1) * (1- z^-L2) / (1 - z^-1)
L1 = 75;
L2 = 53;
fc = 500;
fs = 44100;

den_exp = -exp(-2j*pi*(fc/fs));
num_exp = -exp(-2j*pi*(fc/fs)*L1);
num2_exp = -exp(-2j*pi*(fc/fs)*L2);

%z transform
b1 = [1, zeros(1, L1-2), num_exp];
b2 = [1, zeros(1, L2-2), num2_exp];
a = [1, den_exp];

[hA, w] = freqz(b1, a);
[hB, w] = freqz(b2, a);
hAB = hA .* hB;

plot((w*fs)/(2*pi), 20*log10(abs(hAB)))
xlabel('Frequency')
ylabel('Magnitude (dB)')
title('Recursive filter with cut-off = 500Hz')

The frequency response appears to be off, and I can't seem to figure out why.
This is part of a project where I am trying to build a tone control system on MATLAB. I have read a couple of papers on the RRS, and I understand the idea of combining the transfer functions of high, low, and band pass with selected gains and delays, but I am struggling to implement this on MATLAB.

For example, how does one add a delay block to a transfer function, as shown here?
And how are these transfer functions added together?

I would be grateful if someone could provide some details on how to implement this.

According to your filter equations, your polynomials for $b1$, $b2$ and $a$ should only contain 1, zeros and -1. Your don't need to evaluate the last coefficient (num_exp, and others) at the complex exponential - That is done inside the freqz()
Take a look at this example here
So you should have:
a=[1, -1]
b1 = [1, 0, 0, ..., 0, -1]
b2 = [1, 0, 0, ..., 0, -1]

So $a$, $b1$, $b2$ contain the coefficients of the $z^{-n}$ terms.


FIR filter reverse engineering from step response I am working on a reverse engineering problem - I am trying to reconstruct a filter.

I have:

The step response of the system. (I could calculate the impulse response by derivating)

I know:

The system consists of a low-pass FIR and a decimation unit (decimation by 4, probably).
In the step response, there is always a point in the middle of the step amplitude (the attached image)
The group delay of the entire system is around 68.

The impulse response should be symmetrical - for all obtained impulse responses (I have acquired several step responses) there are two points with similar amplitude - apart from those two points the impulse response is symmetrical.
According to the manufacturer, the filter should have 5% skirt and 104 dB attenuation.

My main task is to repeat the behaviour of the system - convolving new signals with the obtained impulse response works - still - I would like to have a model of my filter.

I would like to synthesise a filter (FIR), which similar characteristics - with group delay around 68.  (the filter should have more than decimation*68*2 taps)

I cannot obtain a filter which would guarantee this middle point. Could you recommend me a tool/approach? (currently, I'm using the MATLAB filter design tool)

EDIT. I would like to generate several filters and select the best filter. (more Bayesian approach)

Thank you in advance

The relation of the discrete-time step response, $s[n]$, and discrete-time impulse response, $h[n]$, is

$$ s[n] = \sum\limits^{n}_{m=-\infty} h[m] $$

and

$$ h[n] = s[n] - s[n-1] $$

from $h[n]$ you can compute the DFT (of a large finite segment of $h[n]$) and get the frequency response of this filter:

$$ H(e^{j \omega}) = \sum\limits_{n=-N/2}^{N/2-1} h[n] e^{-j \omega n} $$

remember, with a finite delay allowed for the filter (this is your 68 samples), you can model both your impulse response $h[n]$ and step response $s[n]$ as symmetric about $n=0$ and define the origin for $n$ there at the 68th sample.  doesn't matter if that makes the filter acausal.  don't sweat that.

also, keep in mind that the operation that turns the step response into the impulse response is a discrete-time differentiator, which means the information about DC is lost.  the DC component in your impulse response that is derived from your step response is the step size of the step response, $s[+68]-s[-68]$.  you can add to both any constant and that will not change the impulse response $h[n]$ at all.

oh, and expect 2×68 + 1 taps.


How does BER performance depend on gray or natural binary? I read that BER performance would be improved if Gray code is used in 16 QAM rather than natural binary coding.  I think that it should not depend on the coding used.  16 QAM simply assigns to each of the 16 combinations of 4 bits, one particular sinusoidal signal with fixed frequency but different amplitude and phase.  When this is transmitted and detected properly at the receiver we can get the 4 bit group transmitted.  The error depends on the detection capability and the noise in the system.  Any arbitrary assignment between the signals and the group of 4 bits should give the same BER as long as the detection performance is the same.  Where am I going wrong?

You are confusing symbol errors and bit errors. What you say is true for the symbol error rate; it only depends on the signal-to-noise ratio (SNR). However, depending on the coding, i.e. the mapping of the bits to the symbols, one symbol error can cause one or more bit errors. The most likely symbol errors are the ones where a symbol is confused with one of its closest neighbors. Gray coding makes sure that those most likely symbol errors cause only one bit error. Other types of coding will result in more bit errors per symbol error. This is why Gray coding minimizes the bit error rate for a given symbol error rate (i.e., a given SNR), if one assumes that symbols are only confused with their closest neighbors, which is very likely for moderate to high SNRs.


Max frequency of cosine in SDR I would like to transmit a 56 MHz real cosine and receive it using SDR with max complex sample rate equals 61.44MSPS and an analog filter with BW = 56MHz (the transceiver is AD9361). Can I do it? How should I generate its samples? (I think I can't since the RF BW is 112 MHz (full BW) around the center frequency, but am I wrong?)
Are these Tx chain and Rx chains the correct way to implement this transmission?
TX

RX

Are these Tx chain and Rx chains the correct way to implement this transmission?

No. When you want to transmit a real cosine in passband, then you need a single complex tone in baseband. Not a real one. Otherwise, you'll have two real tones on the air, one at the mixing frequency + cosine frequency, and one at the mixing frequency - cosine frequency.
But maybe, you wanted that, I'm not sure from your description. However, the fact that you're calling things "chirp" makes this unlikely.
Then, a real-valued sampled signal at a sampling rate of 60 MHz can't represent 56 MHz. The highest frequency it can represent is less than 30 MHz. Nyquist (The Nyquist-Shannon sampling theorem, to be exact)!
Also, you can't throw away the imaginary part in your receiver: you have no knowledge of your phase. So, you're unsynchronized in phase, and then you can't assume that the real signal you've sent is received as real signal. Especially if this, as I might guess, is about channel estimation or radar, you can't make that assumption.
So, I'm sorry but: you really need to refresh your education on sampling theory, and on what equivalent complex baseband and bandpass signals are, first!

I think I can't since the RF BW is 112 MHz (full BW) around the center frequency, but am I wrong?

Yes, the bandwidth is 56 MHz, not twice that.

Also, delete the file sink on the transmitter site: that signal you're transmitting is deterministic, and it's easier to just recreate it than read it from a file; in fact, that file sink hurts a lot, since you need to generate samples at a pretty high rate under hard real-time constraints. You're trying to write 60 MS/s · 32 bit/S = 1.92 Gb/s to a storage device; that takes computational power and bus bandwidth away from your PC, and acts as an artificial slow down.
I've seen few SSDs that can sustain that rate for a significant amount of time under these constraints.
On the RX side, as said, you need to be complex, so you need to write 60 MS/s · 2·32 bit/S = 2·1.92 Gb/s = 3.84 Gb/s to a storage device; that really is challenging. You might want to start with Vector sink first, which stores the samples in RAM, and use the "reserve memory" parameter of that to preallocate RAM for these samples; use a head block to cut off after you've collected enough samples and read the contents of the vector sink in Python. (That "reserve memory option" might only available in recent versions of GNU Radio – you're using an outdated one)


Fair performance comparison betweem LMS & NLMS How can I choose the step size $\mu$, when I'm comparing different algorithms such as LMS, NLMS and transform domain adaptive filters, regarding their convergence speed, to get a fair comparison between them for mean-squared-error learning curves?

Is it necessary for a fair comparison to use the same step size for all algorithms since in most of the literature $\mu$ is not specified for the individual algorithms?

For fair comparison of one algorithm to another, the value of step size does not need to be same. You can adjust the step sizes of both algorithm so that the mean-squared-error learning curves base floor gets same and in this way you will be able to differentiate the performance of algorithm. For base floor I mean the value at which the mean squared error curves coincides with each other at the end of iterations or data samples. You can adjust the step size by hit and trial method.


Combining spectrograms with different windows to getting arbitrary time and frequency resolution Spectrograms convert an amplitude-time signal into a frequency-time signal using Short Time Fourier Transform (STFT). The window size L used in STFT is a design choice. As L increases, you get higher frequency resolution and lower time resolution and vice versa.

If you performed STFT with a large L to get spectrogram A with high frequency resolution then performed STFT with small L to get spectrogram B with high temporal resolution, can you combine spectrogram A and B somehow to get high temporal and frequency resolution?

My gut feeling is no, as changing L is simply "changing how you interpret the amplitude-time data."

Improving the time/frequency resolution balance remains a standard research issue, trying to adapt the STFT, inventing novel/better time-frequency distributions (Wigner, Cohen), using together or combining the above. For instance in music or speech source separation, a couple of STFTs is often used with both short and long windows, to be able to separate   transient from larger signals.

The ways you can combine such SFTF depend a lot on the exact processing or features you want to extract, as detailed by @hotpaw2. Here are a few pointers to references suggesting methods. The first one suggests a point-wise combination of two spectrograms as $$S_c(f,t) = (|S_1(f,t)||S_2(f,t)|)^{1/2},$$ and the second generalizes it to the geometric mean, the minimax or the reciprocal average.  The third encompasses it in a more generic merging approach.

Combined multi-resolution (wideband/narrowband) spectrogram, ICASSP 1991
On Combining Evolutionary Spectral Estimates, 1994
Adaptive Time-Frequency Resolution for Analysis and Processing of Audio, AES 2006


Autocorrelation function and correlation integral I am confused by the definition of autocorrelation function. It is originally defined as the expected value
$$R_{XX}(\tau) = E[(X(t)X(t+\tau)] = \langle X(t)X(t+\tau)\rangle\tag{1}$$
where $\langle\cdot\rangle$ is the ensemble average, and suppose they are all real signals. However, I see another definition defined by the correlation integral (e.g., see Autocorrelation - Wiki), which reads
$$R_{XX}(\tau) = \int_{-\infty}^\infty X(t)X(t+\tau)dt \tag{2}$$
Under assumption of ergodicity, Eq.(1) can be approximated by temporal averaging, i.e.,
$$R_{XX}(\tau) = \lim_{T\to\infty}\frac{1}{2T}\int_{-T}^TX(t)X(t+\tau)dt\tag{3}$$
But I don't see how Eq.(3) turning into Eq.(2), or how Eq.(2) being equivalent to Eq.(1).

So am I misunderstanding or missing anything here, or the one defined by correlation integral is another version of "autocorrelation"?

Any hints would be helpful!

Thanks!

EDIT: Another question is, for Eq.(2), if one sets $\tau=0$, the result is expected to be the variance (or "energy") of the signal. However, if we start from Eq.(2),
$$R_{XX}(0) = \lim_{T\to\infty}\frac{1}{2T}\int_{-T}^T X^2(t)dt$$
For any signal in $\mathcal{L}^2$ space (which has finite energy), the integral will be finite, but as $T\to\infty$, $R_{XX}(0)\to0$; that's not physical. Where did I go wrong?

The definition of auto-correlation depends on the type of signal. For random processes, the auto-correlation function is defined by the expectation given in Eq. $(1)$ of your question.

For deterministic signals, there are two definitions, depending on whether the signal is an energy signal (i.e., has finite energy), or a power signal (i.e., has finite power but infinite energy). In the first case, the auto-correlation is defined by Eq. $(2)$ in your question. In the latter case, the integral $(2)$ doesn't exist, and the auto-correlation is defined by the limit given in Eq. $(3)$.


Fast algorithm for n-dimensional DCT I need to implement an encoder which compresses a 5-dimensional structure of 10 bits values. Each dimension has between 4 and 12 elements. If a dimension ever has more than 12 elements, it is partitioned in half.
So far, I have been using a separable DCT to do that, but my implementation is very slow: I have implemented a naive DCT transform, which has complexity of $O(N^2)$.
I have tried some fast algorithms, but they seem to be slower than the naive approach:
a naive DCT for 6 elements is faster the Fast version, but for $N = 4$ and $8$, the fast approach is definitely faster. Also, I cannot force the friendly power-of-2 sizes.
I remember the Cooley–Tukey FFT algorithm decomposes a sequence of size $N$ into their factors $N_1$ and $N_2$ such that $N = N_1N_2$. I also know that there are specialized algorithms for 2D DCT due to their use in video and image codecs.
My question is: Is there any reference about a fast multidimensional DCT or even a "reversed" Cooley-Tukey
for DCT of which I could transform my 5D array into 1D and use a fast DCT?
Exrta bits:  We are researching not only the DCT-II but the whole family of trigonometric transforms., insights about any of the DST and DCT are welcome!
Edited 9/23
Just as an extra comment, this solution is being implemented in C++.

Since these are very short vector lengths, a "conventional" fast transform is not going to give you any gains.
Your best bet is probably to hand code each possible transform length indivdiually, unrolling all the loops and taking advantage of all "trivial " coefficients manually.
Another option could be to reduce the dimonsionality of the data but that depends how much codec gain you get in each dimension. If there is one dimension that has most of the codec gain, you can keep this and just linearize the other ones.


remove zero padding effect crosscorrealation I am trying to make the cross correlation on the frequency domain and I
use the fft function which takes power of 2 input data so I must pad with zero
my incoming sample signal and localy generated special code to take fft of
this signals.

I compare the first maximum and second maximum peak ratio and then I decide if the code was found or not found. When I use zero padding version this ratio and
phase change.

How can I remove zero padding unwanted effect from this correlation?

Thanks for helps.

so let's say the length of your FFT is $N$.  let's say that you fill half of the buffer with your signal and fill the other half with zeros.

$$ x[n] = 0 \quad \text{  for  } \tfrac{N}2 \le n < N  $$

the DFT of $x[n]$ is $X[k]$.  then magnitude-square and inverse DFT:

$$ R_x[n] = \mathcal{iDFT}\bigg\{ \bigg|X[k]\bigg|^2 \bigg\} $$

in the time domain, this comes out as

$$ R_x[\ell] = \sum\limits_{n=0}^{\tfrac{N}2 - 1 - \ell} x[n] x[n+\ell] $$

we call $\ell$ the "lag" of the autocorrelation.

now suppose $x[n]$ is periodic with some period $P$ which is much smaller than $\frac{N}2$.

$$ x[n+P] = x[n] \quad \text{for } 0 \le n < n+P < \tfrac{N}2  $$

with a lag of zero

$$ R_x[0] = \sum\limits_{n=0}^{\tfrac{N}2 - 1} (x[n])^2 $$

which, if we apply a little bit of statistical sleight-of-hand, we might expect is $\tfrac{N}2$ times the expected power $\overline{x^2}$.  in the limit with $P \ll \tfrac{N}2 $we can sorta do this with periodic $x[n]$.

$$\begin{align}
R_x[0] &= \sum\limits_{n=0}^{\tfrac{N}2 - 1} (x[n])^2 \\
       &\approx \sum\limits_{n=0}^{\tfrac{N}2 - 1} \overline{x^2} \\
       &= \frac{N}2 \overline{x^2}
\end{align} $$

now with a lag of $mP$ where $m>0$ is an integer:

$$\begin{align}
R_x[mP] &= \sum\limits_{n=0}^{\tfrac{N}2 - 1 - mP} x[n] x[n+mP] \\
       &= \sum\limits_{n=0}^{\tfrac{N}2 - 1 - mP} x[n] x[n] \\
       &= \sum\limits_{n=0}^{\tfrac{N}2 - 1 - mP} (x[n])^2 \\
       &\approx \sum\limits_{n=0}^{\tfrac{N}2-1 - mP} \overline{x^2} \\
       &= \left( \frac{N}2 - mP \right) \overline{x^2} \\
\end{align} $$

eliminating $\overline{x^2}$ we have

$$ R_x[mP] \approx \left( 1 - \frac{2mP}N \right) R_x[0]  $$

for $P \ll \tfrac{N}2$.

this suggests an envelope on the autocorrelation of

$$ \big| R_x[\ell] \big| \le \left( 1 - \frac{2 |\ell|}N \right) R_x[0]  $$

this is for the rectangular window.  it can be generalized for other windows.  the envelope results from the two windows having less area in common as the lag $\ell$ increases.


Zero Meaning a Hamming Window in Matlab I am still relatively new to signal processing and MATLAB and am working on a short-time fourier transform of some data using a prexisting code. Within the code a Hamming window is produced and multiplied by the signal prior to an fft.
The function itself works correctly, however I want the peak of the Hamming window to be over the first sample, rather than the midpoint (i.e shifting the Hamming window back by (no of samples/2)). What would be the simplest way to achieve this zero meaning?
This is the current generation and application of the Hamming window, where ncol is the number of time windows being used in the stft.
y = window(:,ones(1,ncol)).*y;

Any help is greatly appreciated.

Method1:
w = w([(end/2+1):end 1:end/2]);

Method#2:
w = circshift(w, ncol/2);


How to determine brown noise variance from samples? To be correct I want to determine variance of white noise before integration (brown noise is white noise after integration). For example, I have a following model (in MATLAB) of noises acting on my sensor:

How can I determine variance (power) of 'white noise 2' from samples (total noise)?

Of course noises are band-limited. I have to take into account following noise model of my sensor in my project but I got stuck on this.

\begin{align}
\tilde \omega(t)&=\omega(t)+\beta (t)+\eta_v (t)\\
\dot \beta(t)&=\eta_u(t)
\end{align}

Can I use power spectral density or Allan variance plot for this (how)?

If the white noise has a Gaussian probability distribution, then it should be possible to find the variance of both noise sources using a power spectral density estimate.

$S_{1}(j\omega) = {\sigma_{1}}^{2}$ and $S_{2}(j\omega) = \frac{{\sigma_{2}}^{2}}{\omega^{2}}$, hence you should be able to minimize
$$
J(\sigma_{1},\sigma_{2}) = \| S_{1}(j\omega) + S_{2}(j\omega) - \hat{S}(j\omega) \|
$$
where $\hat{S}(j\omega)$ is the PSD estimate of your output, in order to find $\hat{\sigma}_{1},\hat{\sigma}_{2}$.

The MATLAB script below shows how this might work in practice.

Ts = 0.1; % sampling time
Nts = 1e6; % number of samples
n = randn(2,Nts); % two Gaussian noise sequences
t = (0:length(n)-1)*Ts; % time vector

sd = 1e3.*rand(1,2); % random standard deviations for noise signals

% filter for the noise
G = sd.*[tf(1,[1 0]), tf(1)];
% zero-order hold
Gd = c2d(G,Ts,'zoh');
% time-response for G(s) by the two Gaussian white noise sources
y = lsim(Gd,n,t);
y = y(1e2:end); % remove initial transient (only want stationary response)

% compute a PSD (of stationary response)
[Pxx,f] = pwelch(y,blackmanharris(length(y)/100),[],[],1/Ts,'psd','onesided');

% scaling using logarithm provides better numerical stability when
% minimizing the function
Jfun = @(p)(log(p(1)*1./(2*pi*f(2:end)).^2 + p(2)) - log(Pxx(2:end)));

p0 = 1e3.*rand(1,2); % random initial values for curve fit
p = lsqnonlin(Jfun,p0); % minimize the cost-function Jfun

sdh = sqrt(p/(2*Ts)); % estimate of standard deviations
est_error = norm(sd - sdh);

disp('Ideal parameters')
disp(sd)
disp('Estimated parameters')
disp(sdh)
disp('Estimation error')
disp(est_error)

fa = 10*log10(p(1)*1./(2*pi*f(2:end)).^2 + p(2));
figure(2), semilogx(f,10*log10(Pxx),f(2:end),fa), grid on
xlabel('Frequency'), ylabel('Power'), legend('PSD','Fit')


Shannon interpolation formula for downsampled data with an "almost ideal" low pass filter Let $x[n]$ be a discrete time signal with DFT given by $X(f)=\sum_n x[n]e^{-2\pi inf}$ supported on $[-1/2M,1/2M]$ with $f\in[-1/2,1/2]$.

I can then down-sample to get $y[n]:=x[nM]$. Then, let

$$\widetilde{x}[n]=\begin{cases}My[n/M],& M|n,\\0,&\text{otherwise}.\end{cases}$$

Then its DFT is given by

$$ \begin{aligned}\widetilde{X}(f)&=\sum_{n\in\mathbb{Z}}\widetilde{x}[n]e^{-2\pi inf}\\&=\begin{cases}M\sum_{n\in\mathbb{Z}}y[n/M]e^{-2\pi inf},& M|n,\\0,&\text{otherwise}\end{cases}\\&=\begin{cases}M\sum_{n}x[n]e^{-2\pi inf},&M|n,\\0,&\text{otherwise}.\end{cases}\end{aligned} $$

Now, let $\hat{x}[n]=(\tilde{x}\ast h)[n]$ be the discrete Hilbert transform of $\tilde{x}$, with $h$ an "almost" ideal low-pass filter with cut-off frequency $f_c=1/M$.

My question is, how do I then apply the Shannon interpolation formula to reconstruct $x(t)$?

Intuitively, I would guess that it would be something along the lines of

$$x(t)=\left(\sum_{n\in\mathbb{Z}}x[n]\cdot\delta(t-n\Delta t)\right)\ast H(f),$$

with

$$H(f)=\begin{cases}\frac{DTFT\{\hat{x}[\cdot]\}(f)}{M\cdot X(f)},&\text{if }M|n, \\ 0,&\text{otherwise}.  \end{cases}$$

Am I correct?

I don't get your downsample step when you downsampled by factor $M$.

Let me go from scratch with the spectrum visualization below, with time domain, continuous frequency domain and discrete frequency domain from left to right.

When we reduce the sampling frequency by a factor $k$, the signal spectrum is copied to new replicas at $f_s/k$. The discrete spectrum is a snapshot of continuos spectrum at $[-f_s, f_s]$ and $[-f_{s,down} = f_s/k, f_{s,down} = f_s/k]$, so it is expanded by the factor $k$.

The DFT works on the discrete frequency domain. The downsample $k$ must be chosen so that the expanded spectrum (in discrete frequency, or the new replicas in continuous version) does not overlap with their copies centered at $-2\pi$ and $2\pi$, or $-1/2$ and $1/2$ if you take discrete instantenous frequency by dividing $2\pi$; otherwise aliasing happens and the reconstruction of $x(t)$ is impossible. If there is not alias, a low pass filter to take the spectrum part centered at 0 is suffice. This filter does not need to be "ideal", just to be sure that the filter takes only the center part.

In your calculation, if I understand well, you are talking about the discrete frequency domain and the original spectrum takes $1/M$ the normalized band (you said $f$ is in $[-1/2, 1/2]$ and your signal has support in $[-1/2M, 1/2M]$). In this case your downsample factor $k$ must be less than $M/2$ and the cutoff frequency of your ideal filter must be at $1/4$ if $k=M/2$ (and yes, in this case $k=M/2$, we need "ideal filter" assumption).


Discrimiation between human speech and computer generated one I'm currently working on voice activity detection (VAD) theme and I just started to wonder - if computer generated speech (as when using voice generators) can be discriminated from human speech with classic VAD approaches (ex. 4Hz energy modulation, zero-crossing ratio, modulation entropy, short time energy, MFCC coefficients).

Or - from mathematical point of view - are they just the same signals (or with very similar characteristics)?

Many current computer speech synthesizers seem to produce less jitter or variability in inter glottal closure periodicity timing than most typical human speech.  (e.g. humans often produce subtle vibrato as well as tremolo, but synthesizers could change to do this as well.)


Why are edges in spatial images represented as edges in their Fourier transform image? Here is a well-known image and its Fourier Transform (magnitude).

If I understand correctly the theory behind the FFT, each pixel in the FFT image represents a certain 2D sine wave with frequency depending on the distance from the center of the image, and orientation depending on angle with the horizontal. Intensity of these pixels indicates the coefficients with which each sine wave is added, which, combined with (hidden here) phase information, gives the original image back if we do an inverse Fourier transform.

While I fail to grasp is that when treating real images such as the ones presented, why are edges in the images visible as edges in the frequential domain ?
In this example, there is a diagonal line in the FFT image (let's forget about the vertical and horizontal line which I think are artefacts based on the way FFT is computed, needs a periodical image, etc). This diagonal line is probably caused the girl's hat. But as I understand it, the line in FFT domain means a sum of sines oriented in the same way but with different frequencies. How does that result in an edge when we convert back ? Since edges are high frequency information, wouldn't an edge be represented by one very bright point in the FFT ?
Does it have to do with adding different sines so that they cancel each other out on some portions of the image ? Does the phase image has anything to do with it ?

A more compelling example might be the following image set :

An edge can be described as a sum of sinusoids of different frequencies with phase all equal to $\pi/2$ (or $-\pi/2$) at the location of the edge.

You might find this paper interesting, it is about how different combinations of sinusoids make up different types of features.

Svetha Venkatesh, Robyn Owens On the classification of image features


Time domain and frequency domain analysis equivalence Let us imagine an LTI system with physically realizable input (ruling out fancy mathematical functions and the concomitant complexities and paradoxes) completely known from -$\infty$ to $\infty$.  We want to calculate the output. We can analyse it in time domain using the linear constant coefficient differential equation and in the frequency domain using the Fourier representation of the signal and frequency response of the system.

The outputs calculated from both the methods should be identical. In time domain, we get (transients + steady state)  response of the system,  In frequency domain, the sinusoids are eternal and hence transients are generally neglected.  But since we are talking about the same system with the same input, the outputs calculated from both the methods should be the same.

Will they be identical in spite of neglecting the transients in the frequency domain?

now, even a steady-state sinusoid can be thought of, in the limit, as a sum of weighted pulses, each with a beginning and some kinda end.  how do all of these "transient" signals add up to a steady-state sinusoid?  but they do:

$$\begin{align}
x(t) &= \lim_{T \to 0} \sum_{n=-\infty}^{+\infty} x(nT) \cdot \operatorname{rect}\left( \frac{t-nT}{T} \right) \\
 &= \lim_{T \to 0} \sum_{n=-\infty}^{+\infty} x(nT) \cdot \frac{1}{T}\operatorname{rect}\left( \frac{t-nT}{T} \right) \cdot T \\
 &= \int\limits_{-\infty}^{\infty} x(u) \cdot \delta(t-u) \cdot du
 \end{align}$$

where $ \operatorname{rect}(u) \triangleq \begin{cases} 1 \quad |u|<\tfrac12 \\ 0 \quad |u| > \tfrac12  \end{cases} $

The Fourier integral models the entire time function as an infinitely dense and infinitely large set of sinusoids.  in a similar manner, those steady-state functions can add to something with a transient.

$$ x(t) = \int\limits_{-\infty}^{+\infty} X(f) \, e^{j2\pi ft} \ df $$


Impulse response convolution and normalization2 when I take inverse Laplace transform of a system transfer function \

Lets say LPF whose TF is

$$\frac{Y(s)}{X(s)} \triangleq H(s) = \frac{W}{s+W} $$

the inverse Laplace/impulse response is

$$h(t) = We^{-Wt}u(t) $$

where $u(t)$ is the Heaviside unit step function:

$$ u(t) \triangleq \begin{cases}
1 \qquad & t \ge 0 \\
0 \qquad & t < 0 \\
\end{cases} $$

Now to see the system response to a square wave $x(t)$ with $|x(t)| = 1$, I need to convolve

$$ y(t) = h(t) \star x(t) $$

Now if you look at $h(t)$, the maximum amplitude of $h(t)$ is

$$\max{|h(t)|} = W $$

Then $y(t)$ is amplified by $W$ times.

So what is happening? how do I normalize this?

Should I normalize this by $\max{|h(t)|}$ or
Should I normalize this with

$$W_{z_1} W_{z_2} \cdots W_{z_n}/(W_{p_1} W_{p_2} \cdots W_{p_n})$$

(product of zeroes)/(product of poles) of transfer function?

Why is this even happening?

I think you're likely forgetting how the anti-derivative of $h(t)$ affects the gain of the convolution operation.  Recall that somewhere in your convolution integral, you'll be taking an integral of the form $\int We^{W\tau} d\tau$.  The Chain Rule requires the  $W$ in the exponent must appear as a $1/W$ factor after integrating.  This factor cancels the $W$ multiplier in $h(t)$ giving unity gain at DC.


Audio signal series' beginning points(about 50) distort after a bandpass filter I am using some bandpass filters to detect a certain freq of the sampled audio series in noisy environment. Everything is OK except this:
I found that for many bandpass filters I tried(butterworth, window based FIR etc.) The audio signals' beginning points(about 50~100) always distort after the filter. For example:
[![enter image description here][1]][1]

I am new to DSP so maybe I didn't express in an scientific way... What I mean by "distortion" is that in some cases these beginning points should have been filtered out but they were not, and in other cases vice versa.After these distorted points, everything becomes fine...
Could u plz help me on this issue?
Thank you very much

To make this question more clear, below is a plot consists of 2 waveforms, the blue one is the original audio signal sampled, and the brown one is the sampled signal passed butterworth filter.It doesn't like filter delay stuff.

That's the filter's output transient, whose length is dependent on the effective length of the filter impulse response. If you are performing block based filtering, then only the very first block should exhibit it and later subsequent blocks should not display this transient, due to properly set initial values.


Help me understand the phase of an FIR filter I'm going over an exercise of designing an simple FIR filter, I have the solution in front of me, but I struggle to understand some parts of it.

I'm asked to have an FIR filter with only 3 coefficients.

I have reached this form:
$$ H(\theta ) = (2\alpha_0 \cos(\theta)+\alpha_1)e^{-j\theta} $$
where $\theta$ is the frequency variable (DTFT).

The magnitude is easy:
$$ \left | H(\theta ) \right |=\left | 2\alpha_0 \cos(\theta)+\alpha_1 \right | $$

But I struggle with the phase, it seems easy: $\phi(H(\theta))=-\theta$ , but I guess this is wrong ? why ?

The solution that is given is:

$$\phi(H(\theta)) = -\theta +\beta$$
$$\beta = 0 \text{ for } 2\alpha_0 \cos(\theta) + \alpha_1 > 0$$

$$\beta = \pi \text{ for } 2\alpha_0 \cos(\theta) + \alpha_1 < 0$$

EDIT:

I think I just answers figured it out, if that expression is negative we can write it as positive with a phase of π. What happens if $\alpha_1$ is larger than $4\alpha_0$ thus the expression is always positive ? does $\beta =0$ ?

Your answer does not take into account the possible sign of the magnitude term in the polar description of the frequency response $H(\theta) = A(\theta)e^{j\phi(\theta)}$. I guess you mix this definition with the more usual definiton of  $H(\theta) = |A(\theta)|e^{j\phi(\theta)}$ where the magnitude will always be positive and therefore phase will be what you predicted, however for your example $A(\theta)$ may assume a negative value and you shall reflect this as an addition of $\pi$ into the phase where $e^{j(\theta + \pi)} = -e^{j\theta} $


How to convert between function forms of single frequency signals I have an equation x(t) = (3+4j)*e^(j*8*pi*t)+(3-4j)*e^(-j*8*pi*t) and am trying to convert between this and shifted sines and cosines. I'm attempting mash this through euler's formula and got a 10 for the peak amplitude. Omega is 8*pi, but I'm struggling with the finding the phase shifts for the sine and cosine and the amplitudes on the individual sines/cosines: A*cos(2*pift)+B*sin(2*pift).

Is this making sense? Can anyone point me in the right direction? Please and thank you!!


Implementing a bandpass filter in R - this code's logic confuses me I am reading this "how to" in R: http://www.inside-r.org/howto/time-series-analysis-and-order-prediction-r

I am confused about the tutorial's implemented band pass filter. In particular, here's the code I'm wondering about. The goal is to extrapolate a time series based only on certain ranges of frequencies determined by peakind.

subsignals <- lapply(c(peakind$freqindex, midindex+1), function(x){
    upperind <- x
    fsub <- f
    notnullind <- ((fsub$freqindex >= lowerind
                    & fsub$freqindex < upperind)
                   |
                    (fsub$freqindex >  (lindex - upperind + 2)
                     & fsub$freqindex <= (lindex - lowerind + 2)))
    fsub[!notnullind,"coef"] <- 0
    lowerind <<- upperind
    Re(fft(fsub$coef, inverse=TRUE)/length(fsub$coef))
})

My question is, why this condition:

 (fsub$freqindex >  (lindex - upperind + 2)
                         & fsub$freqindex <= (lindex - lowerind + 2))

This is basically saying to include (for computation of inverse FFT) those frequencies that are within the narrow range of the top frequency minus upper band frequency +2 and top frequency minus lower band frequency + 2.
So basically just some tale at the top of the frequency. This doesn't make any sense to me. Can someone tell me where I should read up on this or what an intuitive hand-wavey sort of justification is?

Without knowing what values are being passed in, I will assume the +2 is an adjustment for R using 1 based indices.
The band in the upper half of the frequency data represent complex frequencies symmetric about the fft length/2. The motivation is using a property of the DFT that will result in a real value output from the ifft: http://en.wikipedia.org/wiki/Discrete_Fourier_transform#The_real-input_DFT


Monotonic, symmetrical soft-clipping polynomial Okay, I intend to answer this question myself, but I would like other folks to pipe in on it.  The issues or requirements are this:

Input to soft-clipper is $x(t)$, output is $y(t)$.
Output $y(t)$ limited in magnitude to $1$: $|y(t)| \le 1$
Memoryless. $y(t)$ only a function of the present $x(t)$, no past values.  So we may as well call the input $x$ and the output $y$. And we will call the memoryless soft-clipping function: $ y \triangleq f(x) $.
$f(-1)=-1$ and $f(1) = 1$.
$f(-x) = -f(x)$, an odd-symmetry (about $x=0$) function (so polarity of the input is irrelevant).
Monotonically increasing function: $f'(x) \ge 0 \quad \forall x$
For $|x| \le 1$, then $f(x)$ is a finite-order polynomial (so the frequency components generated are finite in frequency) and an odd-order polynomial (so all even-order terms are zero because of the odd-symmetry).  Let's call that odd order $2N+1$ where $N\in \mathbb{Z} \ge 0$.  Being a finite-order polynomial means a limit to the width of spurious frequency components generated which tells us what oversampling ratio suffices.
$f(x)$ to have as many derivatives equal to zero as possible where $x=\pm 1$.  So at $\pm 1$, this splices to a constant function ($f(x) = \operatorname{sgn}(x)$ for $|x| \ge 1$) with as many derivatives being continuous as possible.

So what is $f(x)$?

$$ f(x) = \begin{cases}
-1 & x \le -1 \\
\sum\limits_{n=0}^{N} a_n \, x^{2n+1} \quad & -1 \le x \le +1 \\
+1 & +1 \le x \\
\end{cases}$$

$$ $$

$$\begin{align} y(t) & \triangleq f\big(x(t)\big) \\
& = \sum\limits_{n=0}^{N} a_n \, \big(x(t)\big)^{2n+1} \\
& = x(t)\sum\limits_{n=0}^{N} a_n \, \big(x(t)\big)^{2n} \quad |x(t)| \le 1
\end{align}$$

$f(x)$ is continuous everywhere and as many derivatives as possible are continuous everywhere and the only discontinuity in any derivative is at $x = \pm 1$.

What are the odd-order polynomial coefficients $a_n$?

Okay, so to satisfy requirements 5, 6, 7, and 8, the soft-clipping function will have the form:

$$f(x) =
\begin{cases}
-1 & x \le -1 \\
K \int\limits_{0}^{x} \big(1 - u^2 \big)^N \ du \quad & -1\le x\le 1 \\
+1 & 1 \le x \\
\end{cases}$$

where $K$ is a constant judiciously chosen to satisfy requirement 4.

The 1st derivative of $f(x)$ is

$$ f'(x) = K \big(1 - x^2 \big)^N $$

The 2nd derivative of $f(x)$ is

$$ f''(x) = K N \big(1 - x^2 \big)^{N-1} (-2x)$$

The 3rd derivative of $f(x)$ is

$$\begin{align}
 f'''(x) &= K N(N-1) \big(1 - x^2 \big)^{N-2} (-2x)^2 \ + \ K N \big(1 - x^2 \big)^{N-1} (-2) \\
&= K \big(1 - x^2 \big)^{N-2} \bigg( N(N-1)(-2x)^2 - 2N \big(1 - x^2 \big)  \bigg) \\
\end{align}$$

and, for $n \ge 1$, the $n$th derivative is

$$ f^{(n)}(x) = K \big(1 - x^2 \big)^{N-n+1} \, g_n(x) $$

where $g_n(x)$ is some ($n$-1)th order polynomial function of $x$ and is finite in value.  This can be proven inductively by considering the ($n$+1)th derivative:

$$\begin{align}
\\
f^{(n+1)}(x) &= \tfrac{d}{dx} \Big( f^{(n)}(x) \Big) \\
&= \tfrac{d}{dx} \Big(K \big(1 - x^2 \big)^{N-n+1} \, g_n(x)\Big) \\
&= K(N-n+1)\big(1 - x^2 \big)^{N-n} (-2x) g_n(x) \, + \, K \big(1-x^2 \big)^{N-n+1}g'_n(x) \\
&= K\big(1 - x^2 \big)^{N-n} \Big( (N-n+1)(-2x) g_n(x) \, + \, (1-x^2) g'_n(x) \Big) \\
&= K\big(1 - x^2 \big)^{N-n} \, g_{n+1}(x) \\
\\
\end{align}$$

where $ g_{n+1}(x) = (N-n+1)(-2x) g_n(x) + (1-x^2) g'_n(x)  $.

Because of differentiation, the order of polynomial $g'_n(x)$ is one less than the order of $g_n(x)$, but polynomial $(1-x^2)g'_n(x)$ is one order greater than $g_n(x)$ and so also is $(-2x) g_n(x)$.

When $x = \pm 1$, then the first $N$ derivatives are zero,
$$ f^{(n)}(x) \Bigg|_{x=\pm 1} = 0 \qquad \text{for } 1 \le n \le N $$
making this polynomial maximally flat at $x = \pm 1$.

The integrand is a binomial and can be expressed as a power series using binomial expansion:

$$\begin{align}
\big(1 - u^2 \big)^N & = \sum\limits_{n=0}^{N} \binom{N}{n} \big(-u^2\big)^n (1)^{N-n} \\
  & = \sum\limits_{n=0}^{N} \frac{N!}{n!(N-n)!} \big(-u^2\big)^n (1)^{N-n} \\
  & = \sum\limits_{n=0}^{N} \frac{N!}{n!(N-n)!} (-1)^n u^{2n} \\
\end{align}$$

So the integral can be expressed as an integral of a power series:

$$\begin{align}
\int\limits_{0}^{x} \big(1 - u^2 \big)^N \ du & = \int\limits_{0}^{x} \sum\limits_{n=0}^{N} \frac{N!}{n!(N-n)!} (-1)^n u^{2n} \ du \\
& = \sum\limits_{n=0}^{N} \frac{N!}{n!(N-n)!} (-1)^n \int\limits_{0}^{x} u^{2n} \ du \\
& = \sum\limits_{n=0}^{N} \frac{N!}{n!(N-n)!} \frac{(-1)^n}{2n+1} u^{2n+1} \Bigg|_0^x \\
& = \sum\limits_{n=0}^{N} \frac{N!}{n!(N-n)!} \frac{(-1)^n}{2n+1} x^{2n+1} \\
& = x \sum\limits_{n=0}^{N} \frac{N!}{n!(N-n)!} \frac{(-1)^n}{2n+1} \big(x^2 \big)^n \\
\end{align}$$

When $x = \pm 1$, we get

$$\begin{align}
 \int\limits_{0}^{\pm 1} \big(1 - u^2 \big)^N \ du &= \pm \sum\limits_{n=0}^{N} \frac{N!}{n!(N-n)!} \frac{(-1)^n}{2n+1} \\
  &= \pm \frac{4^N}{2N+1}\binom{2N}{N}^{-1}\\
\end{align}$$

(there is a little improperness in the integral limit.)  The bottom equality is because of this Appendix proof.  So the scaler $K$ must be

$$\begin{align}
 K &= \left( \frac{4^N}{2N+1}\binom{2N}{N}^{-1} \right)^{-1} \\
\\
 &= \frac{2N+1}{4^N}\binom{2N}{N} \\
\\
 &= \frac{2N+1}{4^N} \frac{(2N)!}{(N!)^2} \\
\\
 &= \frac{(2N+1)!}{4^N (N!)^2} \\
\end{align}$$

This makes the entire soft-clipping function to be:

$$ f(x) = \begin{cases}
-1 & x \le -1 \\
\\
x \cdot \frac{(2N+1)!}{4^N \ N!} \sum\limits_{n=0}^{N} \frac{(-1)^n}{(2n+1)\,n!\,(N-n)!} \big(x^2 \big)^n \quad & -1 \le x \le +1 \\
\\
+1 & +1 \le x \\
\end{cases} $$

and it appears that the odd-power (that is the power = $2n+1$) polynomial coefficients are

$$ a_n = (-1)^n \frac{(2N+1)!}{4^N \, N! \,(2n+1) \, n! \, (N-n)!} $$

The polynomials (without splicing to the $\operatorname{sgn}(x) = \pm 1$ saturated components) look like

I think the order $2N+1$ starts at 1 and goes to 9 (or $0 \le N \le 4$)

With the saturation attached, the curves look like

The soft-clipping function is continuous everywhere and all derivatives, up to the $N$th derivative are continuous everywhere and the ($N$+1)th derivative and higher are continuous everywhere except at the splices where $x = \pm 1$.

Here are the same set of curves but with the scaling adjusted so that the slope around $x=0$ (or "gain") remains at 1 (or "0 dB gain").

The point of discontinuity (where the polynomial is spliced to a constant $\pm 1$) is at $x= \pm K$.  I really don't think one needs to get over a $7$th-order ($N=3$) softclipper.  Oversampling by a factor $4$x suffices to deal with this $7$th-order polynomial and prevent any aliasing to the original baseband before upsampling.


Finding preferred pair of $m$-sequence for Gold code I would like to find a preferred pair of $m$-sequence of degree 34 (not a multiple of 4) for Gold code. However, I am lost when I am checking William Wesley Peterson's table in his Error-correcting Codes, Appendix C

I tried first fix the first entry for my first $m$-sequence. To find the other sequence, since $n = 34$, it is even, then
$$2^{\left(\frac{34+2}{2}\right)}+1 = 2^{18} + 1$$
or
$$2^{\left(\frac{34-2}{2}\right)}+1 = 2^{16} + 1$$
But either number is not mapped to any entry in the table that I cannot proceed to find the other $m$-sequence. I assume I should find the other $m$-sequence which has the equivalent root of
$$r^{\left(2^{\left(\frac{34+2}{2}\right)}+1\right)}$$
or
$$r^{\left(2^{\left(\frac{34-2}{2}\right)}+1\right)}$$
To elaborate, the following is what I consider as an easy case of prefered pair of Gold Sequence using Peterson's table and Gold's theorem, presented in this paper. I have no access to Dixon's book therefore I don't know what the correct procedure is.

The author applying the equation in the conjugate form of Gold's theorem in Gold's paper

Printed tables of irreducible polynomials of large degrees are usually incomplete because there are so many of them!  What one can do is use various tricks of the trade to figure out a set of $68 = 2\times 34$ consecutive bits of one m-sequence from the other m-sequence (with known generator polynomial as taken from the Peterson and Weldon table) and then run the Berlekamp-Massey algorithm on these consecutive bits to get the shortest shift register that generates the desired m-sequence.
If ${\bf s} = s_0, s_1, s_2, \ldots $ is one of the m-sequences, then the other m-sequence ${\bf t} = t_0, t_1, t_2, \ldots $ is related to ${\bf s}$ via
$t_i = s_{qi}$ where $q = 2^{16}+1 = 65,537$. See the paper "Cross-correlation properties of pseudorandom and related sequences," Proc. IEEE, vol.68, pp.593-619, May 1980 (unfortunately behind IEEE's paywall, but copies can be found in the Internet) if you are unfamiliar with this idea. So, we crank up the LFSR that generates ${\bf s}$ and run it to produce the sequence ${\bf s}$, recording $$t_0 = s_0, ~~t_1 = s_{65537},~~ t_2 = s_{131074},~~ \cdots, t_{67} = s_{4390979}.$$  Once upon a time, cranking out 4 million+ bits would have taken some effort, but these days, it's no big deal.
What next? Well, the Berlekamp-Massey algorithm can be viewed as a shift-register synthesis algorithm: it finds the shortest LFSR that can generate any given sequence. Thus, applying the Berlekamp-Massey algorithm to the first $68$ bits of $\bf t$ which we have just found gives us the LFSR that generates the m-sequence $\bf t$, and away we go!!


Can I test any frequency with DFT? Ladies, Gentlemen,

Please let me the subject question, for I see in DFT equations testing frequencies sub-multiples of sampling frequency. Regards.

Yes, what you've observed is correct:

The $N$-point DFT only "compares" a signal to the set of complex oscillations $e^{j2\pi f n}$ for frequencies $f$ that fit an integer number of times into the DFT length $N$.

Any other frequency will not be mapped "sharply" to a single DFT bin, but will have some energy in the surrounding bins and spectral repetitions.

We use windows on the input signal of a DFT to shape that "bleeding" a bit.

It is very common to increase the "resolution" by having a greater $N$. To do that, you simply take $K<N$ samples of your input, and attach $N-K$ zeros, so that you end up with a vector of $N$ samples. That way, you can have a finer "grid" of frequencies.


Why is no phase information lost when a received signal goes through a mixer? I'm sorry if the question is too philosophical or makes no sense, but I need to be able to explain why it doesn't make sense.
When we take the in-phase and quadrature-phase base-band components of an intermediate frequency signal, it's very easy to show that any one of the two orthogonal components alone would be insufficient to distinguish the variations in a single-carrier phase from variations in its amplitude. This tutorial makes it quite clear. The Q component, although real, carries the imaginary component of the phasor $A(t)\exp(2\pi f_I\cdot t + \phi(t) )$.

However, when we pass the RF signal through a mixer, there is no Q component lost. No phase modulation will be lost in the base-band signal obtained later. Why not?
A figure to help create context.

The receiver's mixer introduces a phase-shift to the in-phase and quadrature components. In very general terms, assume we have a signal $$a(t)e^{j(2\pi f_0 t + \phi(t))},$$ which is mixed with a local oscillator of frequency $f_{lo}$ and phase $\delta$: $$e^{j(2\pi f_{lo}t + \delta)}.$$ The output signal is

\begin{align}
a(t)e^{j(2\pi f_0 t + \phi(t))}e^{j(2\pi f_{lo}t + \delta)}
&= a(t)e^{j(2\pi (f_0+f_{lo})t + \phi(t) + \delta)} \\
&= a(t)e^{j(2\pi (f_0+f_{lo})t + \phi(t))} e^{j\delta}.
\end{align}

In other words, even though the frequency has been shifted and a phase-shift equal to $\delta$ has been introduced, $a(t)$ and $\phi(t)$ are still present in the signal and that means that all the information carried by the original signal is still there.

Usually the phase-shift is corrected in the digital back-end in the receiver, or you can use a differential encoding (such as DBSK), in which the phase-shift becomes irrelevant.


Why are FIR filters still stable even though they contain poles? How come FIR filters are always stable?
Since they contain poles, shouldn't they be more affected by stability issues than others?

FIR filters contain as many poles as they have zeros.  but all of the poles are located at the origin, $z=0$.

because all of the poles are located inside the unit circle, the FIR filter is ostensibly stable.

this is probably not the FIR filter the OP is thinking about, but there is a class of FIR filters called Truncated IIR filters (TIIR) which may have a pole on or outside the unit circle that is canceled by a zero at the same location.  the simplest example of this is the moving sum or moving average filter.  but, from an I/O perspective, these TIIR filters are FIR.

but i wouldn't naively guarantee "stability".  using control-system language, the TIIR filter are not "completely observable" and may appear stable because its impulse response appears finite in length, but inside the filter states could be going to hell, and with finite numerical precision, that internal instability will eventually show up at the output.

we have to disabuse ourselves from the notion that "FIR filters have no poles".  ain't true.


Get the frequency response curve from FIR filter coefficients + sampling rate What's the fastest way (if possible in the browser thanks to an online tool, or if not possible easily, with Python), to get the frequency response curve (x-axis: Hz, y-axis: dB), when giving just:

the sampling rate, example: 44100
the coefficients, example: 1 and -1, i.e. y[n] = 1 * x[n] + (-1) x[n-1] = x[n] - x[n-1]

?

This would be often useful to quickly have an eye on the frequency response curve of a filter.

Example:

As suggested by a comment, scipy.signal.freqz gives the solution. Don't forget to use this to do the rad/sample to Hz conversion.

Here is a ready-to-use code, posted here for future reference:

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
sr = 44100
w, h = signal.freqz(b=[1, -1], a=1)
x = w * sr * 1.0 / (2 * np.pi)
y = 20 * np.log10(abs(h))
plt.figure(figsize=(10,5))
plt.semilogx(x, y)
plt.ylabel('Amplitude [dB]')
plt.xlabel('Frequency [Hz]')
plt.title('Frequency response')
plt.grid(which='both', linestyle='-', color='grey')
plt.xticks([20, 50, 100, 200, 500, 1000, 2000, 5000, 10000, 20000], ["20", "50", "100", "200", "500", "1K", "2K", "5K", "10K", "20K"])
plt.show()

With logarithmic x-axis:

Another example with a comb filter (COEF = [1] + [0] * 200 + [-1], i.e. y[n] = x[n] - x[n-201]):


Understanding rho and theta in hough I understood the concept behind Hough transform. One thing i am still unsure is the rho and theta range of hough transform. The range is given below.

If theta is in the range [0,pi), then rho values can have negative
values. If theta is in the range[0,2pi), then rho can have only
positive values.

How can rho be negative ?. Is the range standard for all the images ?. I have referred questions on Hough Transform under stack overflow as well as dsp-stackoverflow, but it wasn't helpful. Can someone explain me bit more clearly please ?.

The Hough transform seeks to replace finding colinear points by  finding concurrent lines.  The two problems are  mathematically equivalent. But dealing with all discrete pixels and a large number of digital lines is complicated and computationally expensive.
Most implementations of the Hough transform discretize the set of distances and angles. With respect to some origin, one can describe a line in polar coordinates by:

a positive distance and an angle in $[0,\,2\pi[$
a signed distance and an angle in $[0,\,\pi[$

This is explained in Polar coordinate system/Uniqueness of polar coordinates:

any polar coordinate is identical to the coordinate with the negative
radial component and the opposite direction

History: the rho-theta parametrization universally used today was first described in Duda, R.O.; Hart, P. E. (1972). "Use of the Hough Transformation to Detect Lines and Curves in Pictures".


Speech recognition algorithm in GNU Octave to convert speech to text (and vice versa) Recently, I am building an application in Octave and I saw some algorithms about speech recognition spread around the net, but most of them were implemented in Python, MATLAB, C++ and no one in Octave. I already came to test the MATLAB one in the Octave IDE but there are lots of packages that exists in MATLAB and there aren't in Octave, I have tried already to look for at the net in order to find anyway to run extensions to run Python functions in Octave for example but I've only found simple ones to return list of numbers, these things...
Is there any possibility to implement a speech recognition integrated with google voice that converts speechs in text in GNU Octave?

Is there any possibility to implement a speech recognition integrated with google voice that converts speechs in text in GNU Octave?

Yes.
The Google Cloud command line utility returns results in JSON. Octave is working on its JSON support, but until then you can use something like this.
What we are missing now is triggering the whole process from within Octave.
This can be achieved with the Subprocess control functions.
Usually, you can redirect the stdin, stdout streams towards buffers in your application, which would get you the JSON in a local variable within Octave from which then you can call its parser to make it computable. If you cannot, you would have to create an intermediate file and read any data through that. Which might require some timing control too.
This is the general way of going about it even across languages (e.g. Python) and exchanging JSON is probably the better option than just returning some data over a stream and relying on putting it back together on the other side.
But this is also slow which your users (subjects) might find annoying. There is the overhead of starting up the process and in this case, there will be data exchanged over the network too. Besides, the free version of this google service allows a limited number of accesses per month. So, overall, it might be better to check if you can handle that recognition bit locally. You would still probably have to call an external to Octave package (...indeed, most probably Python with the same way outlined here), but the round-trip will be quicker (and cheaper).
Hope this helps.


Deconvolution in Python I'm trying to use and understand SciPy's deconvolve for a project I'm working on. I'm having some trouble understanding how to use it.

What I would like to do is to take two PMFs from discrete gaussian distributions and recover an unknown distribution using deconvolution. My understanding is that if I deconvolve the PMF from ~N(10, 1) and the PMF from ~N(30, 2), I should recover the PMF for the distribution ~N(20, 1). How can I do this in SciPy. It would be great, if I could see an example. Thanks in advance for the help.

The function is based on Matlab's deconv, so reading that page should help understand it.

Here's a docstring I wrote for SciPy's deconvolve, but haven't submitted yet because I'm not sure it's 100% correct: https://github.com/scipy/scipy/pull/430#issuecomment-13675004

The input to deconvolve is signal and divisor, and your output is quotient and remainder, where signal was originally produced by signal = convolve(divisor, quotient) + remainder.

original = [0, 1, 0, 0, 1, 1, 0, 0]
impulse_response = [2, 1]
recorded = scipy.signal.convolve(impulse_response, original)
print recorded
# [0 2 1 0 2 3 1 0 0]
recovered, remainder = scipy.signal.deconvolve(recorded, impulse_response)
print recovered, remainder
# [ 0.  1.  0.  0.  1.  1.  0.  0.] [ 0.  0.  0.  0.  0.  0.  0.  0.  0.]

I don't understand what you're trying to do with PMFs, though, so I can't give a more specific example.


Is Kalman Filter for LTV systems applicable on Hybrid systems? As all of us know Kalman Filter is designed for LTV (Linear Time Variant) systems. But nobody in the literature applies that on Hybrid systems. In my opinion Linear Hybrid systems are a subset of Linear Time Variant (LTV) systems, So my question is why don't people use Kalman Filter for Hybrid systems with linear modes and resort to robust estimation approaches?

You really have two questions -- one, is the Kalman filter applicable to linear hybrid systems, and the other is why isn't it used.
I can't answer your second question.
As for your first question (is it applicable) -- I think it could be made to be applicable.  In fact, I strongly suspect that there's some fairly obvious combination of a Kalman-Bucy filter for the continuous-time part of the system and a straight Kalman filter for the discrete-time part.
However
I've found that for a lot of real-world problems involving a discrete-time controller and a continuous-time plant, it's easier to just model the continuous-time plant as a discrete-time system as seen by the controller, optimize the control in discrete time, and then just sweep the fact that the plant is actually continuous time under the rug.
This is because if you're going to make a discrete time controller that actually works in practice, the most reliable way of achieving this end is to make sure that the frequencies of any significant plant dynamics are exceeded by the controller sampling rate by a healthy margin.  Do this, and any inaccuracies that you introduce by assuming that the plant responds step-wise to controller commands is less than any inaccuracies from the difference between your best-ever attempt at a plant model and what a real plant is going to do under real circumstances.
However (again)
If you're casting about for a PhD topic, there may be something in here, unless the approach was considered in 1952 and discarded for some good reason.


Stationary vs non-stationary signals? There are nice technical definitions in textbooks and wikipedia, but I'm having a hard time understanding what differentiates stationary and non-stationary signals in practice?

Which of the following discrete signals are stationary? why?:

white noise - YES (according to every possible information found)
colored noise - YES (according to
Colored noises: Stationary or non-stationary? )
chirp (sinus with
changing frequency) - ?
sinus - ?
sum of multiple sinuses with different periods and amplitudes - ?
ECG, EEG, PPT and similar - ?
Chaotic system output (mackey-glass, logistic map) - ?
Record of outdoors temperature - ?
Record of forex market currency pair development - ?

Thank you.

There is no stationary signal. Stationary and non-stationary are characterisations of the process that generated the signal.

A signal is an observation. A recording of something that has happened. A recording of a series of events as a result of some process. If the properties of the process that generates the events DOES NOT change in time, then the process is stationary.

We know what a signal $x(n)$ is, it is a collection of events (measurements) at different time instances ($n$). But how can we describe the process that generated it?

One way of capturing the properties of a process is to obtain the probability distribution of the events it describes. Practically, this could look like a histogram but that's not entirely useful here because it only provides information on each event as if it was unrelated to its neighbour events. Another type of "histogram" is one where we could fix an event and ask what is the probability that the other events happen GIVEN another event has already taken place. So, if we were to capture this "monster histogram" that describes the probability of transition from any possible event to any other possible event, we would be able to describe any process.

Furthermore, if we were to obtain this at two different time instances and the event-to-event probabilities did not seem to change then that process would be called a stationary process. (Absolute knowledge of the characteristics of a process in nature is rarely assumed of course).

Having said this, let's look at the examples:

White Noise:

White noise is stationary because any signal value (event) is equally
probable to happen given any other signal value (another event) at
any two time instances no matter how far apart they are.

Coloured Noise:

What is coloured noise? It is essentially white-noise with some additional constraints. The constraints mean that the event-to-event probabilities are now not equal BUT this doesn't mean that they are allowed to change with time. So, Pink noise is filtered white noise whose frequency spectrum decreases following a specific relationship. This means that pink noise has more low frequencies which in turn means that any two neighbouring events would have higher probabilities of occurring but that would not hold for any two events (as it was in the case of white noise). Fine, but if we were to obtain these event-to-event probabilities at two different time instances and they did not seem to change, then the process that generated the signals would be stationary.

Chirp:

Non stationary, because the event-to-event probabilities change with time. Here is a relatively easy way to visualise this: Consider a sampled version of the lowest frequency sinusoid at some sampling frequency. This has some event-to-event probabilities. For example, you can't really go from -1 to 1, if you are at -1 then the next probable value is much more likely to be closer to -0.9 depending of course on the sampling frequency. But, actually, to generate the higher frequencies you can resample this low frequency sinusoid. All you have to do for the low frequency to change pitch is to "play it faster". AHA! THEREFORE, YES! You can actually move from -1 to 1 in one sample, provided that the sinusoid is resampled really really fast. THEREFORE!!! The event-to-event probabilities CHANGE WITH TIME!, we have by passed so many different values and went from -1 to 1 in this extreme case....So, this is a non-stationary process.

Sinus(oid)

Stationary...Self-explanatory, given #3

Sum of multiple sinuses with different periods and amplitudes

Self explanatory given #1, #2,#3 and #4. If the periods and amplitudes of the components do not change in time, then the constraints between the samples do not change in time, therefore the process will end up stationary.

ECG, EEG, PPT and similar

I am not really sure what PPT is but ECG and EEG are prime examples of non-stationary signals. Why? The ECG represents the electrical activity of the heart. The heart has its own oscillator which is modulated by signals from the brain AT EVERY HEARTBEAT! Therefore, since the process changes with time (i.e. the way that the heart beats changes at each heart beat) then it is considered non-stationary. The same applies for the EEG. The EEG represents a sum of localised electrical activity of neurons in the brain. The brain cannot be considered stationary in time since a human being performs different activities. Conversely, if we were to fix the observation window we could claim some form of stationarity. For example, in neuroscience, you can say that 30 subjects were instructed to stay at rest with their eyes closed while EEG recordings were obtained for 30 seconds and then say that FOR THOSE SPECIFIC 30 SEC AND CONDITION (rest, eyes closed) THE BRAIN (as a process) IS ASSUMED TO BE STATIONARY.

Chaotic system output.

Similar to #6, chaotic systems could be considered stationary over brief periods of time but that's not general.

Temperature recordings:

Similar to #6 and #7. Weather is a prime example of a chaotic process, it cannot be considered stationary for too long.

Financial indicators:

Similar to #6,#7,#8,#9. In general cannot be considered stationary.

A useful concept to keep in mind when talking about practical situations is ergodicity. Also, there is something that eventually creeps up here and that is the scale of observation. Look too close and it's not stationary, look from very far away and everything is stationary. The scale of observation is context dependent. For more information and a large number of illustrating examples as far as the chaotic systems are concenred, I would recommend this book and specifically chapters 1,6,7,10,12 and 13 which are really central on stationarity and periodicity.

Hope this helps.


Cost function for LTI system identification I am currently reading and trying to understand a paper (Kulkarni and Colburn, 2004) that utilizes system identification methods to approximate head-related transfer functions.

The general approach is to

Compute an autoregressive (all-pole-) estimate of the transfer function using the autocorrelation method for linear prediction.
Use the AR estimate as a starting point to compute a pole-zero-representation of the system transfer function iteratively.
Evaluate the result of the estimation process on a logarithmic scale (Error in dB).

For the iterative procedure, the authors are proposing a cost function

$\hat{C} = \frac{1}{2\pi}\int_{-\pi}^{\pi}|H(e^{j\omega})A(e^{j\omega}) - B(e^{j\omega})|^2 d\omega $,

where $\hat{H}$ is the system transfer function, $A$ is the DTFT of the recursive coefficients and $B$ the DTFT of the transversal coefficients.
I understand this approach originates from this paper (Kalman, 1958).

This cost function is then extended for the iterative process as

$\hat{C}_i = \frac{1}{2\pi}\int_{-\pi}^{\pi}|\frac{H(e^{j\omega})A_i(e^{j\omega})}{A_{i-1}(e^{j\omega})} - \frac{B_i(e^{j\omega})}{A_{i-1}(e^{j\omega})}|^2 d\omega $,

where the index $i$ denotes variables corresponding to the $i$th iteration.
The iterative modification originates from this paper (Steiglitz and McBride, 1965).

In order to find a solution on a decibel scale, a weighting function $W$ is introduced:

$\hat{C}_i = \frac{1}{2\pi}\int_{-\pi}^{\pi} |W_i(e^{j\omega})|^2 |\frac{H(e^{j\omega})A_i(e^{j\omega})}{A_{i-1}(e^{j\omega})} - \frac{B_i(e^{j\omega})}{A_{i-1}(e^{j\omega})}|^2 d\omega $,

which is introduced in the paper as

$W_i(e^{j\omega}) = \frac{\log\left( H(e^{j\omega})\right) - \log\left( \frac{B_{i-1}(e^{j\omega})}{A_{i-1}(e^{j\omega})}\right)}{|H(e^{j\omega})A_{i-1}(e^{j\omega}) - B_{i-1}(e^{j\omega})|^2}$.

I understand this weighting function is the squared error in logarithmic scale between true and approximated transfer function, divided by the first cost function.

However, i have trouble understanding the process of arriving at the iterative cost function for several reasons. I would like to ask the following questions:

Why is the first cost function $\hat{C}$ preferred to, say, $\frac{1}{2\pi}\int_{-\pi}^{\pi} |H(e^{j\omega}) - B(e^{j\omega})/A(e^{j\omega})|^2$ in the first place? what does it do?
In the iterative cost function $\hat{C}_i$, what is the purpose of dividing by the recursive proportion of the last transfer function estimate?
For what reason is a weighting function introduced for logarithmic error minimization, rather than just using its numerator as a cost function directly?

I would really appreciate any help or pointers into the right direction.

The chosen cost function is the mean squared error, i.e., the integral over a squared magnitude of the difference between frequency responses. The function

$$E(e^{j\omega})=H(e^{j\omega})-\frac{B(e^{j\omega})}{A(e^{j \omega})}\tag{1}$$

depends on frequency, so you can't minimize it directly, unless you want to minimize it for exactly one frequency $\omega$, which is of course pointless. You can choose several error measures depending on $E(e^{j\omega})$ given in $(1)$. Two common choices are

$$\varepsilon_1=\max_{\omega}W(\omega)|E(e^{j\omega})|\tag{2}$$

and

$$\varepsilon_2=\int_{0}^{\pi}W(\omega)|E(e^{j\omega})|^2d\omega\tag{3}$$

with some positive weighting function $W(\omega)$. Note that $\varepsilon_1$ and $\varepsilon_2$ given by $(2)$ and $(3)$ do not depend on frequency.

The authors of the paper you refer to chose to minimize the weighted mean square error given by $(3)$. However, instead of using the linear difference $(1)$, they chose to minimize the average logarithmic difference, i.e., the average error on a dB scale.

The problem with the minimization of $(3)$ is that for IIR filters, it results in a non-linear optimization problem, which is much harder to solve directly, and which might also have locally optimal solutions that are far from the global optimum. The cost function $\hat{C}$ in your question is linear in the filter coefficients. The point of the iteration is now to solve a sequence of linear minimization problems (which is simple, just solve a system of linear equations) in order to compute the solution of the originally non-linear optimization problem (the minimization of $(3)$). Note that if convergence is achieved, then $A_{i-1}(e^{j\omega})=A_{i}(e^{j\omega})$, so the cost function $\hat{C}_i$ is identical to the cost function of the original non-linear problem. Yet, only linear minimization problems are solved in each iteration.

I'm not sure I completely understand your last question, but the weighting function is there to change to problem from minimizing the average squared difference between frequency responses to the average squared differences between the logarithms of the magnitude responses. The necessary weighting function is unknown, but if the procedure converges - note that there is no guarantee that it does in all cases - the final weight function is such that the mean squared logarithmic error is minimized. Note that this latter problem is highly non-linear, whereas the proposed procedure only solves linear subproblems.


How to determine the frequency of each row in a MATLAB spectrogram I recently learned how to manually compute a spectrogram of an audio signal, which is essentially a matrix. Each column corresponds to a time frame (e.g. 10 ms or something) but I'm asking which row corresponds to which frequency?

See my implementation below:

clc; clear all; close all;
Fs = 44100;
t_max = 3;
T = 1/Fs;

time = 0:T:(t_max-T);
input = chirp(time,1500,1,8000);

window_length_t = 0.01; %10ms window length
window_length_s = round(0.01 * Fs); %window length in samples
if mod(window_length_s,2) == 0
    window_length_s = window_length_s + 1; %make sure we have odd window size
end

signal_framewise = buffer(input , window_length_s , floor(window_length_s/2));
nfft =((window_length_s-1)/2)+1;

out_buffer = zeros(nfft,size(signal_framewise,2));
for jj = 1:size(signal_framewise,2)
    current_frame = signal_framewise(:,jj).*gausswin(window_length_s);
    dtf = fft(current_frame);
    out_buffer(:,jj) = dtf(1:nfft);
end

My theory: The topmost row must be 0.5*Fs (half of the sampling frequency) because of the sampling theorem. All other rows are linear increments to then match to 0.5*Fs. Example: 5 Rows, topmost is 5 Hz, then the first row is 1 Hz, the second is 2 Hz and so forth.

The physical frequency $f$ that a given Discrete Fourier Transform (DFT) frequency bin $k$ corresponds to is $f = \frac{k}{N_{DFT}} \cdot Fs$.

For an $Fs=44100$ Hz, the middle point of a given $N_{DFT}$ point DFT would correspond to $\frac{1}{2} \cdot Fs$, also known as the Nyquist Frequency.

EXCEPT MATLAB, which maps $1$ to $\frac{Fs}{2}$.

For more information, please see here.

Hope this helps.


Purpose of eigenspace of covariance matrix of a blob? Given a blob of an image (representing an object), according to Wikipedia, we can compute the co-variance matrix using the image moments.

I understand that the eigenvectors of that matrix can be used to represent the principal axes of an oriented bounding box (OBB), but what are the guarantees for that particular orientation? Is the OBB the one with minimal area? And, are the eigenvectors always orthogonal?

The eigenvectors will always be orthogonal. There is no guarantee that the eigenvectors will give you the bounding box with minimal area.


Phase response in a filter for an Audio Application Should I be concerned with the group phase delay which a filter produces? Especially in an audio application? Should I aim at using a linear phase delay filter for example?

i will agree with Hilmar that it can depend on the specific application.

if the application is to essentially losslessly store or transmit audio to later retrieve or receive that audio, including conversions of format (and this includes the A/D and D/A and SRC) then i would say that there is no good reason for a process to not be linear phase (which is constant group delay, constant phase delay, and the two constants are the same number).  there are audio format conversion processes (like a lossy codec) that are meant for storage or streaming that, for the sake of data reduction, have non-constant phase delay or non-constant group delay.

so if the application is to change the audio in some manner, it might not be possible to have the entire transfer function to be linear phase because of feedback.

if there is never any feedback and the audio is passed from one process to another, in cascade similar to the figure above, and all processes are meant to not change the perceived sound, then each process can have constant delay and there will be no waveform distortion due to non-constant phase delay or non-constant group delay (non-constant delay would not be a linear phase filter).  so in this non-feedback case, there is no necessity to ever make different frequency components have different delays, and i would suggest to normally make that application linear phase.

sometimes, there are audio application in which a the process is one of several in a feedback loop, such as a feedback exterminator or maybe some reverberation control, it might be in this application that you might be tuning for a specific non-constant phase or group delay.

there are also applications like musical effects where a non-linear phase shift might not be perceptually relevant, but we know that linear phase shift is perceptually irrelevant when the processed sound is not compared directly with the input sound, so that the delay is detected.


Sampling Rate and SNR Assume I have a signal that oscillates between 10KHz and 9Khz every 500msec (Change happens linearly, so at t=0, it is 10KHz and it goes to 9.5Kz at t=250msec and 9KHz at t=500msec and frequency start to increase and goes back to 10KHz at 1 sec and this goes on)

Assume I sample this signal at 50Khz and I do fft to find out what's the frequency. But due to sampling rate, I have to accumulate some samples where the interest frequency is slowly shifting.

If I were to sample at 500KHz or even 5MHz, can I get a better SNR? (The assumption here is that I can sample for a shorter duration -still several cycles of the target frequency in time domain- and this will provide a less spread signal in a given time and better fft snr. )

My gut tells me, this won't help but I like to ask the crowd.

Your instincts are right. You won't increase SNR by increasing sampling rate because you also have more noise samples. The overall signal energy and noise energy does not depend on sampling rate.

But, I believe you are mixing up two problems. In your case the frequency is not constant, therefore Fourier based analysis yield less than optimal results. Your problem is in the signal itself, not the signal to noise ratio. You can have a look at the high-order ambiguity function (try this or that) or the Wigner-Ville transform. Those are more suitable for linearly changing frequency.


Why wavelets based transmitter/receiver digital signal processing aren't common? I have seen this thread: Difference between Fourier transform and Wavelets

AFAIK there is no common usage of wavelets in the real-time DSP world (excluding image and video processing).

I am curious why is it so?

Where on the other hand Fourier analysis is widely common and all the transmitter/receiver block diagrams that I have seen contains an FFT module in it:

Is it related to runtime complexity issues? I guess not, since Wavelets transform can be formalized as convolution, right?

e.g aren't wavelets a good candidate to replace the Short-time Fourier transform method ?

One key difference between the DFT and discrete wavelet transform (DWT) is that the DFT deals with frequency whereas the wavelet transform deals with "scales". (e.g. you estimate the width of a pulse or a series of pulses in a time series).

A radio signal is more or less a continuous signal when observed locally. A constant waveform signal is very sparse when decoded by FFT. So it is very easy and natural to represent an RF signal using the DFT/frequency domain.

One of the common issues one have to deal with in regards to the DFT is spectral leakage. Unfortunately, the wavelet kernels generally have a good amount of spectral leakage also so this issue cannot be improved by using the DWT.

The WT does have some useful properties. Aside from the examples you listed, there are some applications in time series analysis where the WT trumps the DFT. For example, detecting drift of atomic clocks or measuring changes in tides (see Wavelet Methods, by Percival and Walden). Also, a more limited application: the WT can be used for denoising a signal when you want to preserve the spikes as well as the baseline drift (e.g. to visually make a spectrum look cleaner).

In terms of time-frequency analysis, the Gabor wavelet is probably the closest you get to STFT.

Like the DFT/FFT the DWT is an orthogonal transform. However, the continuous wavelet transform (CWT), which is computed using convolution, is not. It's useful as an analysis tool only.


Why is the Hamming and Hanning window the most widely used for speech spectograms? When taking the Short-Time Fourier Transform of the speech signal, windowing is applied to taper away the signal frame at both ends and mitigate introduced artifact at expense of the maximum obtainable resolution. Window design is then motivated controlling the trade-offs or spectral-time resolution, and various figures of merits, like the mainlobe width, sidelode width, sidelobe roll-off per octave is introduced in the literature.

In terms of a speech signal, usually a Hamming or Hanning window is employed, and the only justification I have found for it is that speech signals are narrow band signals. How could I demonstrate that a Hamming window is indeed better than, say, a rectangular window for a speech signal?

In a typical natural speech application, a speaker will exhibits frequencies between 0 and 200 Hz, a little higher perhaps for females. Utterances are usually a lot longer than that to produce required frequency, so we can afford to take more frames. I don't know a required resolution for speech. Thus, is it always better to use a Hamming/Hanning window on speech? Is there a clear reason I should tell someone using a rectangular window on producing spectograms, why he/she should not do it (apart from reproducibility of certain results in academic literature)?

any window (almost any) is better than a boxcar window, unless your spectrum is flat.  speech signals are not flat.

I suspect that Hamming and Hann are commonly used because you can avoid scalloping loss for 50% overlap. One gets a sufficiently favorable tradeoff between time resolution and frequency domain dynamic (side lobe level) range. There is also the fact that these windows are cited in papers that go back to the 70s.

You can also apply the widows with a 3 point convolution (3 coefficients instead of N) in the frequency domain, which might have had some advantages in hardware in the past when memory was a constraint.

Classical fixed windows are chosen according to characteristics like sidelobe level. They aren’t like adaptive windows that use the data itself.  It’s relatively pedestrian processing.

in summary, boxcar windows are for flat spectra.


Pitch Detection in the Browser Does anybody know a good library that does pitch detection in the browser? Anything in Javascript or Webassembly would suffice. My target is mainly singing voices. I already tried pitchfinder, but for me it reports frequencies that are a half tone to low. Are there any alternatives other than implementing an alorithm your self?

Try a samplerate of 48000 instead of 44100.  Or vice versa.


Real world application of signal sparsity? There are theories based on signal sparsity in frequency domain like  Compressive Sensing, Sparse FFT, etc. Throughout searching and studying papers I found out Cognitive Radio is a good example of application of Fourier sparse signals in real world. However, I am curious what are other applications in which the input signal is 1-D and is also sparse? or if there is any?
I feel there should be some important applications where the signal can be modeled as sum of a few sine functions, but I can't find any.

So, my question is in which applications the input signal  sparse in Fourier (only) and is 1-Dimensional?

I would appreciate a lot any suggestion.

Sparsity concept is extensively being used in computer vision and image processing. The Idea is that natural image can be pretty sparse when it is transformed to different bases. this bases can be predefined, e.g. FFT,DCT or can be learned from the image, e.g. sparse coding.

Here are a few well know example of algorithms that uses the sparsity assumption on images:

Jpeg image compression: Uses DCT transformation. assumes that high frequencies are low in natural image, so they are being removed in order to preserve space for compression. see more information on Wikepedia
Sparse coding: is a technique being used that assumes that an image has a lot of redundant information and because of that can be reconstructed from an overcomplete dictionary. each patch in the image should be sparse combination in the dictionary. see instroduction to sparse representation presentation
inverse problems in computer vision like recovering image embedded in noise and recovering sharpe image from blurred image are also using sparse assumption, e.g. wavelete shrinkage is being used for image denoising. See application to wavelet trasform in here