Logarithmic fourier transform(LFT) on audio signal I am trying to analyze the musics as possible as precisely.
Of course I tried FFT, but got some problems.

I found low frequencies have very low resolution than human's hearing. I tried very long time FFT to resolve this problem, but even analyzing with 8192 samples/s in 44100Hz sample rate(Means lack of time resolution), I got not enough resolution on low frequencies.

I found there are few solutions.

Firstly, a quadratic interpolation on FFT bins.
But it seems not a perfect way. Problems of this method are:

1. 'If i want to determine freqs between the freq bins, which three bins should I select to do an interpolation?'

2. 'Even I do this, there are no actual additional informations on result. I know interpolations are kind of tricky method.'

Secondly, extracting each freq bins with desired frequency, so I can extract the bins logarithmically.

But have a critical computational cost problem: (maybe over)N^2.

Thirdly, LFT(Logarithmic Fourier Transform).

This requiers logarithmically-spaced samples and gives me result exactly what I looking for with incredibly fast speed; https://stackoverflow.com/questions/1120422/is-there-an-fft-that-uses-a-logarithmic-division-of-frequency

But I have no idea with that algorithm. I tried to understand the paper and implement it, but it was impossible because of lack of my english and mathematical skills.

So, I need a help of implementation of LFT.

The simplest and most pragmatic solution is to use a normal FFT of a sufficiently large size that you get the required resolution at the lowest frequency of interest. E.g. if you want 1 Hz resolution at the lowest frequency of interest then you will need a 1 second FFT window, i.e. the FFT size would need to be equal the sample rate, e.g. 44100.

Note that even if you could implement a logarithmic FFT then it would still be bound by the laws of physics (information theory) and you would still need a similar length sample window - all you would gain would be convenience (not having to aggregate output bins) at the expense of performance.


Find transfer function given impulse I am given the following discrete time transfer function : $$G_d(z^{-1})=z^{-d}\frac{b_0+b_1z^{-1}}{1+a_1z^{-1}}$$ which has the following impulse response $$g_d[n]=\{0,1,-0.1,-0.05,...\}$$

How can I find the transfer function?($b_0$, $b_1$ etc)


Meaning of Hilbert Transform I understand the Fourier Transform which is a mathematical operation that lets you see the frequency content of a given signal. But now, in my comm. course, the professor introduced the Hilbert Transform.

I understand that it is somewhat linked to the frequency content given the fact that the Hilbert Transform is multiplying a FFT by $-j\operatorname{sign}(W(f))$ or convolving the time function with $1/\pi t$.

What is the meaning of the Hilbert transform? What information do we get by applying that transform to a given signal?

One application of the Hilbert Transform is to obtain a so-called Analytic Signal. For signal $s(t)$, its Hilbert Transform $\hat{s}(t)$ is defined as a composition:

$$s_A(t)=s(t)+j\hat{s}(t) $$

The Analytic Signal that we obtain is complex valued, therefore we can express it in exponential notation:

$$s_A(t)=A(t)e^{j\psi(t)}$$

where:

$A(t)$ is the instantaneous amplitude (envelope)

$\psi(t)$ is the instantaneous phase.

So how are these helpful?

The instantaneous amplitude can be useful in many cases (it is widely used for finding the envelope of simple harmonic signals). Here is an example for an impulse response:

Secondly, based on the phase, we can calculate the instantaneous frequency:

$$f(t)=\dfrac{1}{2\pi}\dfrac{d\psi}{dt}(t)$$

Which is again helpful in many applications, such as frequency detection of a sweeping tone, rotating engines, etc.

Other examples of usage include:

Sampling of narrowband signals in telecommunications (mostly using Hilbert filters).
Medical imaging.
Array processing for Direction of Arrival.
System response analysis.


How the Grayscale images are converted back to RGB? In case of RGB to grayscale conversion we take the average of R,G,B channels and convert them to Grayscale which has songle channel. How to convert back the Grayscale to RGB ??

It is not possible to get color back. A grayscale image can be converted to RGB image, by copying the grayscale data to all three RGB channels,but it will still be grayscale.

Note that in general, RGB is not simply averaged to grayscale, but a weighted average is used to match eye sensitivity.


Analog vs. digital communications in PCM and noise How do you determine if you are using analog or digital communications? For example when dealing with PCM and noise how can one notice any difference in the signals when using the two types of communications? Thanks

This is quite a broad topic.  Let's narrow the scope to toll quality telephone audio.  Further assume both analog and digital channels are bandlimited to the same passband.  I assume you don't want to figure out the difference between a pair of wires and a PCM channel (i.e. a pair of wires has a very slow decay in frequency response where as a digital channel is hard band limited).

Ok a for sort of apples to apples comparison of a roughly 4Khz voice band channel, you can look for sampling artifacts in the channel.  The sampling frequency for telephone audio is 8Khz.  The true passband for telephone quality audio is 200Hz-3600Hz (still roughly speaking here).

Look for quantization noise.  Quantization noise is not random and can be made nearly periodic by choosing the right test signal.  The PCM encoding for standard CODEC ditigal telephone channels uses a type of compression that follows a log scale (u-law or A-law).  The analog reconstruction filter does not completely eliminate the quantization noise.  Because of the compression scheme, the quantization steps increase with increasing signal amplitude.  If you create a test signal that is a ramp running from -Vp to +Vp (where Vp is the maximum signal amplitude the channel can handle), you will see discontinuities that are closely spaced near 0V and spaced farther appart near +Vp and -Vp.  Use a frequency of 1KHz so the signal frequency is harmonically related to the sampling frequency.  This will place the quantization errors nearly at the same location for each period of the ramp signal.

You can also use a spectrum analyzer and look for spurs in the frequency content of a test signal.  In this case, use a 1Khz sine wave.  Once again the quantization errors will remain at nearly the same position on each period of the sinewave.  The quantization errors will appear as distortion of the sine wave and show up as harmonics in the signal spectrum.  These harmonics will be greatest when the sampling frequency is harmonically related to the signal frequency so 1KHz and 2KHz are good choices.  As the test tone drifts away from intger divisions of the sampling rate, the quantization noise spreads and begins to look like a spectral contribution from random noise.  In telephone systems, a 1KHz test tone is generally not used for this reason.  The normal test tone frequency is 1004Hz.

Note that modern CODECS with high compression rates such as G.729 have a difficult time passing periodic signals.  These CODECs are designed in such a way as to interpret all signals as human voice or noise, so the reproduced periodic signal will be highly distored or modulated.


Complex output after inverse FFT of a real signal I have a real one dimensional signal s (light absorbance in a flow cell), which has significant noise and some periodic noise after performing a deconvolution of $S$ from $S_o$. Basically fft($S$) was divided by the fft($S_o$). Let us call the output as $R$ vector. In order to remove the periodic noise, I replaced the noise region in $R$ with 0+0i in a given range. The ifft($R$) is complex. When I plot t and ifft($R$), MATLAB warns that the imaginary parts of complex R arguments ignored  . The output looks as desired, free from noise.

Forum members mention that this can arise due to numerical precision, however the imaginary part in ifft($R$) is on the order of 0.0006i.

Is there is a better window (simple one) than this rectangular window in the frequency domain?

Even in this simpler example which one can test is a simpler code, generates a complex inverse output.

t = [0:1/80:60]'; % Time
x = sin(2 * pi * 2 * t) + sin(2 * pi * 0.05 * t); % Signal
X = fft(x);
X([1:6, 4780:4801]) = 0; % Removing low frequency 0.05 Hz
z = ifft(X); % **z turns out to be complex as well**
figure(1)
plot(t, z)
hold on
figure(2)
plot (abs(X))

Thanks.

Essentially, your code does not respect the inherent Hermitian symmetry of the output of the FFT. Here, your signal is odd-sized $2K+1$. Hence, this FFT yields a complex vector of coefficients $d$ (real) and $a_k$ (generally complex), arranged as:

$$ \left[d,a_1,a_2,\ldots,a_K,\overline{a_K},\ldots,,\overline{a_2},\overline{a_1} \right]$$

If you want to zero-out some coefficients, you should remove then by conjugate pairs (except for $d$), like in

$$ \left[0,0,a_2,\ldots,a_K,\overline{a_K},\ldots,,\overline{a_2},0 \right].$$

Note: for even sizes, a real coefficient $n_{K+1}$ is inserted between $a_K$ and $\overline{a_K}$.

Your line:

X([1:6, 4780:4801]) = 0

does not satisfy this symmetry. Try instead this one (only valid for odd length):

nACBin = 3;
X([1:nACBin+1, end-nACBin+1:end]) = 0; %

and you will get a much lower imaginary residue (here due to numerical issues). Here, it somehow works (see below) because of the very simple context (two sines with lucky interval length), but this is rarely recommended in general:

Especially with short 0 strides, this will convolve your data with a cardinal sine, both broad and wiggling.

As said in comments, you are considering prototyping a filter in the frequency domain using windows. This is a well-covered domain in DSP. However you can avoid those complexities by first looking at filtering in the time domain (Butterworth, Chebyshev, Cauer, Papoulis filters) that already have well-designed implementations. Other filters can be useful for peaks, like polynomial interpolators (Savitzky-Golay, for instance).

Side note: it probably would be more efficient (though a bit more complicated) to embed those frequency constraints into the deconvolution, instead of trying to remove the artifacts afterward.


What really means stochastic in field of signal processing I met two definitions of word stochastic, the first one (cited from wikipedia Stochastic)

  The word stochastic is an adjective in English that describes
  something that was randomly determined

The explanation of stochastic and deterministic what is used in textbooks really make sense according to definition above. The second definition is

  Stochastic means sample by sample

For example SGD - stochastic gradient descent.

This seems to be a little confusing for me. I do not really see any relation between random and step by step.

Do I miss something? Or it is really different meanings for the same word?

Well, getting a bit linguistic, according to the Oxford dictionary:

  stochastic (adj.): Having a random probability distribution or pattern that may be analysed statistically but may not be predicted precisely.

So the definition would be the first one (I don't know where you might have found the second one as you didn't put any source about it). Nevertheless, regarding your specific question, SGD is indeed stochastic. It is an iterative method, but this doesn't mean it is not stochastic. It is iterative and stochastic at the same time: the two concepts are not related to each other.

To put it clearer: if something is blue, couldn't it be big? Well... why not? That's exactly what you are asking here. As size isn't related to colour (absurd example that I think can help here), an approximation can be iterative and stochastic (such as SGD): being one thing doesn't mean it can't be the another one too.


To find the impulse repsonse using the difference equation A causal linear time-invariant filter has transfer function

a) Denote the input signal by x[n] and the output signal by y[n]. Find the difference
equations for the filter.

f) Find the impulse response of the system. (Eliminate j (the imaginary unit) from
your answer.)

I used the partial fractions method and obtained the following solution.

 However it was time consuming. Is there a way that I can simplify this difference equation to get the impulse response of the IIR filter to get the same solution. Any help will be appreciated(as I am not aware of solving by difference equations)

There is no real shortcut for computing the impulse response. Partial fractions is the standard way to do it. However, in the case of a second-order transfer function as in your example, the result can be found in tables of $\mathcal{Z}$-transform pairs. E.g., if you use the last two correspondence of this table, you can pretty quickly write down the result.

By the way, your result can't be correct (at least not for negative $n$), because you know that the system is causal, so its impulse response must be zero for $n<0$.


Raw speech data into spectrogram I am using Delphi trying to make a very simple speech recognition system. The system collects raw data well, but I have a hard time turning the data into a spectrogram:

The system collects an array with $1024$ integer values for each $100$ ms.

I try to FFT this array into $(R, I)$, and plot time on the $x$-axis and $\sqrt{R^2 + I^2}$ on the $y$-axis, but no frequency bands appears.

Does anyone know what I should do in stead? Source code in any language is very welcome.

In principle what you describe is reasonable. Here's an example in Python that implements your idea:

Fs = 44100
T = 100
t = np.arange(0, T, 1/Fs)
N = 1024
t = t[:N*int(len(t)/N)]

f = 0.7*t*t+2000

signal = np.sin(2*np.pi*f*t)  # generate some signal for our analysis

chunks = signal.reshape((-1, N))   # create chunks of the signal of length N
C = np.fft.fft(chunks, axis=1)     # perform FFT of each chunk
C = C[:, :N//2] # keep only positive frequencies

Spec = 10*np.log10(abs(C*C))      # draw the spectrum
plt.imshow(Spec.T);

I can only recommend to try out your DSP algorithms not in Delphi but in some scripting language like Python or Matlab, because you are much more flexible in your analysis and can try out new ideas much more quickly.

To improve the spectrum image, you can use a different window function (i.e. no rectangular window as I use here) and overlap the different chunks. Then, it would be closer to the short-time Fourier transform as noted by MBaz.


Bandpass fundamental and harmonics in Matlab I have a signal in time domain whose sample frequency is Fs=25600. I would like to remove from the fundamental F=285Hz and all its harmonics (2*F,3*F,etc). I tried to use the comb filter in Matlab using this code :

Fs=25600;
N=43;
BW=285;
Apass=200;
[b, a] = iircomb(N, BW/(Fs/2), Apass);
Hd= dfilt.df2(b, a);
x1 = filter(b, a, signal);

Here is the spectrum of the original signal over the frequency interval up to around 400Hz

Here is the result after applying the filter cited above:

I don't get the awaited result. Is there a way to accomplish this in Matlab?

You have a problem with the comb filter design. According to Matlab documentation the following line creates the filter you want:
[b,a] = iircomb(round(25600/285), 2*285/25600/35,'peak');

Looking at the frequency response of the designed filter:
Seems to be able to solve your problem...


Conventional channel coding/decoding I am trying to understand the algorithms of channel coding. Regarding the conventional channel which is, for me, the easiest one to implement, I have a question regarding the output of its decoder.
When using encoder whose rate is 1/2, the output of the decoder is half of its input, for example:
constlen = 7;  codegen = [171 133];
tblen = 32;     % traceback length
trellis = poly2trellis(constlen, codegen);

X = randi([0 1], 1024,1);        %The data to encdoe
Y = convenc(X, trellis);        %Performing of coding
Y_dec = vitdec(Y, trellis, tblen, 'cont', 'hard');  %The decoder

As you see, in last line, Y_dec is of length 1024. However, its input Y has a length of 2048. What I need is to show the whole corrected outputs, I mean I need to have the output Y_dec of length 2048 which includes the data and the redundancy. Is that feasible?
EDIT:
Following Dilip's feedback, I could get back it as below:
constlen = 7;  codegen = [171 133];
tblen = 32;     % traceback length
trellis = poly2trellis(constlen, codegen);

X = randi([0 1], 1024,1);        %The data to encdoe
Y = convenc(X, trellis);        %Performing of coding
Y_dec = vitdec(Y, trellis, tblen, 'cont', 'hard');  %The decoder
Y2 = convenc(Y_dec(1+tblen:end), trellis);
Y3 = Y(1:end-2*tblen) - Y2;     % This is Zero

But as you see, it means that a part of code word can be found which is Y2 = Y(1:end-2*tblen). however the code word size is $2048$, so is there a way to get all the codeword?

Yes, it is possible to determine the path that the Viterbi algorithm found through the trellis. Just apply the encoding algorithm to what you are calling Y_dec and you will get the corresponding codeword of length $2048$. You can then compare it to Y, the transmitted codeword, to see where the channel made errors.
Additional notes:
If the data to be transmitted 1024 information bits and the code to be used is a half-rate $(2,1)$ convolutional code with constraint length $7$, then the convolutional codeword produced by the convolutional encoder is $2060 = 2048+12 = (2\times 1024) + 2\times 6$ bits long where the extra $12$ bits are produced by the encoder as zeroes are pumped into the $6$-bit-long data buffer as the buffer empties out. I don't speak MATLABian gobbledygook at all and so don't know what Y_dec(1+tblen:end) means, or indeed why tblen is involved at all.  If  Y = convenc(X, trellis) produces a codeword vector Y of 2060 bits, then
Y_dec = vitdec(Y, trellis, tblen, 'cont', 'hard') should produce X again. I don't see where there was any noise involved and so the Viterbi decoder should just zoom through the trellis following the same path as the encoder so that Y_dec should equal X. If there is noise, then Y_dec should be the same as X in most of the 1024 bits except for some bursts of errors in that span. In this case, convenc(Y_dec, trellis) will not be the same as Y, but in the absence of any channel errors, Y_dec=X and  convenc(Y_dec, trellis)=convenc(X, trellis)=Y.


Designing digital low pass filter with low pass-band group delay I'm interested in specifying a digital low pass filter with low (and near constant) pass band group delay.  I tried python's spicy.signal.minimum_phase to transform a linear phase FIR low pass filter to minimum phase, but there was a steep price in stop-band rejection level.
Is there a "standard" approach to this problem?

Minimum phase filters will not give you a near constant group delay. You can design a non-linear phase FIR filter with a linear desired passband phase with a specified group delay that is smaller than the group delay of the corresponding linear phase filter. If you use a least-squares criterion, this is equivalent to solving a system of linear equations.
As an example, take an FIR lowpass with $61$ taps. The group delay of a linear-phase filter would be $30$ samples. Now we specify a group delay of $20$ samples instead:

w = pi * [linspace(0,.2,100), linspace(.25,1,375)];   % frequency vector
D = [exp(-1i*w(1:100)*20), zeros(1,375)];   % desired complex frequency response
W = ones(1,475);   % unity weighting
h = lslevin(61,w,D,W);   % least squares design

The function lslevin.m can be found here.
The result looks like this:

Note that the passband group delay is of course not exactly constant, but the desired group delay of $20$ samples is approximated quite accurately over a large part of the passband.
Trying to achieve a low group delay in the passband comes with a penalty in the magnitude response. The figure below shows the magnitude and group delay responses of $4$ length $61$ FIR lowpass filters with different desired group delays in the passband. One is a linear phase filter with a constant group delay of $30$ samples, and the other $3$ filters are non-linear phase low delay filters with specified passband group delays of $25$, $20$, and $15$ samples respectively. Clearly, the magnitude response deteriorates with decreasing desired passband group delay.


Power Spectral Density of Brownian Motion despite non-stationary Note: I originally asked this on Physics Stack Exchange but haven't attracted any interest there so I'm asking here where it may be more relevant.

A white noise process, $\xi(t)$ with delta correlated two-correlation function $\langle \xi(t_1)\xi(t_2)\rangle = \delta(t_1-t_2)$, is clearly stationary  and has a power spectral density which is the Fourier transform of the auto correlation function a la Wiener-Khintchine theorem.

$$
S_{\xi \xi}(f) = \int e^{i 2 \pi f t} \langle \xi(t)\xi(0) \rangle dt = 1
$$

We call such a process white noise because it has a flat power spectral density as indicated by the above equation.

I am concerned here with the power spectral density of a Brownian motion process or a Wiener process, $W(t)$ which is the integral (in some sense) of a white noise process.

$$
W(t) = \int_{t'=0}^t \xi(t) dt = \int_{t'=0}^t dW(t')
$$

See the Wikipedia article on Brownian noise. In that article it is pointed out that since Brownian motion is the integral of white noise we have (using the Fourier transform convention above)

$$
\mathcal{FT}[W](f) = \frac{1}{-i 2\pi f} \mathcal{FT}[\xi](t)
$$

Then the argument follows that the power spectral density is like the absolute value of the Fourier transform so we should have something like

$$
S_{WW}(f) = \frac{S_{\xi\xi}(f)}{|-i 2\pi f|^2} = \frac{1}{(2\pi)^2} \frac{1}{f^2}
$$

That is, the Brownian process has a $\frac{1}{f^2}$ spectrum, or Red or Brown noise.

We can also perform a direct computation of the power spectral density from the definition of power spectral density as the limit of the magnitude of the windowed Fourier transform.

$$
S_{WW}(f) = \lim_{\tau \to \infty} \frac{1}{\tau} \int_{t_1=0}^{\tau} \int_{t_2=0}^{\tau} e^{i 2\pi f (t_1-t_2)} \langle W(t_1)W(t_2) \rangle dt_1 dt_2
$$

It can be shown for a Wiener process that the two-time correlation function is $\langle W(t_1) W(t_2) \rangle = \text{min}(t_1,t_2)$. We can write the integral as

$$
\frac{1}{\tau}\int_{t_1=0}^{\tau}\int_{t_2=0}^{t_1} e^{i 2\pi (t_1-t_2)} t_2 dt_1dt_2 + \frac{1}{\tau}\int_{t_1=0}^{\tau}\int_{t_2=t_1}^{\tau} e^{i 2\pi (t_1-t_2)} t_1 dt_1dt_2
$$

Mathematica (or a diligent calculus student) evaluates this integral to

$$
\frac{2\pi f\tau - \sin(2\pi f \tau)}{4 \pi^3 f^3 \tau}
$$

In the limit of $\tau \to \infty$ the second term vanishes and we are left with

$$
S_{WW}(f) = \frac{1}{2 \pi^2} \frac{1}{f^2}
$$

This is missing a factor of 2 compared to the formula above. I think the reason is related to the fact that I've only considered positive time for the windowed Fourier transform in the second derivation.

Here is my question
I have given two sketches that show Brownian motion has power spectral density $\sim \frac{1}{f^2}$. However, I have also stated the two-time correlation function for Brownian motion is $\langle W(t_1)W(t_2) \rangle = \text{min}(t_1,t_2)$. This shows that Brownian motion is clearly not a stationary process. I know that the Wiener-Khintchine theorem tells us that if a process is stationary then the power spectral density is given by the Fourier transform of the two-time correlation function. But I also had the impression that non-stationary processes didn't have well-defined power spectral densities, or that you could only define a time-windowed power spectral density. However, Brownian motion seems to be an exception that rule. It seems to be an example of a process which is non-stationary (so you can't use the Wiener-Khintchine theorem, for example) yet still has a well-defined power spectral density. Am I correct in this characterization? Can anyone give me more information about non-stationary processes that still have power spectral densities?

Please see N.J. Kasdin, "Discrete Simulation of Colored Noise and Stochastic Processes and $1/f^\alpha$ Power Law Noise Generation", Proceedings of the IEEE, 83(5), (1995), 802-827. On page 806, Kasdin discusses the factor of 2 issue and states that "It is this author's belief that the factor of two in (24) is an error due to a broad-band bias introduced by the rectangular window." With a Hanning window (with 8/3 normalization), Kasdin gets an expression (equation 25) that converges to $1/\omega^2$ for large T. He also says "Since the formal spectral density of Brownian motion doesn't exist, the better question is which provides the most accurate information about the second order properties and possible linear models (see Section II-D) for Brownian motion?" Hope this helps a bit!


Impulse Invariant method for digital filter design One of the known methods for discretizing analog filters is impulse response invariant. We get the impulse response in time domain, discretize it and then get the Z transform.
What I am trying to understand is why the freq response of the resulting digital filter has a freq response magnitude scaled by (1/T) T:sampling time?
Matlab, using c2d command, modifies it by multiplying by T so that the freq response is similar to the analog filter, but this is not the result of the Z transform I described earlier.

This is just as it turns out when you do the math. The discrete-time Fourier transform (DTFT) of the sampled continuous-time impulse response $h(t)$ is
$$H_d(e^{j\omega T})=\sum_nh(nT)e^{-jn\omega T}\tag{1}$$
With
$$h(nT)e^{-jn\omega T}=\int_{-\infty}^{\infty}h(t)e^{-j\omega t}\delta(t-nT)dt\tag{2}$$
this can be written as
$$\begin{align}H_d(e^{j\omega T})&=\sum_n\int_{-\infty}^{\infty}h(t)e^{-j\omega t}\delta(t-nT)dt\\&=\int_{-\infty}^{\infty}\left[h(t)\sum_n\delta(t-nT)\right]e^{-j\omega t}dt\\&=\mathcal{F}\left\{h(t)\sum_n\delta(t-nT)\right\}\\&=\frac{1}{2\pi}H(\omega)\star\frac{2\pi}{T}\sum_k\delta\left(\omega-\frac{2\pi k}{T}\right)\\&=\frac{1}{T}\sum_kH\left(\omega-\frac{2\pi k}{T}\right)\tag{3}\end{align}$$
where $H(\omega)$ is the Fourier transform of $h(t)$, and $\star$ denotes convolution. From $(3)$ we see that the DTFT of the sampled impulse response equals the sum of shifted spectra of $h(t)$, scaled by $1/T$.
If we assume that $H(\omega)$ is approximately band-limited and that $T$ is chosen sufficiently small such that aliasing becomes negligible, we obtain the approximation
$$H_d(e^{j\omega T})\approx\frac{1}{T}H(\omega),\qquad |\omega|<\frac{\pi}{T}\tag{4}$$
For the step-invariance method, we use samples of the step response instead of samples of the impulse response, and we obtain a relation analogous to $(3)$ between the DTFT $G_d(e^{j\omega T})$ of the step response of the discrete-time system, and the Fourier transform $G(\omega)$ of the continuous-time step response:
$$G_d(e^{j\omega T})=\frac{1}{T}\sum_kG\left(\omega-\frac{2\pi k}{T}\right)\tag{5}$$
In order to obtain the frequency response $H_d(e^{j\omega T})$ we multiply $(5)$ by $1-e^{-j\omega T}$, because the impulse response is obtained by computing a first-order difference of the step response:
$$H_d(e^{j\omega T})=\left(1-e^{-j\omega T}\right)G_d(e^{j\omega T})=\frac{1-e^{-j\omega T}}{T}\sum_kG\left(\omega-\frac{2\pi k}{T}\right)\tag{6}$$
For frequencies that are small compared to the sampling frequency, i.e., for $|\omega T|\ll 1$ we obtain from $(6)$
$$\begin{align}H_d(e^{j\omega T})&\approx\frac{1-(1-j\omega T)}{T}\sum_kG\left(\omega-\frac{2\pi k}{T}\right)\\&=j\omega \sum_kG\left(\omega-\frac{2\pi k}{T}\right)\tag{7}\end{align}$$
If we again assume that aliasing can be neglected, we arrive at
$$H_d(e^{j\omega T})\approx j\omega G(\omega)=H(\omega),\qquad |\omega|<\frac{\pi}{T}\tag{8}$$
From $(8)$ we see that, unlike for the impulse invariance method, the step invariance method doesn't involve a scaling of the continuous-time frequency response.


Software PLL tracking of carrier frequency in bandlimited transmission My understanding is that a PLL(either in hardware or software) can recover the carrier frequency(or be used to correct the frequency offset) of a system. My question is how does this work when the transmitter is sending a pulse-shaped signal?

For a an approximate sinusoid the operation of the PLL makes sense but how does tracking work for carriers transmitting data?

As an example, suppose you're receiving a raised-cosine pulse, what do you do?

Here is a link to an example carrier tracking implementation for a QPSK and QAM system, applicable to full-response raised cosine pulse shaping. This post also has links to many other concerns with demodulating the waveform, including timing recovery and setting the parameters in the tracking loop implementation.

Recovering signal for psk

A x4 carrier recovery approach may be the simplest way to see how the pulse shaping effects the carrier recovery. For QPSK without pulse shaping you can recover the carrier by raising the waveform to the fourth power, and then dividing the resulting 4x carrier by four to complete the recovery. This is because phases add when you multiply, so if we have four possible phases in our modulated signal:

$$e^{(j\omega t + 0)}$$
$$e^{(j\omega t + \pi/2)}$$
$$e^{(j\omega t + \pi)}$$
$$e^{(j\omega t + 3\pi/2)}$$

Raising to the 4th power results in the following for all the cases above.
$$e^{(j4\omega t + 0)}$$

If we observe our waveform over time and it instantly goes from one phase location to the other (no pulse shaping, infinite bandwidth), then the 4th power result will be a pure tone at 4x the carrier frequency with no other components.

However if we slowly transition from one phase to the other, we will still get a pure tone at 4x the carrier frequency, but there will be significant energy in sidebands close to this frequency as well. The 4x tone will still be the strongest in vicinity of this 4x signal and is easily tracked with a PLL to complete recovery. (This is evidenced by a positive mean in the baseband equivalent analytic signal).  It is because this signal is the strongest that a PLL can track it (typically if the signal being tracked is +6 dB higher than the sidebands a PLL can lock to it), and it is because of the additional sideband energy that we need to the PLL to complete the recovery, effectively acting as narrow bandpass filter to provide a clean carrier.

I added some graphics below for the simpler case of BPSK with and without pulse-shaping that will hopefully provide an intuitive answer. Personally I find it much simpler to think and work with describing frequency using the complex exponential ($e^{j\omega t})$ rather than with cosines and sines, and for the same reason I avoid including an actual carrier but model at the baseband equivalent analytic signal. However I think in this case showing an actual sinusoidal carrier will more immediately get across the point I am trying to make. For the same reason, this is very illustrative in showing how the pulse shaping affects a tracking loop but it is not how I would typically do carrier recovery (see the link above for a typical software recovery loop I would implement).

This first image shows a sinusoidal carrier with an example (1 0 1) data pattern. Below the data pattern is shown the result of BPSK modulation with no pulse shaping applied (infinite bandwidth available).  Observe in the third plot showing the modulated signal, that the phase is changing from 0 to $\pi$ radians as the data transitions. This abrupt change back and forth of the phase would prohibit a simple phase tracking PLL from tracking the phase: if it could even lock onto one phase state (because the data doesn't transition for some time for example), once the data transition occurs the very significant phase error that would occur would cause a simple tracking loop to break lock (unless the phase detector had 0°/180° ambiguity -- hint, hint at another recovery approach!). Another way of describing this, with a random data pattern that is 50% ones and 50% zeros, the carrier would be completely suppressed and therefore does not exist, so there is no signal at that frequency for any loop to lock onto. Notice too that in the example plot I give, the data transitions are not commensurate with the carrier. I did that purposely just to show that regardless of the data rate and where transitions occur, squaring this signal will result in a pure tone (when no pulse shaping or other bandwidth limiting is used) at exactly twice the rate of the carrier.

The fourth plot is simply the waveform in the third plot squared with no other filtering applied. This waveform is completely in phase with the original carrier, at exactly twice the frequency.  In this case we could simply pass this signal into a frequency divide by two function (quite simple) and we would completely recover the carrier. No PLL is needed, but if there were any other noise sources present, a PLL would help increase the SNR of our recovered carrier in the presence of noise. However our waveform itself and the recovery process, in this case, contributes no additional noise or distortion.

What is occurring by squaring the signal is we are doubling the frequency (and the phase). This is clear from the trigonometric relationship:

$$cos(\alpha)*cos(\beta) = \frac{1}{2}cos(\alpha+\beta)+ \frac{1}{2}cos(\alpha-\beta)$$

$$cos(2\pi f_c t + \phi)cos(2\pi f_c t + \phi) = \frac{1}{2}cos(4\pi f_c t + 2 \phi)+ \frac{1}{2}cos(0)$$

$$ = \frac{1}{2}cos(4\pi f_c t + 2 \phi)+ \frac{1}{2}$$

So if the phase is only in two states: 0 (0°) and $\pi$ (180°), doubling this results in 0 and 2$\pi$ which is also 0.

Now observe what occurs in the same situation with pulse shaping as shown in the plot below. Here the only difference is a pulse shape has been added to each symbol. Notice that the result is to only add an amplitude variation (the pulse shape) to the squared signal but it has no effect on the phase. This is an important observation as to why a normal tracking loop could still track the phase in this situation, and also illustrates why the loop that is used should be able to "fly-wheel" in the absence of updates (as most tracking loops in both carrier and timing recovery). The phase in the doubled signal is in phase with the carrier at exactly twice the frequency but the amplitude is going up and down. This does of course impact the SNR of the recovered signal, as there are fewer phase error updates for tracking but unlike the case of the modulated signal itself, the phase does not abruptly change in a manner that would be disastrous for the loop to track. This also shows the importance of using a phase tracking loop in this case, which will ride through the amplitude variations and create a clean doubled carrier with no amplitude variation.

This is perhaps clearer by comparing the two final cases of the multiplied outputs directly as shown below. Here we clearly see that the pulse shaping only affects the amplitude and not the phase, hence a phase tracking loop should be able to track and recreate a cleaner replica of the doubled carrier signal, which we would then follow with a divide by two to completely recover the carrier.

The plots above should provide a satisfactory intuitive explanation in the effects of pulse shaping on carrier tracking loops, but not to suggest the favored approach to doing carrier recovery. In a digital or software implementation instead of squaring, PLL and divider, I would suggest decision directed and other similar digital Costas-Loop approaches as I have shown in the link above.

Another intuitive explanation is given by observing the above signals in the frequency domain, and noting how the carrier is suppressed prior to squaring and then recreated at twice the carrier frequency after squaring (or raising to the fourth power for QPSK, to the eight power for 8-PSK etc).


FFT of color images incorporated into an Object Recognition method? Is there any instance where Fourier transforms of color images have been used in conjunction with other object recognition method?

Any instance of usage of Fourier transforms in color images? I think there is great potential in the usage of Fourier transforms of color images.

Okay I found some ideas about using Fourier Transforms in Object Recognition. What I am unable to understand is however, its downsides. Why was the method discarded.

For those who are interested in reading about how Fourier Transforms were applied in Object Recognition , look for Fourier Mellin Transforms.

The Fourier-Mellin transform is a useful mathematical tool for image recognition because its resulting spectrum is invariant in rotation, translation and scale. The Fourier Transform itself (FT) is translation invariant and its conversion to log-polar coordinates converts the scale and rotation differences to vertical and horizontal offsets that can be measured. A second FFT, called the Mellin transform (MT) gives a transform-space image that is invariant to translation, rotation and scale.

The steps:
The Fourier Transform

The Discrete Fourier Transform (DFT) is given by the following expression:

This is often computed using the Fast Fourier Transform algorithm to speed up the time needed for the calculation.

Cartesian to Log-Polar conversion
The FFT is projected onto the log-polar plane by the coordinate transform shown below.

For the conversion from Cartesian coordinates to Log-Polar coordinates the following equation is true:

$$r = \sqrt{x^2+y^2}$$

The origin $(m_o, n_o)$ should be at the center of the image matrix to ensure the maximum number of pixels is included. If the image consists of a square $N \times N$ matrix then the coordinates of the origin are:

The maximum sampling radius for the conversion can now be calculated as:

If an inscribed circle is chosen as the conversion boundary, some pixels that lie outside the circle will be ignored. If a circumscribed circle is chosen, all pixels will be taken in account, but some invalid pixels will also be included (pixels falling inside the circle but outside of the image matrix).

Since the pixels in Cartesian coordinates cannot be mapped one-to-one onto pixels in the Log-Polar coordinate space, an average of the surrounding pixels needs to be calculated. The standard methods to do this includes nearest neighbor, bilinear and bicubic resampling.

The relationship between the polar coordinates ($\rho$, $\theta$) used to sample the input image and the polar coordinates of the log-polar image (r , $\theta$) can be described by:

$$(\rho,\theta)    =  (e^r, \theta) $$

To map the input image pixels imagein(xi,yi) onto the output image pixels imageout(rm, $\theta$m), the co-ordinates xi, yi are computed using:

where $$(\rho,\theta)    =  (e^r, \theta) $$

The next step is to get a transform-space image that is a rotation and scale invariant representation of the original image. This is done using the Mellin transform which can be expressed as:

Converting to polar coordinates, we have:

Now using $$(\rho,\theta)    =  (e^r, \theta) $$

we get:

The Fourier Mellin transform has been used for Image Recognition in R. Moller, H. Salguero, E. Salguero: “Image recognition using the Fourier-Mellin Transform”, LIPSE-SEPI-ESIME-IPN, Mexico


What is the physical significance of the Fourier coefficients in the DFT of audio signals and how can they be best displayed in a spectrogram? I am trying to display the spectrogram of an input audio signal (a newcomer to all things spectral here). The amplitudes of my audio signal can be assumed to be normalized between -1 and +1.

Suppose for the frequency $f$ in the DFT's output, the corresponding coefficient is $c$, in general a complex number. I would like to display the "amplitude" information contained in $c$ as a corresponding colored region on the spectrogram. Ideally to plot the spectrum, I would need a number that ranges from 0 to 1. 0 will represent the lowest amplitude I will record, and 1 the highest. There are a few ways I can think of to achieve this:

Since the $|c| \le N$, where $N$ is the number of samples, I can normalize the magnitude of $c$ with respect to $N$ ($v = |c|/N$) and display the resulting $v$ as a spectral colour.
Same as 1, but normalize $|c|$ wrt to $M$, the maximum magnitude of all Fourier coefficients. But this has the disadvantage of having to scan the entire data before the spectrum can be displayed, and will be a problem for large files. Also, this means that the spectrum display is relative to the intensity levels present in the current audio signal, and not valid in an "objective" sense.
Convert $|c|$ into decibels, say $d = 20 \log_{10}(|c|/b)$, where $b$ is a base magnitude. I will then normalize $d$ from say $l$ (where $l$ is a lowest dB value I will register) to $u = 20 \log_{10}(N/b)$. I can then use $d$ in displaying the spectrum. But in this case, what should $b$ be?

My question: what is the exact physical significance of the magnitudes of the Fourier coefficients in the DFT of audio signals and how can it be best displayed? I have tried all the above approaches, but am wondering which of them is the most proper or if I am missing any others.

This is a two year old question, but I think I add a little bit for those who might find it doing searches.

When you are talking about a spectrogram, there are really two steps involved that have to be considered.  First is converting from the signal to DFT bin values, and the second is displaying those bin values in a pixel display.

Ideally, a spectogram of a pure tone of increasing frequency should make a solid, non varying, upward sloping curve.  Inherently DFTs don't work this way.  When you are at a whole integer frequency (in terms of cycles per frame), there is a single bin value, the rest are zero.  When you are halfway between integer frequencies, the bin values are spread out which is commonly known as leakage.  My recommendation is to use a VonHann window on your signal before you do the DFT.  Depending on your scaling, this is the equivalent of subtracting the average of the two neighbor bin values from each bin value.  This has the effect of "tightening up" the leakage, and "spreading out" the sharp peaks.  This makes the display look a lot more consistent independent of the value of the frequency compared to the frame length, which is what you want.

The phase values in the DFT are framing dependent, and therefore don't help much with the display.  Instead you should concentrate on making an intensity display based on the magnitude of the bins.  You are not limited to a linear mapping, or a log mapping.  You are going to want to use a lookup table anyway to speed up processing, so you can pretty much define an monotonically increasing function that works the way you want.  Use can use a monochromatic color scale, or changing color one.  The standard green for low, yellow for mid level, and red for intense is easily comprehended.

The second choice you need to make is the vertical axis scale.  The two logical choices are linear, which the DFT is inherently. or logarithmic corresponding to musical scales.  In either case, you are going to need to map a bin, or a set of bins, to a pixel, or a set of pixels, depending on the magnification at that part of the scale.  Again, in practice, you will build lookup tables for this mapping.

The most proper choices in the above considerations are the ones that give you the visual results you desire.  You are not going to be making calculations based on the display values so there is not problem with "distorting" the values (i.e. non-linear mappings) to give better visual results.

There actually is a physical significance of the values of the DFT coefficients.  Discovering this was a major aha moment in my understanding of the DFT.  A DFT bin value is actually a weighted average of a set of Roots of Unity, where the signal values are the weights and they are "wrapped" around the unit circle.  The bin index is the step size in roots, and thus also then number of times the function is stretched and wrapped around the circle.  For instance, suppose you have a signal with five bumps in it.  When is stretched five times and wrapped around the circle five times, all the bumps align and throw the center of mass in that direction.  On the other hand, when it is strectched by a factor of four, and wrapped around four times, the bumps are evenly spread and cancel each other out.  This is how a DFT "plucks" out the frequencies that are in a signal.

I have written a blog article "DFT Graphical Interpretation: Centroids of Weighted Roots of Unity" to explain this.

This is long enough, hope it helps some one.

Ced


Noise reduction from very noisy audio I am trying to write an algorithm that would automatically segment a piece of audio with bird calls recordings. My input data are 1 minute-long wave files and on the output I would like to get separate calls for further analysis. Problem is that signal-to-noise ratio is quite terrible due to environmental conditions and poor quality of a microphone (mono, 8 kHz sampling).

I would be most grateful for any advice on how to proceed further with noise reduction.

Here is an example of my input, one minute audio recording in wave format: http://goo.gl/16fG8P

This is how the signal looks like:

Band-pass filtering, in which I am keeping only anything in between 1500 - 2500 Hz, does improve situation, but still it is far from expectations. In this spectrum still a lot of noise is present.

I have also plotted long-term (over 32-sample intervals) average energy and removed some clicks from it. Here is the result:

With all the remaining noise I have to set a very low threshold to the onset detection algorithm to pick last 10 seconds of bird calls. Problem is if I tweak it in such a way then in next recording I can get load of false positives.

Moving average filter helps a bit with wind noise. Any other ideas? I was thinking of "Spectral Subtraction", but here it seems to me I have sort of chicken and egg problem - to find noise-only area I have to segment the audio and to segment the audio I need to remove the noise. Do you know of any libraries that have this algorithm or some implementations in pseudo-code? Methinks Audacity uses such a method to remove noise. It is very effective, but it is left to the user to mark noise-only area.

I am writing in Python and it is a free, open-source project.

Thanks for reading!

At the end what has proven to be the best solution was onset detection based on either high frequency or energy content. Before it could work I had to use high-pass filter to cut out first 1 kHz, since it contained too much noise.

Once I had noise-only area I could use its profile to reduce noise from rest of the sample.

One library I found particularly useful was Aubio. It has a good set of examples and provides a lot of algorithms to choose from for onset detection.


Are circular buffers relevant to transposed direct form IIR filters? I'm currently working on a project where I'm programming IIR digital filters. I have already finished the direct form 1 and & 2 implementations.

When doing direct form 1 and 2, I was able to leverage the use of circular buffers for storing data. This allowed me to easily insert new values into a tap line and save the cost of shifting data in memory.

However, I've just started working on the transposed direct form 1 filter and am having a difficult time in identifying whether or not a circular buffer is needed. Unlike the direct forms 1 and 2 where the memory was simply shifted, it seems like the transposed forms require all the previous values to be rewritten.
Essentially, the block diagram appears that every memory state will change depending on some current value.

Referring to the transposed direct form 1 diagram illustrated here, it shows that each 'state' will be dependent on a previous state as well as some instanteanous value. The instantaneous value in this case is v[n].

So my questions:

Do/can transposed direct form filters benefit from circular buffers?
If not, does this mean that transposed realizations are computationally 'more expsensive' due to the fact that every state (memory) has to be recalculated for every filter update?

Circular buffers are only useful if the contents of the delay elements are not changed when shifted. So you're right that circular buffers cannot be used for the transposed structures. However, this doesn't necessarily mean that the transposed structures are less efficient than the direct forms, it just means that the data aren't just shifted through a delay line. The complexity for computing the new state vector in the transposed structure is comparable to the complexity for computing the output from the states in the direct form.


FFT analysis for Vibration Signal I am working on balancing an air-spindle. For the unbalance analysis I use an accelerometer (NI device). I have the voltage signal from the accelerometer corresponding to the vibration of spindle at a particular frequency (rpm) saved in an excel file. To analyze the unbalance from this vibration signal I use the fft function in MATLAB. My data is sampled at a sampling frequency of $100,000\textrm{ Hz}$. I am using the same example code given in the fft documentation. The result of fft is as shown.

.

The present signal consists of signals of various frequency. How can I filter this result in order to get only the signal of frequency equal to my rotational frequency ($40\textrm{ Hz}$). I would need this signal to calculate the vibration phase and magnitude by comparing it with the tachometer signal. My MATLAB code looks like this:

%For vibration analysis of signal without any trial mass.

filename = '2400RPM.xlsx';
sheet = 1;
xlRange = 'C40:C516039';
x = xlsread(filename,sheet,xlRange);
T=1/100000;
L=length(x);
Fs=1/T;
t=(0:L-1)*T;
Y = fft(x);
mag1 = abs(Y/L);
mag = mag1(1:L/2+1);
mag(2:end-1) = 2*mag(2:end-1);
ph1 = rad2deg(Y/L);
ph = ph1(1:L/2+1);
ph(2:end-1) = 2*ph(2:end-1);
f=Fs*(0:(L/2))/L;

%PLOTTING RESULTS
%--------------------------------------

subplot(2,2,[1,2])
plot(t,x);
title('Vibration Signal: 2400RPM');
xlabel('Time (seconds)');
ylabel('Amplitude (voltage)');

subplot(2,2,3)
plot(f,mag);
title('Magnituge Plot');
xlabel('Frequency (Hz)');
ylabel('Amplitude');

subplot(2,2,4)
plot(f,ph);
title('Phase Plot');
xlabel('Frequency (Hz)');
ylabel('Phase (degree)');

Any help is much appreciated.

NOTE: (Updated noting the sampling rate and desired frequency) Your post describes extracting a 40 Hz signal at a 100,000 Hz sampling rate. Your filtering implementation will be significantly easier (and more effective) if you apply course low pass filtering and decimation to a lower rate before applying the 2nd order tuned resonator derived from an exponential averager as shown in detail below. Please review this post Fast Integer 8 Hz 2nd Order LP for Microcontroller. I suggest decimating to a 200 Hz sampling rate which is a decimate by 500 (in stages such as 5-10-10).

Once decimated, consider using an exponential averager implementation, following this post for converting the lowpass filter to a bandpass, tunable to any frequency including your desired 40Hz:

Transfer function of second order notch filter

The difference between the notch filter in the post above and an exponential averager, is that the exponential averager has the following transfer function:

$$H(z) = \frac{\alpha z}{z+\alpha-1}$$

I think you will see the similarity and be able to follow the process in the linked post above. A diagram of the implementation as an "Exponential Averager" or low pass filter is shown below. This would create a maximum at DC and would need to be converted to a bandpass implementation similar to what was done in the linked post.

UPDATE: In pursuing the transformation myself I came up with the following result:

$$ H(z) = K\frac{z^2}{z^2-2\beta \cos\omega_n z+\beta^2} $$

Where:

(update, adding the closed-form equation for K)

K is a normalizing constant that is dependent on $\omega_n$ and $\alpha$ given as:

$$K= \alpha \sqrt{\beta^2-2\beta cos(2 \omega_n)+1}$$

$\omega_n$: center normalized radian frequency, between 0 and $2\pi$

$\alpha$: bandwidth constant, bandwidth gets tighter as $\alpha$ approaches 0.

$\beta = 1- \alpha$

The resulting implementation as a tuned 2nd order resonator is shown below.

Below are two examples with different bandwidths at center frequency $\omega_n=\pi/3$. The first plot is with $\alpha =0.1$ and the second plot is with $\alpha = 0.001$.

Note that the very long delay this second filter would have as evidenced by the very steep phase slope through the passband (group delay = $d\phi/d\omega$), as required for any tight bandwidth filter. This would also mean a relatively long convergence time. Given the ability to tune bandwidth and center frequency independently, you can envision how this can enable an optimized implementation for certain applications that is initially wide band for fast convergence while the bandwidth tightens over time to maximize SNR for the signal of interest. This structure can also be used as a bandpass loop filter or in a bandpass Sigma Delta implementation for bandpass sampling.

Below is an equivalent form of the exponential averager that only requires one multiplier as Rick Lyons has shown at the site linked below. Here I have chosen $\alpha$ to be similar to what Rick had done, in that $\alpha$ is a small number approaching 0 to make the bandwidth tighter (this is in contrast to my other post linked above for the notch filter implementation where I used $\alpha$ approaching 1 to be a tighter bandwidth, so swapping $\alpha$ with $1-\alpha$ between these two posts. I think the formula is cleaner if you actually use $\alpha$ as the number approaching 1, in that it eliminates my need to use $\beta$ as I have done above without the formulas appearing too cumbersome.

https://www.dsprelated.com/showarticle/182.php


Interpretation of Convergence plot in terms of mean square error I am trying to learn various learning algorithms used in training of neural network. When plotting the Mean Square Error in dB with respect to the number of iterations, if the points are around 0 dB then can we say that the curve converges or does the points need to go to minus dB? If the curve goes towards minus dB with the increase in iterations, then what can be inferred? How good or bad is the technique?

0 dB corresponds to a MSE of 1. So, if your method gives MSE that approaches 0 dB as the number of iterations increases, then your method has an error floor; it cannot drive the MSE to any value smaller than 1 no matter how many iterations it takes (or how many data points you have) etc. A good method should drive the MSE to 0 (corresponding to $-\infty$ dB), that is, the MSE should decrease continually as the number of iterations increases, and not saturate at a value of 1 (0 dB). We can then argue as to whether the approach to 0 MSE ($-\infty$ dB) is exponential (MSE $= O(e^{-n})$) or polynomial (MSE $= O(n^{-k})$) including as special cases quadratic (MSE $=O(n^{-2})$) and
linear (MSE $= O(n^{-1})$), or sublinear
(e.g., MSE $= O(n^{-1/2})$ or MSE $= O((\ln n)^{-1})$ or similar)
where $n$ is the number of iterations (or number of data points).
Naturally, the faster the decrease, the better (and usually
more expensive) the method.

On the other hand, if the MSE of the best currently
known method is approaching
an error floor of 0 dB, there is a fantastic research opportunity (and
possibly a get-rich-quick opportunity too): devise a method that
drives the MSE to 0 instead of 1 and gain instant fame
(and if you patent the method and found a start-up company
to further develop it, possibly fortune RSN as well).


Are Chebyshev Type 1 filters very unstable at very high orders? I am currently experimenting with a digital implementation of a Chebyshev Type 1 filter of mine, for audio purposes. I have written the code myself so I was wondering if the behavior I've noticed is due to bad code or something else. My implementation works perfectly fine up to an order of 30. I know that this is ridiculous but that's not the point, it's just for educational purposes.
Anyway, when using a filter of order 30 or higher part of the frequency response that should be flat has resonant peaks. For a lowpass filter, it meant boosted low frequencies for a highpass filter it resulted in boosted low-mids. The highpass filter did not behave this way when the cutoff was above the resonant peaks, by the way, and changing the ripple didn't change the trend.
Is this normal and a result of the filter getting more and more unstable? Or does this indicate that there is some part of the code not working properly?
I am using double precision where I can, but the framework I am working in wants the sample values with 32-bit precision so I can only do the filter coefficients and filter memory with 64-bit precision and have to reduce precision when writing to the audio buffer. Maybe that's an issue?

Edit: To add to my OP, I am using cascaded biquads to realize the filter.

With infinite-precision arithmetic they will be stable. However, since you don't have infinite-precision arithmetic, you will have quantization issues even if you use 64-bit precision. These quantization issues can make your filter unstable. Even if your filter is stable, perhaps you will not get the frequency response wanted because of these quantization issues.

First solution : split your order-30 IIR filter into 15 order-2 IIR filters and cascade them. This should fix your problem.

https://www.dsprelated.com/freebooks/filters/Series_Second_Order_Sections.html
https://www.mathworks.com/help/signal/ref/tf2sos.html

Alternatively : Instead of cascading, you can parallelize the order-2 IIR filters.

https://www.dsprelated.com/freebooks/filters/Parallel_First_and_or_Second_Order.html

Finally, do you really need an order-30 chebyshev filter? I find it hard to believe to you really more than an order-8 or order-10 IIR filter.

Edit : You edited your post to mention that you already use biquads. Have you tried an simpler filter, like an order-4 IIR? Make sure it works for an order-4 IIR filter (i.e. 2 cascaded biquads) and then try increasing the order. Also, are you sure that your factorization algorithm is correct? Maybe, there's a simple bug hidden somewhere and it will be easier to spot it with a reduced-order filter.


Object identification based on both shape and intensity -- i.e., identify the stars and their diffraction spikes at the pixel level Here's an example image to make the problem very clear:

I've used OpenCV's blob_log() as well as the rank.otsu threshold approaches. The problem with the blob detection is that the blob radii are quantized, e.g., 3, 7, 11, etc. The otsu local adaptive threshold works fairly well but needs severe hyperparameter tuning -- in other words, it works in some locations but not others.
I don't know how I can ask the following question: "find all circular blobs in the image. If the blob is large (meaning it is overexposed), expand the search to include diffraction spikes (at known fixed angles) plus the Airy disk."
My goal is to remove the stars and their optical effects from the image, then use an inpainting algorithm to fill in the missing pixels. I have several inpainting candidates to try.
There is noise in the image, and I'm fine with missing any stars that are indistinguishable from noise.
This is one channel -- the other image channels do add additional information, and the actual removal of stars/effects will be based on the logical or of all three analyzed channels.
I am aware that there are other methods for star removal that pop up in google search. The best automatic one I've used is starnet, e.g. in https://github.com/nekitmm/starnet . The issue I have with it is it doesn't separate out the star mask but goes ahead and does the inpainting -- unfortunately, very poorly in large regions. And it doesn't handle spikes well, per the author's commentary and my observations.
I wanted to avoid a machine learning approach. Instead, I wanted to solve the problem directly from first principles. It seems feasible to build a structural model of the star/spike/airy disk with some effort(e.g., as central disk at max intensity becomes larger, spikes grow larger and so does Airy disk.) But how do I match this model to the objects in the image?
Thanks in advance.

The desired output is a 1 bit mask image. Bit set to 1 if this is a star pixel, 0 otherwise. I don't have a problem converting a list of (y,x,r) values to a bit mask, and already have code to do that.
The blob detectors I have used (a) don't compute accurate blob sizes (b) generate blobs for some of the noisy areas.

Regarding finding the isolated blobs being a snap, give me an idea of what you consider is a "snap." My very first attempt looked for local maxima and attempted to grow a circular region. It was very slow and had problems with noise.

The region-growing approaches all failed miserably. A different 2d filtering approach turned out to be more promising. The basic idea is to:

Capture all of the overly bright areas via thresholding, e.g., against around 248 in the bottom image. That gives the first mask = image > 248, e.g.
Apply an emboss filter at 45 and 135 degree angles (filter kernels are [-1 -1 0][-1 0 1][0 1 1] plus 128 for the 45 degree angle.)
Threshold values more than 45 away from grey (128). mask = (embossed < 128-45) | (embossed > 128 + 45), e.g.
Or together all of the three thresholds.

That gives the following result:

That approach has two magic numbers, yes. But it works fairly well for such a simple algorithm.
I'm aware the blobs are a bit too big in places. Also, the blue and red channels of the original image don't work quite as well due to more noise in them. Time to experiment further.

Here is a first stab.  It is essentially a doubly smoothed contour map.


Detection Bandwidth for Noise Power Calculation I am trying to calculate the $P_{noise}$ in an experimental setup, such as given
below. The DUT is excited at a predefined frequency, $f$ , with
power $P_\text{in}$. The total power is measured by the RF digitizer. We
can assume that the generator and the sample are ideal so that only noise in the
system , $N_{a1}$, is originating from the Low Noise Amplifier (LNA).

Now, I know that $P_{Na1} = kT_e \Delta f$ and $T_e$ can be found in LNA
specs. What confuses me is the choice of bandwidth.

If we look at the RF front end of the detector;

We see that RF input is first shifted to IF band, then adjusted (to cover ADC full range I suppose) and filtered before ADC.
ADC has $f_{sample} = 250\,\text{MHz}$. Also, the device allows $f_{sample}$ down to 10 kHz.

In view of these informations, here are my questions:

What determines the $\Delta f$ for $P_{noise}$ calculation ? Is it $f_{c}$ of the BPF at the detector front-end ?
What is the role of the sampling frequency? Would it be correct to say $\Delta f =  f_s$ or $f_s / 2$ ?
What is the role of the measurement time? If I measure during 1 s and average the data, can we say that the detector will be accepting noise from 1 Hz window of noise PSD ?
(additionally,) Given the fact that ADC has $f_s = 250\,\text{MHz}$ yet it's possible to sample at $\text{kHz}$ regime, can we say that data is being averaged in case of downsampling ?

The noise of interest would ultimately be over the bandwidth of your signal and specifically all filtering done prior to ultimately making a decision over which symbol was transmitted- which would typically be the bandwidth of your matched filter. So given a noise density in $dBm/Hz$ you can translate that to the total noise in the signal to noise ratio of interest. This is only significant at the point where you make your decision unless you are dealing with waveforms that have negative SNR prior to correlation such as GPS--- in this case you would want to watch the total noise power at the input to each amplifier an ADC to ensure you are not causing saturation and clipping (as the result of too much gain).

For example, the thermal noise floor before we do anything with additional electronics at room temperature is -174 dBm/Hz (given by kT). If we had a 1 MHz signal bandwidth, then this noise level without added contribution would be -174 dBm/Hz + 10Log10(1MHz) = =114 dBm (given by kTB). The added noise by an LNA is the "Noise Figure" and specifically it is how much this noise floor is elevated if the LNA was terminated with it's matched load (meaning no gain prior to the LNA). If the noise figure was 3 dB, which means a 2x increase in power, this would mean the LNA is adding just as much noise as thermal noise is. So this means the LNA self-noise, referred to the LNA input is also -174 dBm/Hz If the gain of the LNA was 10 dB, then the noise power at the output of the LNA would be -174 dBm/Hz + 3dB + 10 dB  = -161 dBm/Hz. If the noise figure of the LNA was 0 dB (meaning perfect, no added noise) the noise power at the output of the LNA would be -164 dBm/Hz.

In your comments you discuss a concern with fluctuation and "removing the noise coming from the LNA". This would have nothing to do with the noise figure of the LNA but sounds like an oscillation, instability or intermodulation spur issue in your receiver (or an interference signal being received, which you can rule out by disconnecting your antenna and terminating the antenna input to the receiver). The noise of the LNA that is in band cannot be "removed". The only significance of the noise of the LNA is to set the sensitivity of the receiver. For this the LNA is just one contribution (a dominant one!) but the noise figure of the overall receiver would be computed as a cascaded noise figure calculation for EVERYTHING up to the point where the symbol decision is made- this includes the LNA, the noise figure of every analog component factored by its gain, the aliasing effects of the mixer stages and ADC, the phase noise and jitter of the clocks, and the quantization noise and further possible aliasing effects in the digital part of the system.

I would suggest you do the following: Determine with measurement the overall noise figure of your receiver to see if you have an issue with this. First measure the gain of your overall receiver by connecting a calibrated (known power level) test tone to the receiver input. Sweep this tone over the full signal bandwidth to ensure there isn't significant variability across the pass band.  Measure this gain at the latest possible point in your digital receiver so that you can capture all effects up to that point in the system. Ensure the input test signal is low enough in power such that the receiver is linear by increasing the input signal 5 dB and ensuring the output increases 5 dB (I would increase to the point where this isn't true and then back-off 10 dB as you want as strong of a signal as possible).  This is the gain, $G$, in dB and can be in any units of dB you choose (dBm, dBFS, dBc etc).  Then terminate the input with a matched load and measure the noise power density in dB/Hz using the same units of dB. This is $N$ in the equation below.

Your noise figure in dB can then be estimated using:

$NF = N- G - 174 $

This would only be as accurate as the gain variation across your bandwidth, which is why this is actually typically done with "white" noise sources that spread the test signal evenly across the bandwidth of interest. But if you confirmed a minimum passband ripple as I suggested, this is valid approach. This noise figure HAS to be higher than the LNA noise figure given the receiver is a cascade of all noise figure computations. A well designed receiver with an interest in maximum sensitivity would keep all additional contributions below an extra 1 dB, but I would not be surprised if a low cost receiver had an extra 2 to 3 dB of degradation.


revert back to stft from abs and phase I have the following code on python
stftspectro = librosa.stft(audio, n_fft=fft_length, hop_length=stride, win_length=window, window='hann')
abs = np.abs(stftspectro)
angle = np.angle(stftspectro)

Can i revert back to stftspectro from below abs and angle. I have no knowledge about sound signals but spectogram from stft is great for audio properties as i have read, and i am applying some deep neural networks that rely on such features.
Moreover, doing the inverse stft will produce the same sound:
stfttoaudio = librosa.istft(stftspectro, hop_length=stride, win_length=window, window='hann')

so collecting both abs and angle is helpful in my thought.

As Marcus Müller said, there is an error, when converting abs and angle values back to complex numbers. It should be like this:
\begin{align}
x = |x|(\cos \phi + j \sin\phi)
\end{align}
Note, that the result convstft of this conversion will be close, but not equal to the original values in stftspectro, most probably, due to inaccuracy of floating point computations.
Here is a reproducible example:
import numpy as np
import librosa

# Creating some audio signal
y, sr = librosa.load(librosa.ex('choice'), duration=10)
tempo, beats = librosa.beat.beat_track(y=y, sr=sr)
audio = librosa.clicks(frames=beats, sr=sr)

# Applying short-time Fourier transform with default parameters
fft_length = 2048
stride = 2048 // 4
window = 2048

stftspectro = librosa.stft(audio, n_fft=fft_length, hop_length=stride, win_length=window, window='hann')
abs_val = np.abs(stftspectro)
angle = np.angle(stftspectro)

# Converting abs and angle back to complex number
convstft = abs_val*(np.cos(angle)+1j*np.sin(angle))

# Checking results
# Choosing an arbitrary element index
n=50
m=45

# Is this element similar?
print("stftspectro element:",stftspectro[n][m])
print("abs_val element:",abs_val[n][m])
print("angle element:",angle[n][m])
print("convstft element:",convstft[n][m])

# Are stftspectro and convstft the same?
A=stftspectro
B=convstft
eps = 0.0001+1j*0.0001
print("Are they similar: ",np.array_equal(A,B))  # test if same shape, same elements values
print("Are they close: ",np.allclose(A,B,atol=eps)) # test if same shape, elements have close enough values

Also, I wouldn't use abs as a name for a variable, since it's the name of the python method abs().


Downsampling oversampled image (smallest size which retains image 'detail' within acceptable limits) I have a large collection of oversampled images (family photos of varying 'sharpness' scanned at a uniformly high dpi) and would like to downsample them to optimise storage space whilst retaining detail within some reasonable threshold.

I am hoping to find (or build) a tool that can help me identify what I (a confessed layman) would call the 'natural resolution' of the scanned image.

For example a 70 year old black & white 6x4" print scanned at 600dpi might be downsampled to 253dpi and look 'sharp' (for some definition of sharp). The downsampled image might then be upsampled to 600dpi again with a deviation from the original that falls within an acceptable threshold. Downsampling the original to 252dpi then upsampling might result in too much deviation, therefore the tool may decide that 253dpi is the sweet spot for that image. Ideally the tool could determine the sweet spot without having to iteratively downsample and upsample.

There are many complications (noise, dust, sharp edges that represent the edges of the photo instead of the photo itself) so the tool may require user supervision. At the moment I'm just looking for pointers as to what techniques and terms I should be familiarising myself with. I have put many hours into googling this issue but all results are either overly simple pc mag stuff (e.g. 'always scan at X dpi') or impenetrable research papers.

Edit (clarification, 2014-01-11)

The paragraph above about downsampling then upsampling is just to illustrate my point about determining loss of detail. In my actual workflow I only intend to downsample (no subsequent upsample step). The goal is to find an algorthm that on a per-image basis will determine the lowest resolution I can downsample to without sacrificing image detail (or sacrificing within an acceptable limit).

There are many options yet it is highly dependent on the properties of the image being resampled.  Images downsampled with bicubic interpolation are regarded as smoother and with fewer interpolation artifacts. You may also try nearest neighbor interpolation and bilinear interpolation. Make sure a low pass filter is applied before further processing with downsampling strategies, otherwise the aliasing artifacts may be induced.

Besides, you may also consider vectorization algorithm  in which a resolution independent vector representation is created before rendered to the image at a desired resolution. it is mainly used in up-sampling but the vectorization method is helpful in preserving the feature connectivity in down-sampling as well.

Another interesting blog on comparison between different image re-size approaches.

EDIT

Regarding the design of low pass filter served as a pre-processing method before downsampling, you can consider this way:

Suppose a low pass filter $g$ works on the image $x$, and the filtered image is $u = x*g$. After that, simply decimate $u$ by 2, you get $v$. Now you want to recover the original image. You can zero-upsample the $v$ into $w$ where $w_{2n} = v_{n}$, and $w_{2n-1} = 0$. Then another filter $h$ is applied, trying to obtain $\hat{x} = w * h$. After Z-transform, $U(z) = X(z)G(z)$, $\hat{X}(z) = W(z)H(z)$.

Since $W(z) = 1/2(U(z) + U(-z))$, we have $\hat{X}(z)=1/2(H(z)G(z))X(z) + 1/2(H(z)G(-z))X(-z)$

A good downsampling is equivalent to a good recovery to some extent. In order to best recover the original image ($\hat{X}(z)=X(z)$), you are trying to make $H(z)G(z)$ as close to $2$ as possible, while $H(z)G(-z)$ as close to $0$ as possible. The design of low pass filter $G(z)$ is based on this criteria.


Extending Background Image Suppose I have the following image. I want to extend/expand the grey background behind the model. The extension should only be upward and sideward. I don't want to extend downward because it will be ugly.

The purpose is to make the image into size (w, h) where w is greater than or equal to current image width and h is greater than or equal to current image height.

One way I could think is, to copy the pixel in the (0, 0) and fill the image extension using this pixel value. However, as you might notice, the grey background at bottom right is slightly different than the rest. This will make noticeable colour difference.

What are better algorithms that is usually used to do that? The resulting image should not have noticeable colour difference and looks natural.

What you're looking for basically is "Inpainting".

In Photoshop is called "Content Aware Fill".

You should look for articles for Inpainting.

Trivial option would be using "Replication" of the border and later enforce some smoothness (Either by Gaussian Blur or even Poisson Image Editing).

P.S.
We'd be happy if you could join our Community Proposal at http://area51.stackexchange.com/proposals/86832 for Image Processing and Computer Vision.
We need more people to follow, up vote suggested questions with less than 10 up votes and suggest new questions.
Thank You.


Scale and Rotation invariant feature descriptors Can you list some scale and rotational invariant feature descriptors for use in feature detection.

The application is for the detection of cars and humans in video captured by a UAV, using a multi-class classifier.

So far I have been looking at SIFT and MSER (which is affine invariant). I have also looked at LESH, LESH is based on the local energy model, but is calculated in a way that is not rotationally invariant, I have been trying to think of a way to make use of the Local energy, to build a rotationally invariant feature descriptor, I read here What are some free alternatives to SIFT/ SURF that can be used in commercial applications? ,that " if you assign orientation to the interest point and rotate the image patch accordingly, you get rotational invariance for free", but don't know if this is even relievent or how i could apply this to my problem, any help would be appreciated, thanks

As far as alternatives to SIFT/SURF go, the question you linked provides very good answers.

There were two more questions I could read out:

"how could I build a useful (e.g. rotation invariant) feature descriptor"?
"regarding the statement from the linked question, how does he accomplish free rotational invariance?"

Building feature descriptors

This is a valid research topic. Good feature descriptors are not something just anyone can build in an afternoon. People publish articles when they successfully model feature descriptors with desirable properties. This is a reason that currently only a handful of state-of-the-art descriptors are used, and that is also what I advise you to do: find a feature descriptors that's good for your needs.

Achieving rotational invariance "for free"

You can determine the dominant gradient, or orientation at an image patch (your feature area). Then, rotate the image patch so that the gradient is always looking in the same direction, e.g $0$ (upwards). E.g. if you had a |black->gray->white| and a |white->gray->black| image, their dominant gradients would be pointing left ($-90$) and right ($90$) respectively, and when you rotate them by this amount, you get the same images.

This way you will always calculate the descriptor on an image patch with the same dominant orientation (the rotated patch), and thus you achieve rotational invariance.


Isolate recurrent pattern vs. evolving background sound Here is a sound example for what follows.

Let's assume we have a signal

$$s(t) = r(t) + e(t)$$

where:

$r(t)$ is a signal which is recurrent with a given period, i.e. in my example $r(t) = $ the exact same snare sample, which is played every 1 second, starting at 00:00:00,500
$e(t)$ is an evolving sound (in my example it's the synth pad soft sound)

I assume the signal is mono, and that no compression or other master effect has been applied.

Question:
How to recover $r(t)$ from the mix $s(t)$, by analyzing which "pattern" is recurrent in $s(t)$? (the period 1 second is considered as given)

In one word, how to isolate the snare sound from the mix?

I would like to do it by slicing the signal each 1 second:

slice1.wav
slice2.wav
slice3.wav
slice4.wav
slice5.wav
slice6.wav
slice7.wav

and analyze what part of signal is recurrent (= the snare) and which is not (= the synth)  among all slices ?

Is there a separation algorithm based on analyzing recurrent parts / non-recurrent parts? Using FFT or something else?

So, assuming the linear superposition you suggest, the sound you want to exctract lives in a periodic subspace of your signal space. This subspace is uniquely characterised by the exact spacing between the sound instances, which I assume you know in the following. If you don't know it, you can find it from calculating the autocorrelation sequence of your composed signal.

Without further knowledge, the only linear method for finding your sample is to project the composed signal onto the periodic subspace it lives in. There are two obvious ways you can do this.

The first is taking the Fourier series of the composed signal and nulling all coefficients that do not belong to an integer multiple of the frequency corresponding to the period length of the repetition. Summing the modified series will give you the desired signal.

The second is to just accurately cut the sample into segments that align perfectly with respect to the repetition period and just sum all those segments up into a single sample.

The result of these two operations will give you a representation of the composed signal the same periodic subspace. However, since the subspace is relatively large compared to the total signal, there will be components of the composed signal that are contained in this subspace by mere chance. So the separation will be far from perfect. You can improve the separation by taking more repetitions, make the full signal space much larger than the periodic subspace.

If you need a better method you won't get around characterising the nature of the signal you want to extract. That's where perception models, source models, statistical signal models, etc come in. I cannot really go into this in this answer. It's a topic that fills many books and is in general very difficult.


O-QPSK constellation interpretation - GNU Radio I am currently trying to recover the I&Q of a O-QPSK signal generate from the standard XBEE arduino module using a B200. I can get somehow get a "look-like" of a constellation. Here an screenshoot:

I highlight the zone where there is a symbol recover. When I run the acquisition I can see symbol at the places. Here we just see one because there is a limited number of points displayed. However what is the rest? I got a SPS (sample per symbol) of 2. A part of the rest is certainly the transistion, but I am quite sure there is phase noise also. How can I attenuate those points. My flowgraph is just a UHD source feed into the polyphase clock source which is sent to a costa loop and after that there is my QT constellation sink.

[EDIT]
A other picture with less points on it:

The phase recovery was finaly perform in Matlab. Because the signal we are studying (Zigbee) send a preamble before the communication start, we can compensate for the group delay and the channel's distortion with a anti-causal analysis of the expected signal and the received signal.

The GNU-Radio's blocks  weren't efficient at recovering O-QPSK's signal. However with fact that we don't real time treatment of the signal, Matlab was finaly most suitable way to resolve the O-QPSK demodulation of the Zigbee signal.


Is this graph of Wigner-Ville distribution OK? I am studying some code for smoothed Wigner-Ville distribution.

update 3- solution, graph was over amplified in a 2d overdrive/oversaturation effect. if i divide by 20 only a line is visible.

Update - indeed there is a bug, the current image is from 128 samples of a $44.1\textrm{ kHz}$ signal where sample[i] = sin[i/10.0]; i.e. a $4410 \textrm{ Hz}$ sine wave, to fit inside the 128 samples.

As suggested by JazzManiac, here is the code which is causing errors, it's actually P Kootsookos Cohens.C code, with a minor change, instead of a .raw input audio TFD file i am using a sine wave array in get_my_signal function... A trivial and simple patch, wouldnt you agree? that can't be the error?!

I have a professional C programmer who has been coding for 30 years who provided this version of the code, it's just setting the argc and argv options in advance, to spwvd, hamming, and default window (was 128, is now 1024).

Is smoothed Wigner-Ville the same as smoothed pseudo Wigner-Ville?...

The original author (i.e. me) gets this for each of the three test cases:

The graph is over amplified. If the image in blue is divided by factor of 20, the transients vanish and only a line of the detected frequency is visible.


Quantifying the energy of a voiced frame What should the energy, pitch frequency, and zero crossing rate of a voice/unvoiced/silence frame in signal be?

My professor said that it is experimental, but I do not have any idea what the factors of these are and where I should begin from.

I would appreciate any guesses based on experimental data or any limits on the quantities.

I've tried a zero crossing rate between 2000 and 94000, energy between 0.0001 and 0.7, and
pitch frequency between 180 and 500.

Are you being asked to write a simple vocoder?

Your question is somewhat vague and confusing, but it seems like you are being asked to determine parameter thresholds for deciding when a speech frame is voiced (has fundamental periodic component and harmonics), unvoiced or frictional (noise generated), and silent (do not encode).

It looks like you would be using the pitch frequency estimates to determine if the frame is voiced/unvoiced, and signal energy to determine if it is a silence frame.

If he said it's experimental, I think he is asking you to look at the energy and pitch content of an audio signal and see if you can see differences that will help you choose thresholds.  For example, get some audio of someone speaking.  When the person is vocalizing (i.e. vowel sounds like "aaaahh" and nasally sounds like "mmmmmmmm") the spectrum should show periodic components.  For unvoiced sounds (i.e. short "th" "t" "s" "f" type sounds) the spectrum will look more like noise.  And when there is silence, the signal should be relatively low energy.  Take what you find and use that information to determine your thresholds.


What does the exponential term in the Fourier transform mean? We know that Fourier transform $F(\omega)$ of function $f(t)$ is summation from $-\infty$ to $+\infty$ product of $f(t)$ and $e^{-j \omega t}$:

$$ F(\omega) = \int\limits_{-\infty}^{+\infty} f(t) \ e^{-j \omega t} \ dt  $$

Here, what does the exponential term mean?

It's a complex exponential that rotates forever on the complex plane unit circle:
$$e^{-j\omega t} = \cos(\omega t) - j \sin(\omega t).$$
You can think of Fourier transform as calculating correlation between $f(t)$ and a complex exponential of each frequency, comparing how similar they are. Complex exponentials like that have the nice quality that they can be time-shifted by multiplying them with a complex number of unit magnitude (a constant complex exponential). If the Fourier transform result at a particular frequency is a non-real complex number, then the complex exponential of that frequency can be multiplied by that complex number to get it shifted in time so that the correlation to $f(t)$ is maximized.


What technique can I use to remove useless information from an image? Is there any way to remove useless information (unimportant details) from an image?

If, for example, we have a specific image, so this contains both useful and useless information. If this image can be recognized directly without the useless information, the unnecessary information can be removed (e.g., illumination, some noises, variations, etc.).

Although you did not spent even a minute in researching the question I will post an answer to it. There are multiple ways; I will try to demonstrate it using wavelets in Mathematica.
So, first of all we need an image.

img = ExampleData[{"TestImage", "Mandrill"}]

Then we apply the DiscreteWaveletTransform using the HaarWavelet

dwd = DiscreteWaveletTransform[img]

We plot the transformation

WaveletImagePlot[dwd]

We threshold it, i.e. removing the noise according to some criteria

wtd = WaveletThreshold[dwd]

Again plot the transformation and the filtering

WaveletImagePlot[wtd]

We now inverse the transformation

iwd = InverseWaveletTransform[wtd]

And visualize the differences between the two images

ImageDifference[iwd, img]


Extracting ridges in automatically in image I am trying to extract ridges from an image, calculated using the eigenvalues of the Hessian matrix (see answer in Best way of segmenting veins in leaves?). In this image I have zeroed out the negative eigenvalues.

Visually I can see the regions I want, but how would I be able to sensibly extract them automatically? If I LogNorm() the colorbar it is not so clear

Initially I thought if the noise in the image can be estimated, I can mask the image accordingly, but the only estimate to the noise I can find is its standard deviation (assuming Gaussian noise) in skimage.restoration.estimate_sigma and not its mean.
Masking via threshold is also possible, but do the magnitudes of the eigenvalues have physical meaning?
I'm new to image processing so any help would be greatly appreciated. I have attached the data for my image also (https://drive.google.com/file/d/1Wdvqtv4rUmLtSNiq-p0M5BwSKxUoW76T/view?usp=sharing).
Thanks!

This should be doable with ssqueezepy's extract_ridges; try varying penalty and bw (see their docstrings). As last resort, feeding a cropped image that excludes region without ridges may work better, as the algorithm assumes the ridge spans the entire frame. You can automate this by finding indices at which column energies fall below a set threshold, e.g. np.where(np.sum(np.abs(stft_transform)**2, axis=0) < 0.1).
'Cropped' example:


Evaluating different inputs into a given system We need to find which signal wont be the same after going through this sytem.

My approach is this one:

I took the filters frequency and did 2pi * 1/4 and gave me this "limit" 1.5707

(2pi because the periodicity and 1/4 because the filter "bandwith"?)

Then i proceed with the inputs evaluating them this way:

a) 4 * pi/7 = 1.7951
b) 4 * pi/9 = 1.3962
c) 4 * pi/8 = 1.5707
d) idk

The result "a", which is out of the "limit" is the one that does not follow the rule y[n]=x[n] and thus, the correct answer. (It is the actual correct answer)

Did i have luck there with the math or the approach is correct?

The output will be equal to the input $(y[n]=x[n])$ if no aliasing occurs as a result of downsampling. Which means $X(e^{j\omega})=0$ for $\frac{\pi}{4} \leq |\omega| \leq \pi$.
a) $X(e^{j\omega})$ has impulses at $\omega = \pm \frac{2\pi}{7}$ $=>$ aliasing exists $=>$ $y[n]\neq x[n]$
b) $X(e^{j\omega})$ has impulses at $\omega = \pm \frac{2\pi}{9}$ $=>$ no aliasing $=>$ $y[n]= x[n]$
c) $X(e^{j\omega})$ has impulses at $\omega = \pm \frac{\pi}{8}$ $=>$ no aliasing $=>$ $y[n]= x[n]$
d) You need to sketch $X(e^{j\omega})$ to see if aliasing exists  hint: the impulse response of your filter is $h[n]=\frac{sin(\pi n/4)}{\pi n}$


Calculating similarity between bit strings/signal I have a project I'm working on where I'm trying to efficiently compare sequences of binary data to determine how similar they are. Basically, a sequence might look like this:
1001010010100011110...

I spent a while on my own trying to come up with a way to do it because it's fun for me to work through problems like this on my own, and actually made some progress. Unfortunately, everything I came up with never quite worked; a solution that accounts for substitutions (bits that flip) is easy, but what keeps getting me is handling bits that are inserted or deleted, i.e.:
100101 -> 1000101

Trying to find a solution online keeps leading me here, and I keep coming across the term cross-correlation. I also see lots of interesting links on the wikipedia page for Time Series for calculating similarities. Unfortunately, at the moment, a lot of the math is unclear to me, and I can't efficiently go through and distinguish what the appropriate technique might be for my situation.

So I have two different approaches to my problem that I'm interested in learning about:

1) What method would you recommend for calculating the similarity between two bit sequences (typically 50 to 150 bits in length) that is robust to substitutions, insertions, and deletions?

2) I'm not sure this is possible, but is there some kind of distance metric(s) I can calculate independently for each signal that will allow me to group them together with actually doing direct pairwise comparisons?

For 1, I feel like cross-correlation might be the right approach, but I haven't gotten past just plugging sequences in R and seeing how it works for identifying lag. I know edit distance will work, but ultimately will want something more efficient. Eventually, I will be trying to group billions of these together.

Thanks

You should indeed look into sequence alignment methods used in bioinformatics. Those can penalize both for insertions/deletions and mismatches. The binary alphabet having a size of 2 and that of DNA having a size of 4 does not require the use of different algorithms if the element-wise comparison results in either "match" or "mismatch", and a mismatch is penalized for. Check out the Needleman–Wunsch algorithm, which is based on dynamic programming.

Some DNA/protein sequence alignment algorithms such as MAFFT use Fast Fourier Transform (FFT) for accelerated cross-correlation of residue properties to approximate match calculations.


Why can't I add variable covariance noise to this 2D kalman filter? The goal is to include variable measurement noise to my kalman filter.
To begin with I kept the noise covariance matrix as follows:
$$
R=\begin{bmatrix}1&0&0&0\\
0&1&0&0\\
0&0&1&0\\
0&0&0&1\end{bmatrix}$$

This gave good results, the Kalman filter did a good job following the input data.

Next I tried to calculate the measurement noise covariance myself:
https://pastebin.com/7ETMC0Fy

Which gave such results:

Could someone tell me why when I try to calculate my measurement noise covariance matrix myself the result is unsatisfactory? Why doesn't the output of my Kalman filter follow the input data at all?

Instead of:

x = detrend([1, 5, 10, 15, 20,  1, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80]);
y = detrend([1, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80]);

do:

x = [1, 5, 10, 15, 20,  1, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80];
y = [1, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80];
dx = detrend(x);
dy = detrend(y);

and then:

  myMeanX = mean(dx);
  myMeanY = mean(dy);
  varianceX = 0;
  varianceY = 0;
  for i=1:length(dx)
       varianceX = varianceX + (myMeanX - dx(i))^2;
       varianceY = varianceY + (myMeanY - dy(i))^2;

  end

and

R(1,1) = varianceX/length(dx);
R(2,2) = varianceY/length(dy);

and I'm not sure why you're doing this:

%set everything to zero except diagonal

P should be as-calculated. No need to zero out entries.

The code you show must be wrong because:

[31.5202     0     0    0
   0       0.0461  0    0
   0       0.0458  0    0
   0        0      0    0.0051]

shows you are setting off-diagonal entries in R.

Compare:

varianceX = varianceX + (myMeanX - x(i))^2;

with

varianceVx = varianceVx + (myMeanVx - VxMeas(i));

Notice anything²?

Also, it would be good if you could actually print out the values of your estimated R matrix.

² $\tiny \mbox{If you don't, please read it again.}$


Receiver for low data rate I am working on satellite receiver burst mode
I am trying to build receiver. My data rate is typically 10kpbs.
As the data rate is low and carrier frequency offset may also be around 10kbps

so i have implemented reciever please let me know the blocks mentioned below are ok

coarse frequency recovery
carrier recovery
timing recovery
matched filter

i have seen many designs where timing is done before carrier recovery.

what is the advantage of doing carrier before timing? is it because offset is more then do carrier recovery first otherwise vice versa.
Will the nested loops help
but quantitatively how much symbol rate to freq offset can be recovered like this.
Once recovered the burst length can be as long as possible.

This is very implementation specific and there is not one rule that one should be done before the other, and joint estimation another option where there isn't a clear order of events.

A consideration is how well your carrier tracking approach works under timing offsets, and similarly how well your code tracking approach works under carrier offsets. I recommend analyzing/simulating this specifically for your implementation since it is very application specific and you didn't provide any info about your waveform or algorithms for timing and carrier recovery. For example, the Gardner Timing Error Detector works very well with relatively large carrier offsets but has higher noise than the Mueller and Mueller detector (M&M), but the M&M being decision based does not perform well with relatively much smaller carrier offsets.

One suggested approach to quantifying this for your specific implementation/application: determine your pre-decision signal SNR with fixed SNR received signal input as you vary both carrier offset and timing offset. Also you can look at the noise of your symbol timing and carrier error estimators as another metric: determine what the noise is on the symbol timing error as your carrier offset is varied (When you are close to your track point for tracking considerations, and when you are further away for acquisition considerations). Similarly, determine what the noise is on the carrier error as your timing offset is varied. If you do this work, please do share it with us here.


Calculating power and energy of a signal in matlab I am a little new to matlab and would like to get some help on how to write a program or a function that can calculate the energy/power of a signal depending on whether the signal is an energy or a power signal.
Here is the program the I tried:

function [ ] = power_or_energy( s,t)
%UNTITLED Summary of this function goes here
%   s is the input signal
%   t is the independent variable

s2=s.^2;
energy=int(s2,t,-inf,inf);

if energy~=inf && energy~=-inf
    disp ('the signal is an energy signal with energy equal to ' );
    disp(energy);
end

power=limit((int(s2,t,-t/2,t/2))/t,t,inf);

if power~=inf && power~=-inf
    disp('the signal is a power signal with power equal to ');
    disp(power);
end

But the problem here is that when i want to check it for a function like

t=0:0.01:10;
s=sin(t);
power_or_energy(s,t);

then i get the following error:
"Undefined function 'int' for input arguments of
type 'double'."

So how to i actually calculate the energy for such signals?

Also the above function works fine with the below input:

syms t;
s=sin(t);
power_or_energy(s,t);

Please help me in writing a function that can be used to find the energy of any input signal.

The Matlab function int needs a symbolic expression as an integrand, and a symbolic variable as an integration variable. In your example s and t are both vectors, and you can't use int in this case. If you want to compute the power or energy of discrete time signals, then you need to use the corresponding definitions:

$$E_x=\sum_{n=-\infty}^{\infty}|x_n|^2\\
P_x=\lim_{N\rightarrow\infty}\frac{1}{2N+1}\sum_{n=-N}^N|x_n|^2$$

In practice you will need to consider finite length signals unless you have an analytic formula for $x_n$ and some software package that can deal with the analytic evaluation of infinite sums. Obviously, finite length signals always have finite energy and $P_x=0$ according to above definition.