Filtering noise and identifying phones boundary in speech recongition I am planning to do a self project on speech recognition for a vocabulary detection system using arduino. I know CMUsphinx or bitvoicer is probably the way to go. But I would like to code the project from the ground up as the purpose of the project is learning.

I have been reading different papers like

http://revistaie.ase.ro/content/46/s%20-%20furtuna.pdf (dynamic time warping)

http://citeseer.ist.psu.edu/viewdoc/download;jsessionid=DEC7F40E5AE89F1508A3B6D0625CC1DA?doi=10.1.1.131.2084&rep=rep1&type=pdf
(HMM)

But I have encountered the following difficulties

I have a hard time reconciling the entire process with higher level algorithm. From what I have read so far the entire process is the following

i. input signal

ii. filter all noise, getting all signals that is high energy volume and low frequency. (How do I know what to compare against? What is consider high energy volume, and low freqency?)

iii. divide the signal into phones using FFT (where do I get dataset for phones, what constitutes a boundary between two phones, and how similar does the detected phones have to match the dataset?)

iv. matching the list of phones to each word using HMM or DTW

For beginners like me I feel that steps 2-3 is the hardest part due to the difficulties I listed.

Can someone point me a direction?

Thanks

system using arduino

Arduino performance is not enough for any serious speech recognition. You should better play with your mobile phone instead. With mobile you can build at least some interesting applications. Or try Raspberry Pi.

  From what I have read so far the entire process is the following

You probably need to read more papers (including the source code of existing software like CMUSphinx) to understand the details because the description of the process you have is not very accurate either. Ideally you need to read a textbook like "Fundamentals of speech recognition" by L. Rabiner.

  i. input signal

This is two lines of code.

  ii. filter all noise, getting all signals that is high energy volume and low frequency. (How do I know what to compare against? What is consider high energy volume, and low freqency?)

This is not very accurate description, there is no such thing as "low frequency". There is a voice activity detection (VAD) task which is pretty complex and has many dedicated algorithms. The simple form is VAD based on energy of the signal. You compare background energy level with current energy level and if the difference is higher than threshold you count audio as speech. This is not very accurate VAD but it is ok for many cases (no big noise, clean audio). It can be extended to frequency subbands and mixed with noise suppression by spectral subtraction. Noise suppression combined with VAD is even more complex task to solve.

For simple start you can implement energy-based VAD (google can give you many examples of implementation)

  iii. divide the signal into phones using FFT (where do I get dataset for phones, what constitutes a boundary between two phones, and how similar does the detected phones have to match the dataset?)

For DTW comparison with the prerecorded audio you do not need to care about phones at all. For HMM-based recognition, HMM decoding algorithm finds boundaries between phones implicitly but actually it doesn't even need them. Speech is considered as a sequence of subphone states, not as phones. Phones are not very stable to model with HMM. You can find description of the algorithm in the Rabiner's paper linked above.

The database of phone properties is extracted from transcribed recordings. There are many public databases available (some are mentioned on cmusphinx website).

  iv. matching the list of phones to each word using HMM or DTW

There is no such thing neither in HMM nor in DTW. You match the audio to the word model, not taking individual phones into account. You will have to implement a separate training algorithm.

  For beginners like me I feel that steps 2-3 is the hardest part due to the difficulties I listed.

Since your process is not very accurate the difficulty estimation is not very meaningful too. HMM training and decoding are quite complex tasks, VAD or FFT are simple.

The topic is not very easy, it is complex. You can not expect to implement calculus algorithm without knowing the theory of calculus, like what is integration and so on, so the speech recognition. There are many tutorials around, this one is ok to get started too:

http://kastnerkyle.github.io/blog/2014/05/22/single-speaker-speech-recognition/

Maybe it will be helpful for you. However, take into account that any simple tutorial is usually not very professional and doesn't really cover things. Speech recognition can not really be explained in 10 pages of text. If you really want to understand a topic, read a textbook or attend a course on ASR.


Phase spectrum range I computed spectrum analysis of sine wave (without windowing and zero-padding) and I got plots as follows:

As far as I know they are correct. When I learnt about DSP, I got the huge problem with understanding what is the phase spectrum. Recently I realised that "phases" are translations along time axis of each sinusoidal component. So I don't see any point to get them out of range $\langle -\pi, \pi \rangle$. I computed my spectrums with python and to get the phase spectrum I used unwrap function from numpy. As I understand this documentation correctly, it is impossible to get the result larger than $\pi$ or less than $-\pi$ unless I set the discount parameter. It does make sense for me.

Could you explain me:

whether I understand the phase spectrum correctly?
why the values of phase spectrum are out of range $\langle -\pi, \pi \rangle$?
how does the function unwrap work?

Yes you can think of phases as describing translations.

If you (kind of passively) wrap the phase to ⟨−π,π⟩ you may be introducing discontinuities that are not truly there considering that −π and π are equivalent. Such artificial discontinuities will hinder calculation of things like group delay.

Unwrapping the phase properly is not trivial. One way to do it properly is to factor the system into such small component systems (connected in series) that you can track the phase frequency response of each without any wrapping ambiguities. The phase frequency response of the composite system will be the sum of those of the component systems. Calculated this way, the phase naturally get values outside of ⟨−π,π⟩. Computationally, the factorization can be done by finding the z-plane zeros of the system. My guess is that numpy.unwrap does nothing of the sort but simply looks for big jumps of the wrapped phase as function of frequency and tries to make those jumps shorter by adding or subtracting 2π. That will easily fail do do a proper job if the frequency response has two zeros at the same frequency, as that gives a jump of 2π that phase unwrapping ought to preserve.


How can I get maximum complex value from an array of complex values(524289) in python? I have an array of complex values. Actually this array contains fast Fourier transformation fft values. now I want to fetch maximum complex value from the array. Thanks

"Maximum" is ambiguous when it comes to complex values.  The complex value itself doesn't have a min or max.  Which is greater, 1+0j or 0+1j?  What about 1+0j and 0+.5j?  The answer to these questions determines what exactly what you want to do.

You can get the maximum real part (max(1+0j, 0+10j) == 1), maximum imaginary part (max(10+0j, 0+1j) == 1), maximum absolute value (max(10+0j, 0+1j) == 10, max(1+0j, 0+10j) == 10), complex value with the maximum real part (max(1+0j, 0+10j) == 1+0j), complex value with the maximum imaginary part (max(10+0j, 0+1j) == 0+1j), or complex value with the maximum absolute value (max(10+0j, 0+1j) == 10+0j, max(1+0j, 0+10j) == 0+10j).  All are possible with numpy arrays.

The default in numpy (arr.max() or np.max(arr) when arr is complex) is the complex value with the maximum real part.

import numpy as np

arr = np.random.random(1000)+np.random.random(1000)*1j  # generate example data

maxreal = arr.real.max()  # maximum real part
maximag = arr.imag.max()  # maximum imaginary part
maxabs = np.abs(arr).max()  # maximum absolute value

maxcompreal = arr[arr.real.argmax()]  # complex value with maximum real part
maxcomp = arr.max()  # complex value with maximum real part, same as above
maxcompimag = arr[arr.imag.argmax()]  # complex value with maximum imaginary part
maxcompabs = arr[np.abs(arr).argmax()]  # complex value with maximum absolute value


Block normalization of histogram of oriented gradient I am a beginner in image processing. I have some questions regarding HOG.

https://www.learnopencv.com/histogram-of-oriented-gradients/

From this website, it mentioned that each element inside the HOG descriptor will be divided by the sum of all elements in a block. However, how to do block normalization if overlapping the blocks? I do not understand because some elements in the vector will be elements of more than one blocks.

Picture of block normalization

For example, if we have a histogram [v11 v12 v13 v14...v19] ,[v21 v22 v23 v24...v29], [v31 v32 v33 v34 v35...v39] and [v41 v42 v43 v44..v49] in the red block, we can normalize the first element of the first vector like this: v11/(v11+v12+...+v21+v22+v23+..+v31+v32+...+v41+v42...+v49)

However, how can I normalize the vectors in the red block and yellow block?

Thank you very much.

I stucked at this question too, because overlapped blocks will be divided twice, right?

and so I found these :

1 http://www.geocities.ws/talh_davidc/#cst_extract

2 https://stackoverflow.com/questions/32417531/hog-what-is-done-in-the-contrast-normalization-step

to verify my answer, I done this:
on 1, they using image size HOG126x63(), block size is 18x18, with 50% overlapping blocks means it had 9px stride.
126/9, 63/9 = (13, 6) blocks.

they mentioned that final vector would be 6318 = 78*81 == 78*9*9 (78*number_of_cells_in_one_block*number_of_bins)
turns out, 13*6 = 78, so final vector is
(blocks_in_image_col *blocks_in_image_row *number_of_cells_in_one_block *number_of_bins)

conclusion is, one block is counted as one set vectors (2*2*9 as in image) , moving the blocks will produce one set of vectors again. there won't be cells normalized twice (it will make the histogram spike at first and flat on remaining due to recursive dividings)


Denoising approach for a combination of several ADC voltage channles I have 2 ADC channels of constant voltage measurement with small amount of high frequency noise only and no low frequency oscillations. The final signal should be the simple sum of those and denoised with the moving average in the end. Does the moving average denoising of each separate channel before final summation be more noise effective than applying summation first and do the final moving average later?

channel_1_noisy = channel_1_clean + noise_1
channel_2_noisy = channel_2_clean + noise_2

sum_noisy = channel_1_noisy + channel_2_noisy
sum_denoised = moving_average(sum_noisy)

or sacrificing computer time for 2 moving average filtering operations?

channel_1_denoised = moving_average(channel_1_noisy)
channel_2_denoised = moving_average(channel_2_noisy)
sum_denoised = channel_1_denoised + channel_2_denoised

Having more channels will the effect be exaggerated?

According to your processing chains (filter --> sum vs sum --> filter) the two approaches, by being LTI systems, should produce the same results under infinite precision arithmetic. You can show this mathematically by using the properties of convolution operation.

Note that a moving average filter is just an ordinary LTI system.

Note also that nonlinearities such as quantization operations may affect the above fact's validity.


Compute RMS value of non-stationary signal I would like to know what is the most accurate way to compute a measured signal's fundamental RMS value as a function of time.

The following is true for the considered signal:

signal represents a voltage
sample rate: fs = 6 kHz
capture time: t = 100 s
fundamental frequency: f1 = 60 Hz
signal is not a pure sinusoidal but has additional frequency content
amplitude of the signals fundamental changes slowly over time

The general way to obtain high resolution envelope of a signal is called envelogram which is the magnitude of the analytic signal obtained from your original signal(the analytic signal obtained by using Hilbert transform).

Simply take FFT of whole signal to obtain the spectrum, then replace the value of spectrum over negative frequencies with zero (keep the value of spectrum over positive frequencies as before). Now taking inverse FFT will give you a complex signal which its magnitude is the envelope of your signal which is the instantaneous power or RMS.
By applying a band pass filter around your fundamental frequency you could obtain the RMS of your fundamental frequency.

Another way is to use time-frequency analysis which is similiar to the method you suggested (it uses a moving window instead of splitting signal)


Is it possible to infer the characteristic curve of a recording device only from (many) recordings? ... with a certain degree of confidence / corridor of accuracy.

What do I mean?
Case in point: There is the huge archive of digitized Edison wax cylinder recordings.

They have annotations, usually year of recording, sometimes producer of recording.
Assuming that, with a good likelihood, grouping the intersection of e.g. one year of recording and a geographic area or even company yields sets of recordings very likely made with the same model / iteration of the technology within a set, given how slow-paced development was back then:

Could one infer the rough frequency response curve of an imaginary "average recording device" of those times from such a narrowed-down set of recordings?

My not exactly DSP expert naive train of thought went like this:

if I have a great enough number of as diverse as possible audio recordings, their sources (i.e. before entering the brass horn to be recorded) should have an average spectrum roughly resembling white noise. This probably means using broadly set up orchestral recordings, not so much "solo singer with background accompaniment", which is weighted too much to vocals in the spectrum.
with a good likelihood of being made on the same model of recording apparatus, an "average spectrum" calculated e.g. with FFTs over the whole set of recordings as described above, should yield the rough response curve of the imaginary recording device

Now, if premise 1. is already wrong, then I'm out of luck.
And 2. could be wrong if the spread of characteristics between individual exemplars of devices of a series is too great - although I guess that's a bit of a question with regards to "what do I need".
Anyway, this was just to show I made some effort thinking about this. If it doesn't work the naive way I am imagining:

Is there a way to do that, only from recordings, without being in possession of an actual device of that type (for measurement)?

The amount of distortion by e.g. saturation within the mechanical recording process can probably not be inferred just from recordings?

Yes, that's definitely possible. The key would be to find recordings where it's reasonably easy to find an acoustic reference. Great candidates for this would be classical or big band recordings. Let's say you have wax cylinder of Beethoven's 5th, you can try to create a reference by analyzing a dozen or so recent recordings of the same piece. The standard deviation would give you a "margin of error" for the reference and than you can simply divide the spectrum of the wax cylinder by the reference spectrum. Repeat for as many recordings as you can find good references for.

  should have an everage spectrum roughly resembling white noise

No. It's certainly not white: closer to pink but with more high end roll off. Hence you need a known good reference to determine exactly what it is.


Can I compute Fourier series without looping through all frequencies? I need to compute Fourier series of an audio stream.
But DFT/FFT is slow.

Are there any ways to compute Fourier series of a signal without using the Fourier transform to check if whether a frequency is present or not?

Make sure you do not mix up FFT and DFT, they gives identical results, but are different.
DFT complexity is about $O(N^2)$ while FFT complexity is $\frac{34}{9}N\log_2(N)$. This makes a huge difference is computation time.

I do not know any faster algorithm for computing a complete spectrum than the FFT. If you need a subset of the spectrum (only few frequency bins), you could look at Goertzel Algorithm.

According to wikipedia, as a rule of thumb, Gortzel algorithm will be faster than FFT if :
$$ M < \frac{5N_2}{6N}\log_2(N_2)$$

Where

$M$ : Number of frequency bin you need
$N$ : Number of point in your dataset
$N_2$ : The first power of 2 bigger than N.


Frequency/harmonic analysis of accents I am looking to find out how accents could/would be detected by a programme, I haven’t as of yet undergone an experiment into detecting or observing anything noticeable about recorded accents. I am not sure a frequency/harmonic analyses will reveal anything. Trying to find any literature review on the subject is sadly lacking, so I was hoping to find out from the community where one might start or what would you look for in any given accent. Something quantifiable could or should surely be observed dependant on an accent?

For instance say I have 10 people from town A and 10 people from town B, all males with roughly the same age Town A and Town B both speak English but have different accents. If I ask each of them to repeat the same word, we should notice a difference in the frequency/harmonic analysis of Town B folk from Town A folk.

What is your best guess at where this could be found? Is frequency or harmonic analysis going to reveal anything between the two groups?

For instance if we take specific vowels we could possibly see a certain threshold reached in amplitude level and frequency (pitch) and could then work out the resonance frequencies of the vocal tract for certain accents based on key vowels? In speech science and phonetics the use of formant is used to mean an acoustic resonance of the human vocal tract. It is often measured as an amplitude peak in the frequency spectrum of the sound, using a spectrogram or a spectrum analyzer, though in vowels spoken with a high fundamental frequency, as in a female or child voice, the frequency of the resonance may lie between the widely-spread harmonics and hence no peak is visible.

If you notice that using a formant table you can construct vowels with a vocal formant filter, could this be reverse engineered to essentially "detect"?

My worst fear however is that any given vowel is going to have either a lack of accent specific information, or the analysis done will have too much sporadic information between individuals.

Could anyone point me in the right direction here im quite stuck?

I would guess that the most prevalent difference between accents (of the same language) is how the vowels are pronounced. Humans distinguish vowels mainly by the frequency of the formant, i.e. regions of high energy in the harmonic spectrum.

A good way to display this is the "formant map" which plots the frequency region of the second formant against the frequency region of the first. See for example: http://www.phy.davidson.edu/fachome/dmb/digitalspeech/formantmap.gif

In this 2D map every map has a specific region. For different accent the centroids of this regions probably shift around. That is certainly true for different languages. An interesting approach would be to measure this map for, say, an Australian vs. an American English speaker and look for systematic differences.


Are Z-transform time shifting and differentiation properties always compatible? I'm currently studying the Z-transform, and I'm having issues in understanding the time shift and differentiation properties, to be precise: calculating a Z-transform explicitly, and obtaining it by using those properties yields different results.

For example, a common Z-transform is the transform of the discrete unit step function:
$$
\mathcal{Z}\left\{ u\left[n\right]\right\} =\frac{1}{1-z^{-1}}=\frac{z}{z-1}
$$
The time shifting property of the Z-transform states (leaving aside its effects on the ROC for the moment) that $\mathcal{Z}\left\{ x\left[n-k\right]\right\} =z^{-k}\mathcal{Z}\left\{ x\left[n\right]\right\} \quad\forall k\in\mathbb{Z}
 $, thus it implies that:
$$
\mathcal{Z}\left\{ u\left[n-1\right]\right\} =z^{-1}\frac{z}{z-1}=\frac{1}{z-1}
 $$
Since $\frac{d}{dz}\frac{1}{z-1}=\frac{1}{\left(z-1\right)^{2}}$, the differentiation property, which states that $\mathcal{Z}\left\{ n\, x\left[n\right]\right\} =-z\frac{d\mathcal{Z}\left\{ x\left[n\right]\right\} }{dz}$ gives:
$$
\mathcal{Z}\left\{ n\, u\left[n-1\right]\right\} =-z\frac{1}{\left(z-1\right)^{2}}=-\frac{z^{-1}}{\left(1-z^{-1}\right)^{2}}
$$
which appears to be wrong, since the correct result should be:
$$
\mathcal{Z}\left\{ n\, u\left[n\right]\right\} =\frac{z^{-1}}{\left(1-z^{-1}\right)^{2}}
$$
What am I doing wrong? Is there perhaps any case in which those properties cannot be applied?

First

$$\frac{d}{dz}\frac{1}{z-1}=-\frac{1}{(z-1)^2}$$

(with a minus sign in front). So you have

$$\mathcal{Z}\{nu[n-1]\}=z\frac{1}{(z-1)^2}=\frac{z^{-1}}{(1-z^{-1})^2}$$

Now everything is OK because

$$nu[n-1]=nu[n]$$

because for $n=0$ also the right-hand side vanishes, and for $n>0$ both are equal anyway.


Why if I place the zero of the lead compensator at lower frequencies, do I obtain a resonance peak? I am studying control systems, and I am studying the lead and lag compensator. I have seen than if I use a lead compensator for the closed loop transfer function,
$T(s)=\frac{s+2}{(s+1)(s+8)},$
with the designed lead compensator,
$lead(s) = \frac{3s+1}{0.1s+1}$
and I use it as a pre-compensator, so that the lead compensator is outside the loop,
$lead(s)\cdot T(s),$
and plot the frequency response, I see that I have a frequency response in which there is a resonance peak, which increases as I decrease the frequency at which the zero is present.

In this case, the red line is a lead compensator defined as,
$\frac{3s+1}{0.1s+1},$
and the green line is the system in which I have used the following compensator,
$\frac{6s+1}{0.1s+1}.$
Can somebody explain to me why?

The first compensator $\frac{3s+1}{0.1s + 1}$ has a high-frequency gain of 30.

The second compensator $\frac{6s+1}{0.1s + 1}$ has a high-frequency gain of 60.

The difference is 6.02 dB, it looks pretty close to the difference between the green and red line.


Time-reversal and delay commutativity and resulting signals I'm reading Schaum's DSP book, and in Fig 1-3 they demonstrate why shifting and reversal are order-dependent, showing a couple of simple systems, 1) delay followed by reversal, and 2) reversal followed by delay.

The equations they present as follows:

1) (Time delay -> Reversal)

$x(n) \rightarrow x(n-n_0) \rightarrow x(-n-n_0)$

2) (Reversal -> Time delay)

$x(n) \rightarrow x(-n) \rightarrow x(-n+n_0)$

I'm fine with shifting and reversal not being commutative due to the different results above, but shouldn't the results be opposite? I.e.:

1) (Time delay -> Reversal)

$x(n) \rightarrow x(n-n_0) \rightarrow x(-(n-n_0)) = x(-n + n_0)$

2) (Reversal -> Time delay)

$x(n) \rightarrow x(-n) \rightarrow x(-n - n_0)$

I didn't find this in the errata for the book which is why I'm thoroughly confused. Which is correct?

Your counterexample to the book's assertion
is confusing between two different uses for $n$.  There was a question earlier
in which some user (endolith? datageist?) gave an answer containing
a detailed description of what exactly
this confusion is and how to interpret the results correctly. My cursory search
has not found this great answer, and so I will just exhibit
a simple example.

Suppose that
$x[n] = \begin{cases}1, &n=0,\\2, & n = 1,\\0, &\text{otherwise,}\end{cases}$

Delaying $x[n]$ by one unit gives
$x_d[n] = \begin{cases}1, &n=1,\\2, & n = 2,\\0, &\text{otherwise,}\end{cases}$

and time-reversing this gives
$x_{d,r}[n] = \begin{cases}1, &n=-1,\\2, & n = -2,\\0, &\text{otherwise,}\end{cases}$  which equals $x[-n-1]$ exactly as the book says it does. Don't believe this?
Here is a table of values for the two functions:
$$\begin{array}{|l|c|c|c|c|c|c|c|c|c|}
\hline\\n&\cdots&-3&-2&-1&0&1&2&3&\cdots\\
\hline\\
x[n]&0&0&0&0&1&2&0&0&0\\
\hline\\x_{d,r}[n]&0&0&2&1&0&0&0&0&0\\
\hline\\
-n-1&\cdots&2&1&0&-1&-2&-3&-4&\cdots\\
\hline\\
x[-n-1]&\cdots&0&2&1&0&0&0&0\\
\hline\\
\end{array}$$

OK, thus far?  On  the other hand,

Time-reversing $x[n]$ gives
$x_r[n] = \begin{cases}1, &n=0,\\2, & n = -1,\\0, &\text{otherwise,}\end{cases}$
and delaying this by one time unit gives
$x_{r,d}[n] = \begin{cases}2, & n = 0,\\1, &n=1,\\0, &\text{otherwise,}\end{cases}$
which equals $x[-n+1]$ exactly as the book says it does. Don't believe this either? Here
is the table for it.
$$\begin{array}{|l|c|c|c|c|c|c|c|c|c|}
\hline\\n&\cdots&-3&-2&-1&0&1&2&3&\cdots\\
\hline\\
x[n]&0&0&0&0&1&2&0&0&0\\
\hline\\x_{r,d}[n]&0&0&0&0&2&1&0&0&0\\
\hline\\
-n+1&\cdots&4&3&2&1&0&1&2&\cdots\\
\hline\\
x[-n+1]&0&0&0&0&2&1&0&0&0\\
\hline\\
\end{array}$$

So how can we establish these results without setting up tables?

Well, $x_d$ is a sequence whose $n$-th term $x_d[n]$ equals $x[n-n_0]$
for all choices of $n$; in other words, take whatever is inside the square brackets
after $x_d$, subtract $n_0$ from it, and stick it in as the argument/index
for $x$. Now, for
any sequence $y$, its time-reversal $y_r$ is a sequence for which
$y_r[k] = y[-k]$ for all choices of $k$ (negate the argument of $y_r$ and
stick it in as the argument for $y$. Thus, the time-reversal of
$x_d$ is given by

$$x_{d,r}[k] = (x_d)_r[k] = x_d[-k] = x[(-k)-n_0] = x[-k-n_0]$$
as the book says, and not $x[-k+n_0]$ as the OP claims it should be.

Similarly, $x_r[n] = x[-n]$ and so
$$x_{r,d}[k] = (x_r)_d[k] = x_r[k-n_0] = x[-k+n_0]$$
as the book says it should be, and not $x[-k-n_0]$ as the OP
claims it should be.


Capacity computation in Costa sheme I can not understand the following sentences, could you please anyone guide me ?
"Comparing with the embedding in audio signals, the usable bandwidth is greatly reduced,
from 22.05 kHz to 4 kHz. According to Costa’s result, the corresponding capacity will be
reduced to l/5 th accordingly."

The costa formula is Cw,Costa=(1/2)*Log2(1+(var_w/var_N))

I also compete it in MATLAB

((1/2)(log(1+(1/4))/log(2)))  /   ((1/2)(log(1+(1/22.05))/log(2)))    = 5.0311
is it right ?

You really didn't give sufficient context and background information to provide a complete answer, but I'll try anyway. You didn't say so, but I assume that $\sigma_W^2$ is the power of a watermark and that $\sigma_N^2$ is the power of additive noise. I furthermore assume that the "usable bandwidth" is the bandwidth for generating the watermark. If we additionally assume that the watermark signal is white with spectral density $W_0$, then its power inside the usable bandwidth $B$ is given by

$$\sigma_W^2=2W_0B$$

If $B_1=22050Hz$ and $B_2=4000Hz$ we get for the capacity ratio

$$\frac{\log_2(1+\frac{2W_0B_2}{\sigma_N^2})}{\log_2(1+\frac{2W_0B_1}{\sigma_N^2})}$$

which is not a fixed ratio, but it depends on the constants $W_0$ and $\sigma_N^2$. If they are not given, then there is no way to claim that this ratio is $1/5$.

However, if you can assume that $\sigma_W^2/\sigma_N^2\ll 1$ then the capacity can be approximated by

$$\frac{1}{2\log (2)}\frac{\sigma_W^2}{\sigma_N^2}=\frac{1}{2\log (2)}\frac{2W_0B}{\sigma_N^2}$$

and with the same assumptions as above the capacity ratio is equal to the ratio of bandwidths $B_2/B_1=0.18$, which is indeed close to $1/5$.

If any of my assumptions are not correct, please clarify your question.


What is a speaker behaviour for out of range frequencies? I was watching this video. The experiment is to create a square signal and record what is outputted by a speaker. We see two diracs signal (+ and -) that looks like the derivative of the signal. Can anyone explain this phenomena? My intuition is that it is caused by some really high frequencies that are out of range for the speakers designed for human. But why the derivative?

Imagine an amplifier/speaker combination with no DC blocking, so the speaker cone moves exactly in proportion to the input voltage.
Now imagine a microphone that measures the actual absolute pressure of the atmosphere, again with no DC blocking.
Even with such a speaker and microphone, because of the finite size of the speaker, the air will act as a highpass filter.
It's complicated, and I'd have to do some research myself to be able to put numbers to it, but roughly speaking the actual response of the sound pressure to the speaker position is going to start rolling off at about the point where the wavelength is the same order of magnitude as the speaker's diameter.  There will be much better information out there on the web, but I suspect you can search it as well as I can.
This is why woofers are bigger than tweeters, and its why subwoofers are constructed so (seemingly) oddly.  Subwoofers, in particular, are significantly smaller than a wavelength of the sound they're producing, so the designers need to play all sorts of acoustic tricks to make them work, and they're horribly inefficient, which is why they often need -- or come with -- their own amplifiers.
When the speaker pushes out suddenly in response to the square wave, it tries to generate a permanent rise in pressure -- but for wavelengths longer than the box, the pressure wave more or less goes around the box and into the hole in the back, and very little sound is produced.  What's left is that spike you see.


Need a reference on period of LFSRs I need a reference for the fact that each LFSR with XOR operation has a period of $2^m-1$ excluding zero state and every LFSR with XNOR operations has the same period excluding all ones state.

I want to cite it when I'm writing about the usage of LFSR in a paper, I've found some references but they avoid a proof and all of them include only the case of XOR operations.

Since LFSRs are mostly used in signal processing I ask it here.

Edit:
I know that the feedback polynomial should be primitive.

Not all LFSRs have the property of maximum cycle length. In fact, those who do have this property are based on a primitive polynomial of order n. It is a basic property of the extension field $GF(2^n)$ which is "spanned" by this primitive polynomial. The LFSR merely iterates through all of the powers of the polynomial primitive root.

Open a book on algebraic structures and cite the relevant part.


Noise estimation SNR matrix I have a signal matrix which is a $256\times 192$, and I want to calculate the SNR considering that my $259\times 192$ matrix  is an average of a $256\times 192\times 330$ matrix, where $330$ is the number of frames, after a reshape I have obtained a $256\times 192$ matrix. So my code is:

signal = 256x192x256;

signal_average = mean(signal,3);
noise_estimation = signal_average - repmat(mean(signal_average,2),1);
signal_power = mean(abs(signal_average).^2);
noise_power = mean(abs(noise_estimation).^2,2);
SNR = 10*log10(signal_power./noise_power);

Is it correct to use this approach in order to obtain a matrix with different SNR for each frame?

In your case you probably want to calculate the SNR as mean over standard deviation.

signal=rand([256,192,330]); %demo data
SNR = mean(signal,3)./std(signal,[],3);
SNRdb = 10*log10(SNR);

this way you obtain different SNR values per pixel. 256x192 pixels in 330 frames.

To get the values for each of the 330 frames instead you must first reshape your matrix.

signal=reshape(signal,[],330);  %pixel values merged
SNR = mean(signal)./std(signal);
SNRdb = 10*log10(SNR);


Basic tools for digital filter design? What tools to use for practicing elementary filter design?
Is MATLAB all there is?
Do I need some specific toolboxes?
What functions do I need?

I'm starting from the ground up in digital filter design and I thought that I need to decide on a program that lets me experiment with different filter designs (not pre-built functions, but rather self-written). Make plots and perhaps .wav outputs and such.

I'm looking for something higher level than C++, because I think the "design language" should allow for more rapid prototyping than the "implementation language".

The tools should facilitate the evaluation of filter designs.

There are lots and lots of software that can aid you in designing digital filters. MATLAB is probably the most used software, at least in the university sector. The  DSP toolbox and the signal processing toolbox probably cover all of the well-known methods for digital filter design, as already mentioned.

Alternatives to MATLAB (that are free) include: Julia, Octave, Scilab, and SciPy (Python with libraries for technical computing). There are a lot of others, but these are the ones that I know have high quality libraries/methods for filter design.

For Julia, you can use the Filters.jl package, for Octave the signal package, for Scilab the Signal Processing toolbox, and for SciPy you have scipy.signal.

A lot of filter design is done using optimization, and if you want to customize the cost-function or make other tweaks to e.g. the least-squares method or the Parks-McClellan method, you could have a look at a high-level optimization library such as CVX (works with MATLAB and Julia ++), or JuMP for Julia.

These are just the tools that I have a fair bit of familiarity with...


Why Is the DFT of a Signal Symmetric About Its Central Bin? When i take an N point DFT of a signal, it comes out to be conjugate symmetric about the point N/2 . Could someone please tell how to understand this intuitively or mathematically ?

First of all it is only conjugate symmetric if the input signal is Real Signal.

The reason for that is the DFT is built by Complex Exponential which can be decomposed into Sine and Cosine.
Since the Sine governs the Complex part and is $ N $ periodic one could see it creates conjugate symmetry of the DFT.

You just need to follow the definition of the DFT.


Differential Encoding for M-PSK Confused about a few topics w.r.t. differential encoding/decoding. First, consider the below BER curve taken from wikipedia ( https://en.wikipedia.org/wiki/Phase-shift_keying ) comparing BPSK/QPSK and DBPSK, DQPSK
I ran some simulations in matlab and found that by including standard differential encoding of my QPSK that my error rate was approx 2x higher than normal QPSK (which seemed correct given the error propagation of the differential coder). However according to this curve, the penalty for DQPSK appears much worse, am I missing something or is there a mistake in this graphic?
EDIT: I'm going to answer this first question - I've conflated terms. The performance in the graph for DQPSK i.e. differential QPSK, is different than differentially encoded QPSK (which is the system I've described in my post). In differential QPSK, I believe the transmitter I described using modulo-M differential encoding is still correct (in that it encodes the information into the phase transitions of the constellation) - however the key difference is that DQPSK at the RX side is assumed to be operating differentially and so detection of the data is done by comparing the current and prev phase of the constellation.
It's clear now my questions are not about DQPSK but instead about differentially encoded QPSK. In this case, we demodulate using a coherent receiver (i.e. including a costas loop) and the differential encoding is only used to resolve the 2*pi/M phase ambiguity. Thus we only get a small BER penalty from the error propagation in the differential encoding.
I'm still interested in the subsequent questions below on the topic of differential encoding schemes.
Secondly, I want to understand two forms of differential encoding for M-PSK that I am aware of. The traditional way as I know it is shown below (assuming perfect timing/freq sync already) where it takes a stream of input bits, groups them into log2(M) chunks, re-maps them via a grey code, differentially encodes them (modulo M), and then maps to symbols (not shown in the flowgraph but the constellation here must be natural, i.e. no mapping). In other words this encodes the data into the transitions of the constellation as opposed to the absolute phase. For QPSK this allows us to correct the 90 degree phase ambiguity problem as seen in the costas loop for example.

The alternate method I am aware of for example with M=4 involves taking the stream of input bits, deinterleaves them (even and odd bit indices) and differentially encodes each stream separately modulo-2, interleaves them back up, groups into chunks of 2-bits, re-maps them via a grey code, and then maps to symbols. This essentially is treating it as two independent DBPSK streams. This is as shown below:

Questions:

Is there a common name for the alternate encoding scheme or other reference material on it? I've found none.
EDIT: I have since found that the interleaved scheme does not appear to handle 90 degree ambiguity (only 180deg) as I had originally thought. In fact I came across a paper that seemed to state as much as well https://ieeexplore.ieee.org/document/1632094 I still am unclear why exactly this is the case however.. I had originally thought that since each stream could separately handle a 180 degree offset, combined they could handle 90 as well.
Is this alternate method acceptable for higher order M-ary PSK?

I don't fully understand what the OP means by DQPSK but suspect that he believes that the D in DQPSK stands for "differentially coherent demodulation of" QPSK, that is, DQPSK means differentially-coherent demodulation of a plain vanilla QPSK signal, and that in his terminology, differentially encoded DQPSK is differentially-coherent demodulation of differentially-encoded QPSK. Most people would take DQPSK to mean the latter because the former interpretation of DQPSK (differentially-coherent demodulation of a plain vanilla QPSK signal) is a really bad idea.
Looking only at the output of a QPSK transmitter, one sees an RF signal whose RF phase has constant value over successive $T$-second epochs, but can change abruptly at the transition from one epoch to the next from one value in $\left\{\frac{\pi}{4},\frac{3\pi}{4}\frac{5\pi}{4}\frac{7\pi}{4}\right\}$ to another value in $\left\{\frac{\pi}{4},\frac{3\pi}{4}\frac{5\pi}{4}\frac{7\pi}{4}\right\}$ (or stay the same, too). For some details, see this answer of mine. Assuming that the signal is a plain vanilla QPSK signal, the (coherent) receiver needs to work hard (using PLLs or Costas loops) to acquire both the carrier frequency and the carrier phase of the incoming signal so that the receiver's local oscillator can be synchronized with the incoming signal and coherent demodulation is possible. With such synchronization, the QPSK system essentially works as two parallel BPSK schemes operating on phase-orthogonal carriers, and the I and Q branches of the receiver make decisions independently as to whether their separate matched-filter outputs correspond to a transmitted $1$ or a transmitted $0$. Thus, for a given value of $E_b/N_0$ (the "bit energy to noise density" or BEND ratio as I like to call it), the BER for coherent demodulation of QPSK is the same as the BER for BPSK. Note however, that the QPSK signal is transmitting $2$ bits in $T$ seconds and so the signal power is doubled.
Next, let us consider  differentially coherent demodulation of the plain vanilla QPSK signal. Here it is not necessary to determine the RF carrier phase -- synchronization in just the carrier frequency suffices --  but all that can be determined is whether the RF phase of the current epoch is the same as the RF phase during the previous epoch, or that the current RF phase differs from the previous RF phase by $\frac{\pi}{2}$ or by $\pi$ or by $\frac{3\pi}{2}$. Now, if we know what the previous dibit was ($00, 10, 11, 10$ as the case may be), we could deduce what the current dibit is. But, this notion is fraught with the possibility of not just error in deciding what the current dibit is, but error propagation. If we get the current dibit wrong, the next dibit will be wrong too because its value will be determined from the wrong value of the current dibit, and the next dibit decision will be wrong, and the next ad infinitum until another dibit error of just the right kind fortuitously restores everything back to normal. Moral: Don't try to do differentially-coherent demodulation of plain vanilla QPSK; it doesn't work well at all regardless of what the OP believes.
If one is doing differentially-coherent demodulation of differentially-encoded QPSK, then the information is being carried not in the dibit values during the epochs, but in the change in dibit value between epochs. For example, the data bits for the current epoch are determined by the change in phase that occurred during the transition between the previous epoch and the current epoch. Such a change of phase is precisely what the differentially-coherent QPSK receiver is designed to detect (cf. the previous paragraph). It is worth emphasizing that differential encoding in QPSK is not the same as doing differential encoding of the binary data streams of the I and Q branches; it is different. For details of thhe differential encoding, see this answer. For more information ow soft outputs from the matched filters are processed in a DQPSK receiver, see here.


$\mathrm{FFT}$ generates non-linear noise I followed the tutorial for generating the frequency spectrum of a signal here, and for testing I generated a sine signal with the frequency of $440 \textrm{ Hz}$. That resulted in the time domain in the following picture:

And when I convert it into the frequency domain, I get:

My question now is:

Why is the background noise shaped like it is?
I would have expected a more even noise, or rather a linear shaped noise. Where does this shape come from?

Edit concerning the technical data:

Sample rate is $44100\,\textrm{Hz}$, which should be enough
FFT length is identical to the signal length, i.e. 1323000 points.

First of all, look at the power of what you call background noise:

It's between -150 dB and -120 dB below your peak's power. Which means that the magnitude of these FFT points is somewhere in $[10^{-7.5};10^{-6}]$ of the main peak magnitude.

That screams "numerical inaccuracy" all over the place:

Assuming your FFT was implemented using 32 bit IEEE754 floating point numbers, you'd have to realize that for each point of the FFT output, there's at least an FFT-length amount of additions taking place with a very finite numerical precision.

Now your FFT has an enormous length; $1323000\approx 1.3\cdot 10^6$ additions bring the most stable algorithm to its knees. For a bit of details on error analysis on commonly used FFT algorithms, refer to [1].

Have a look at the DFT matrix to understand the concept of how these values come to be: For the frequency of your spectral peak, signal-period ($\frac{1}{f_\text{signal}}$) spaced samples are added up; since they always have the same value (definition of "period"), the result is large; for all the others, the should cancel each other out. For example, the bin at $2\cdot f_\text{signal}$ always adds up the same absolute value, but with alternating sign. For other bins, you just take values from the sine and add them up. What happens now is that due to limited floating point math, when you add up a large and a small number, and then subtract the large number again, you won't get the same small number you had – accuracy reduces with magnitude!

Now, the FFT is an algorithm that is pretty resilient with respect to that (because it doesn't try to add up all FFT-length values at once), but there's only so much stability that an algorithm can achieve with a given machine $\epsilon$.

  I would have expected a more even noise, or rather a linear shaped noise.

Why that? If anything, you're observing a process where the imperfections of the FFT are correlated with the input signal. Something like that won't have "white" PSD! Remember, the PSD of a process is the Fourier transform of its autocorrelation; in your case, it seems there's more autocorrelation at higher frequencies.

Do your understanding a favor and add a bit of as-good-as-feasible white, gaussian noise to the input, maybe with a variance of $N=10^{-3}$ of the sine's amplitude, and do your mag squared FFT again. You'll see the noise floor where you expect it, and it will be flat, because the (dominant) source of noise is actually white, and not strongly correlated.

Also, your $1.3\cdot10^6$ FFT really isn't all that helpful, unless you actually need a frequency resolution of $\frac{f_\text{sample}}{N}=\frac{4.41\cdot 10^4}{1.323\cdot10^6}=33\,\text{mHz}$ (which I dare to doubt). Try with shorter pieces of your signal!

[1] Tasche, M., and H. Zeuner, Roundoff error analysis for fast trigonometric transforms, in Handbook of Analytic-Computational Methods
in Applied Mathematics, G. Anastassiou (ed.), Chapman & Hall/CRC, Boca Rota, 2000, 357–406.


Bandpass Filter Implementation in C I've been googling for a while and came across the audio cookbook and Bjorn's blog explaining it further Link Here.

Background info: I've modified the rtl-sdr open source code to tune to a VOR station (transmits heading, typically used for aircraft). I tune to 113.1 MHz which is the the frequency the nearby VOR station transmist data. Inside that 113.1 MHz frequency, there are 3 signals;

Morse Code at 1020Hz
Reference Signal (9960Hz +- 480Hz) (FM this signal, inside this there is another 30Hz signal)
Variable Signal (30Hz)

I'm trying to isolate these signals. Now when I tune to the station (113.1MHz) and perform amplitude modulation, I can hear the Morse code signal from the speaker with lots of noise in the background. I followed the audio cookbook guide along with Bjorn's implementation to write a BPF function (pasted below) to isolate the Morse code signal (1020Hz) to hear only the morse code with very little to no noise in the background. However, I just about hear nothing, it's like someone just really lowered the volume almost to mute. I'm starting with this because if I can actually hear the morse code signal, then that proves that the BPF works and now I can use that function to isolate the other two signals.

The values I used to compute a0 to b2: I also did normalize the a0 to b2 values.

Sample Rate: 24KHz
Center Freq: 1020Hz
BW (in octaves): 1 (don't really understand this)
gain: 1 (assuming 1 should be ok but again not really sure)

BPF Function in C Language:

int band_pass_fir(struct morse *_m, int16_t *_signal, int _len)
{
  // _m->signal (int16 array) is my desired output
  // _signal (int16 array) contains the actual amplitude modulated data received from the signal
  // _len is the length of _signal

  int i = 2;

  const float b0 = 0.1205498139,
              b1 = 0,
              b2 = -0.1205498139,
              a1 = -1.7626236142,
              a2 = 0.8273910712;

  while (i < _len)
  {
    _m->signal[i] = (int16_t)(b0 * _signal[i]) +
                    (int16_t)(b1 * _signal[i-1]) +
                    (int16_t)(b2 * _signal[i-2]) -
                    (int16_t)(a1 * _m->signal[i-1]) -
                    (int16_t)(a2 * _m->signal[i-2]);

    _m->signal[i-2] = _m->signal[i-1];
    _m->signal[i-1] = _m->signal[i];

    i++;
  }
}

Questions:

To isolate this 1020Hz signal, is it safe to assume I need a BPF? or do I need some other type of filter?
Assuming whatever type of filter I need to isolate the 1020Hz would be the same type of filter I would need to isolate the 30Hz signal (variable signal) ? Is that correct?
Do you see if I'm doing anything wrong with my above method?

There are a few potential issues:

Your code doesn't compute a full output vector, you just keep the last three samples (often referred to as "state variables". Your output should be a vector of the same length as the input
This is only a second order filter. Depending on the exact properties of your signal, this  may not be enough to isolate the component that you want. I'd recommended a Butterworth bandpass filter. In this case you can independently control the lower cutoff and the higher cutoff of your bandpass and also the steepness of the slopes. There are butterworth design tools in "Octave", see http://octave.sourceforge.net/signal/function/butter.html
If you need to go to higher order filters, they should be implemented as "cascaded biquads" or "cascaded second order sections". There are tons of code examples to be found through google


How to smooth the echo duration change to get rid of audio clicks? I implemented an echo effect using the following diagrams:

To summarize, an echo effect is basically a positive feedback effect that is implemented with a delay line zD which delays the input signal by D samples. My effects works fine. I get the echo effect.
My issue is I gave the user to adjust D sample amount during runtime just like in Pioneer equiments you can set echo as 4beats, 2beats, 1beats etc. When a user request a beat duration change I immediatly change the D parameter in my audio plugin. The problem is sometimes this causes click due to sudden audio level changes. I am trying to find a way to counter this but so far I can not think of anyways to do it. How may I ensure that the audio level is continous(doesnt spike and result in audio clicks) when user changes the parameter D?

Maybe, when $D$ changes, you can crossfade from the audio at the old $D$ to audio at the new $D$.
If you ran a pitch detector alongside of the delay (and stored the detected period length in a corresponding array) you could crossfade one period at a time until you got to within a period of your new $D$.  Then slide very slowly to the exact $D$.


What's the minimum decibel value? I'm trying to scale wave amplitude (0-1) to perceived loudness (still 0-1). The decibel scale seems to be unbounded. For my audio, the amplitude 1 represents a decibel $d$ - i.e. the audio is being played at volume $d$. But what would the decibel be for 0 amplitude? I know 0dB is the minimum a human can hear, but what's the actual minimum. If I can obtain this, then I can successfully convert wave amplitude to perceived loudness (while staying between 0 and 1)

By any reasonable metric, the dB value for zero amplitude is minus infinity.


Resample an array in real-time (application: pitch-bend) I'm working on a sampler software (that plays .WAV files, when notes are played on a MIDI keyboard).

In order to implement pitch-bend feature (i.e. you have a "pitch bend wheel" on the synth and you can turn this wheel to retune the sound in realtime, classical effect used in funk music),

we have to:

Retune: this is possible by resampling the array of the .WAV.
With Python, there is scipy.signal.resample or probably the same in numpy.
But I heard there are algorithms optimized for audio repitching. What are the audio-optimized resampling algorithms?
[Secret Rabbit Code : http://www.mega-nerd.com/SRC/. Why is this important? How to implement it in Python, with numpy for example?]
Retune in realtime: indeed, we can play with the pitchbend-wheel very fast, thus, hundreds of resampling factor could be needed per second! How to do that?
How to resample an array with different resampling factors, in realtime?

Audio-optimized resampling is done by using anti-alias filtering and interpolation kernels with good stop-band attenuation and low pass-band ripple/distortion.  On current PC and mobile processors (the kind that can easily compute hundreds of transcendental functions per sample time), these filter and interpolation kernels can not only be precomputed and cached in poly-phase lookup tables, but for certain window types, (re)computed in real-time per sample.

Thus, resampling can be done in real-time by computing the needed interpolation kernels and interpolating each sample as needed at the required sample delay from the last for the current sample rate.

Note that changing the sample rate too quickly may result in FM modulation skirts that may need additional filtering to reduce aliasing artifacts.


Purpose of Phase Information I am learning Fourier Transform from many days but till now I am not able to understand what does phase angle image show us or tell us? They say that MAGNITUDE tells "how much" of a certain frequency component is present and the PHASE tells "where" the frequency component is in the image. I understood about the MAGNITUDE part but not able to understand PHASE part. Please explain what it means by "PHASE tells "where" the frequency component is in the image". Thanks

I don't think of the phase as providing information related to position or location.

A simple example might be illustrative. In 1-D, consider these two signals, x and y:

t = 0:0.001:0.1;
x = cos(2*pi*20*t) + cos(2*pi*30*t);
y = cos(2*pi*20*t + pi/4) + cos(2*pi*30*t - pi/6);

The magnitude spectrums of x and y are the same, but they're not the same signal; if you plot them, you'll see that they're different. In the frequency domain, you need both the magnitude and the phase to fully describe the signals. It doesn't make sense to refer to the harmonics as "being somewhere".


3 dB cut-off frequency of first-order IIR high-pass filter As descipbed in this great answer (thank you Matt L.), the suggested exact formula was very precise.
I need a first order high-pass with similar precision.
EDIT: (clarification)
The low-pass filter discussed in the link above was this one, derived from the simple RC filter.
What I'm looking for is the discrete version of the simple high-pass RC filter with a formula for getting the exact alpha term (that's what Matt L. did) for an amplitude response of -3dB at a given cutoff frequency.
I've tried and measured the low pass version and it worked as expected.

A discrete-time first-order high pass filter with unity gain at Nyquist and a zero at DC is described by the following difference equation:
$$y[n]=\frac{1+\alpha}{2}\big(x[n]-x[n-1]\big)+\alpha y[n-1],\qquad -1<\alpha<1\tag{1}$$
Its transfer function is given by
$$H(z)=\frac{1+\alpha}{2}\frac{1-z^{-1}}{1-\alpha z^{-1}}\tag{2}$$
Evaluating the squared magnitude of $(2)$ on the unit circle $z=e^{j\omega}$ and equating it to $\frac12$ ($-3$ dB) results in the following relation between $\alpha$ and the $3$ dB cut-off frequency $\omega_c$:
$$\begin{align}\big|H(e^{j\omega_c})\big|^2&=\frac{(1+\alpha)^2}{4}\frac{\left|1-e^{-j\omega_c}\right|^2}{\left|1-\alpha e^{-j\omega_c}\right|^2}\\&=\frac{(1+\alpha)^2}{4}\frac{2\big(1-\cos(\omega_c)\big)}{1-2\alpha\cos(\omega_c)+\alpha^2}\stackrel{!}{=}\frac12\tag{3}\end{align}$$
Eq. $(3)$ results in a quadratic equation for $\alpha$ with the solution
$$\alpha=\begin{cases}\displaystyle\frac{1-\sin(\omega_c)}{\cos(\omega_c)},&\omega_c\in(0,\pi)\setminus \frac{\pi}{2}\\0,&\omega_c=\frac{\pi}{2}\end{cases}\tag{4}$$
(where the requirement $|\alpha|<1$) has been taken into account).
For $\omega_c=\pi/2$ we obtain $\alpha=0$ and the corresponding filter is a simple $2$-tap FIR filter. All other cut-off frequencies $\omega_c\in(0,\pi)$ result in IIR filters.
The figure below shows the magnitude responses of $9$ high-pass filters with specified cut-off frequencies $0.1\pi,0.2\pi, \ldots,0.9\pi$. The corresponding values for $\alpha$ were computed according to Eq. $(4)$.


Reading the Wavelet transform plot I am having trouble understanding on how to read the plot plotted by a wavelet transform,

here is my simple Matlab code,

load noissin;
% c is a 48-by-1000 matrix, each row
% of which corresponds to a single scale.
c = cwt(noissin,1:48,'db4','plot');

So the brightest part means the scaling coffiecient size is bigger, but how exactly i can understand this plot what is happening there ?
Kindly help me.

Imagine for one second, that you just plotted your daubechies-4 wavelet, as you can
see here in red.

Now imagine that you take this waveform in red, and simply do a cross-correlation with your signal. You plot that result. This will be the first row of your plot. This is scale-1. Next, you dilate your Daubechies-4 wavelet, (that is, you simply make it 'stretch' it in time, by some factor). Then, you again do a cross-correlation of this new waveform with your signal. You then get row two of your plot. This is scale-2. You keep doing this for all scales, which means you keep taking your original 'mother' wavelet, and you keep dilating, then cross-correlating, dilating, and cross-correlating, etc, and you just plot the results one on top of the others.

This is what the CWT plot is showing you. The results of performing a cross-correlation of your signal with a wavelet at different scales, that is, at different dilation (stretch) factors.

So let us interpret your image. In the first row, you can see that you have weak amplitudes in your cross-correlation. That means that it is telling you, almost nothing in your signal correlates, (or 'matches') your wavelet, when it is at scale-1, (when it is at the smallest scale). You keep stretching your wavelet and correlating, and still it is not matching anything in your signal, until you reach say, scale-31. So when you stretch your wavelet 31 times, and performed a cross-correlation, you start to see some bright spots, which means you are getting good cross-correlation scores between your stretched wavelet, and your signal.

If you look however at the very top, we have the brightest spots. So for scale-46, you made that row by stretching your original wavelet 46 times, and then cross-correlated it with your signal, and then that is your row-46. So you see a lot of nice bright spots. You can see that at positions (x-axis) ~25, ~190, and ~610, I have bright spots. So that is telling you, you have some feature in your signal, that very closely matches your wavelet that is stretched 46 times. So you have 'something' at those locations that closely match your wavelet at this scale.

(Of course, in your case, you have used noise, so the positions that I have talked about are random - that is, there is nothing really 'interesting' going on. Do a CWT with a sine pulse and what I am saying can be made clearer to you.)

In summary, the CWT is simply showing you all possible correlation scores between your template/matched filter (in this case daub-4 wavelet), at different positions, (x-axis) , also at different stretched factors, (y-axis).

Hope this helped.


Fast way to calculate a few points in convolution? Say I have two 1000-point sequences. I can compute the entire convolution fairly efficiently using FFTs (keeping in mind that zeroes have to be padded on to avoid doing a circular convolution).

However, I'm only interested in about 500 points in this convolution. To be more specific, say the indices of the convolution are 0 to 1998; I'm only interested in 1000 - 1499 inclusive.

Is there some way to calculate only these points, without calcualting the entire convolution and throwing most of it away?

I'll write names of time domain sequences in lower case and their frequency-domain counterparts in upper case, for example $X = \operatorname{DTFT}(x)$, where $\operatorname{DTFT}$ stands for discrete-time Fourier transform. I chose this transform, because it is unaffected by zero-padding of the time-domain sequences, allowing equations that mix frequency-domain counterparts of sequences of different lengths. With the upper-case frequency-domain symbols, time-domain convolution is denoted by multiplication, and a time-domain delay of $N$ samples is denoted by $z^{-N},$ where $z = e^{i\omega},$ with $\omega$ as continuous frequency. No frequency-domain thinking is required here. You can mentally translate the equations while reading to time-domain convolutions, delays, and additions. Addition in frequency domain is addition in time domain.

The frequency-domain counterparts $X$ of $Y$ of time-domain sequences $x$ and $x$ each of length $2N,$ with $N=500,$ can be written in partitioned form $X = X_0 + X_1z^{-N}$ and $Y = Y_0 + Y_1z^{-N},$ where $\{X_0, X_1, Y_0, Y_1\}$ are frequency-domain counterparts of time-domain partitions $\{x_0, x_1, y_0, y_1\}$ each of length $N.$ The partitions are non-overlapping sub-sequences of $x$ and $y$, stripped of delay which is then reintroduced in the equations. Convolution of $x$ and $y$ can be written in frequency domain as:

$$\begin{align}&XY\\
= &(X_0 + X_1z^{-N})(Y_0 + Y_1z^{-N})\\
= &(X_0 + X_1z^{-N})Y_0 + (X_0 + X_1z^{-N})Y_1z^{-N}\\
= &X_0Y_0 + X_1Y_0z^{-N} + X_0Y_1z^{-N} + X_1Y_1z^{-N}z^{-N}\\
= &X_0Y_0 + (X_1Y_0 + X_0Y_1)z^{-N} + X_1Y_1z^{-2N}\end{align}\tag{1}$$

Now all convolutions between sequences are convolutions between partitions of length $N$. Each result of such a convolution is of length $2N-1$. Denoting the frequency-domain counterpart of the first $N$ samples of the convolution result by $H()$ as in "head" and the frequency-domain counterpart of the delay-stripped sub-sequence formed by the remaining $N-1$ samples by $T()$ as in "tail":

$$\begin{align}= &H(X_0Y_0) + T(X_0Y_0)z^{-N}\\
+ &\big(H(X_1Y_0) + T(X_1Y_0)z^{-N}\big)z^{-N} + \big(H(X_0Y_1) + T(X_0Y_1)z^{-N}\big)z^{-N}\\
+ &\big(H(X_1Y_1) + T(X_1Y_1)z^{-N}\big)z^{-2N}\\
= &H(X_0Y_0)\\
+ &\big(T(X_0Y_0) + H(X_1Y_0) + H(X_0Y_1)\big)z^{-N}\\
+ &\big(T(X_1Y_0) + T(X_0Y_1) + H(X_1Y_1)\big)z^{-2N}\\
+ &T(X_1Y_1)z^{-3N}\end{align}\tag{2}$$

Now we have neatly split the full convolution result into length $N$ partitions. You were interested in the partition of length $N$ beginning at index $2N,$ or in frequency domain:

$$T(X_1Y_0) + T(X_0Y_1) + H(X_1Y_1)$$

The time-domain counterparts of these three terms can can be calculated by three FFT convolutions of length $2N$ compared to your single FFT convolution of length $4N.$ I don't know if that is of any use. You'd still be throwing away a total of $3N$ samples of convolution results, as originally.

But wait! The time-domain counterpart of $T(X_0Y_1) + H(X_1Y_1)$ can be calculated by a circular FFT convolution of length $2N$ between $x$ and $y_1.$ First, writing that as linear (not circular) convolution:

$$\begin{align}&XY_1\\
= &X_0Y_1 + X_1Y_1z^{-N}\\
= &H(X_0Y_1) + T(X_0Y_1)z^{-N} + H(X_1Y_1)z^{-N} + T(X_1Y_1)z^{-2N}\\
= &H(X_0Y_1) + \big(T(X_0Y_1) + H(X_1Y_1)\big)z^{-N} + T(X_1Y_1)z^{-2N}\end{align}$$

Wrapping that around to make it circular, $z^{-2N}$ wraps to $z^{0} = 1:$

$$H(X_0Y_1) + T(X_1Y_1) + \big(T(X_0Y_1) + H(X_1Y_1)\big)z^{-N}$$

The second half of that has what we wanted, $T(X_0Y_1) + H(X_1Y_1).$ We still need to calculate the time-domain counterpart of $T(X1, Y0)$ by another length $2N$ FFT convolution so we are down to two FFT convolutions of length $2N,$ with only a total of $2N$ samples of the results wasted.

We can also confirm that user Fat32's suggestion of using a length $3N$ FFT convolution works, by wrapping the latter part of Eq. 2 around whereby $z^{-3N}$ becomes $z^{0} = 1$:

$$\begin{align}&H(X_0Y_0) + T(X_1Y_1)\\
+ &\big(T(X_0Y_0) + H(X_1Y_0) + H(X_0Y_1)\big)z^{-N}\\
+ &\big(T(X_1Y_0) + T(X_0Y_1) + H(X_1Y_1)\big)z^{-2N}\end{align}$$

The middle line has what we wanted. So a single length $3N$ FFT convolution is also an option, with only a total of $2N$ result samples wasted.


Phase delay between two signal using FFT I am working in calculating the phase delay between two signals. The signals are sent through the same transmitter and received in two different receivers.
They have many frequencies, but with a primary one. I need to calculate the phase delay at the maximum value. The FFT magnitude gives me the main frequency component. I can calculate the phase difference at any exact frequency bin around the max freq if I use the same freq for both signals.
B
elow are the FFT magnitude and phase.

As you can see the phase at maximum freq is always located at a peak. Any small change because of noise gives a big error.
When I unwrap the phase the difference is very big.
Is this the correct approach? Any idea of how to deal with this?
Thank you

Do an fftshift before the FFT (rotate your sample vector by N/2), this relocates the phase reference to the center of your sample window (instead of at a potential circular discontinuity around the ends).  This also often requires less unwrapping, as well as making phase interpolation (between FFT result bins) and phase differences more intuitive.


An intuitive explanation as to why a matched filter is time reversed? It easy enough to study correlation, and matched filters. But the challenge I see unmet anywhere to date and which I struggle to meet myself is a simpler presentation, to with more difficult task I suspect.
Can you explain, in lay terms, to satisfy the intuitions of a listener, why the matched filter is a time reversed copy of the signal we hope to detect? It is not a mathematical proof, or explanation I'm after (they abound) but other some way of explaining it (to lay listeners) that generates a "a ha" experience, a feeling that it makes sense, and I get why we time reverse the signal for a matched filter ...
This is something I have yet to find, read or master myself.

Picture the transmitted signal as a «signature». You want to find some process that maximize the probability of detecting that signature even when there is noise.
What do you do to find some signature buried in noise? You make a template that exactly match the known signature, and you slide it back or forth in time, noting how much the actual signal deviates from the ideal template in any one spot. The time-shift that gives a «large enough» correspondence is your assumed signal location.
The time-reversal of the template is just about reversing the reversal of that argument in the convolution operation. Think of it as non-reversed correlation if you prefer (leaving out complex numbers right now).
-k


How to compute an old style shore line drawing? I wish to represent a shore line, old-style, as in the following image (city on the left, water on the right):

UPDATE: new image for another example of what I'm after, namely the thin lines (given the thick lines). [Credit: https://www.atlasobscura.com/articles/early-map-of-chicago ]

My input is an ordered list of XY points that represent the shore line. I have no clue whatsoever how to compute the additional lines. I have trouble finding resource on the subject through search engines since I'm quite unsure of the correct English terms to describe such an operation.

I've thought of applying a simple dilatation, but this is not the way to go here (in particular, it does not work for concave parts of the shore line).

Any ideas?
(and any ideas how to tag this question?)

Define a level function curve (L) like this:

$$ L(x,y) = C_0 \cdot x + C_1 \cdot y + C_2 \cdot x^2 + C_3 \cdot xy + C_4 \cdot y^2 .... $$

It has no constanst term.  Now, find the best fit to the level curve $L(x,y)=1$.  This will yield the values for the coefficients.

Normalize the coefficients so the gradient is of a magnitude you like.

Find the level curves of $L(x,y)=t$ where $t$ is near one.

Finding the best fit is done with Linear Algebra.  Finding the level curves is an interative search algorithm.


How can I calculate the frequency using different color bands? Let's say that I have a 1 dimensional signal that has data which comes from one of 2 different spectral bands. The measurements come evenly mixed between the two bands. I'm interested in a method of extending the frequency calculation by combining information from both signals.

To show an example, let's just say that I have a Red band and a Green Band. The detector arrays would look like this: RGRGRGRGRG. Let's say that I had 512 R, and 512 G, just to give a concrete example. As the detectors are neighboring each other, and most signals have some component of red and green, it seems likely that a common frequency calculation could be made by combining the bands, and possibly double the frequency range which can be detected. How could I best go about doing this?

From a practical perspective, there's good news and there's bad news.  The good news is that, in a sense, consumer grade digital cameras [at least] already do something like what you're describing, so there's a large amount of information available about the process, which is called demosaicing.  The bad news is that, in "standard" designs, the various RAW color channels actually have lower resolution than the physical sensor array, and the color channels in non-RAW images have already been reconstructed.

The reason for this is that the physical sensors only detect light intensity (not color), and therefore each available physical "pixel" is paired with a color filter, according to some pattern. For example, the Bayer Filter is a common pattern which allocates 50% of the available spatial resolution to the green channel, but only 25% each to the red and blue channels.  The various spatial patterns are called Color Filter Arrays (CFAs).  In effect they model the idea of separate RGB sensors, just in a more cost-effective manner.  If there is some way of exploiting "spectral crosstalk" from the materials used in a given CFA, it would be handled in the context of demosaicing.

That means, if you're dealing with a camera which uses that type of sensing scheme, the most reasonable thing you can hope for is access to the camera RAW data (i.e. un-demosaiced), from which you can apply a higher-quality demosaicing algorithm to reconstruct the image.  From the wiki:

  When one has access to the raw image data from a digital camera, one
  can use computer software with a variety of different demosaicing
  algorithms instead of being limited to the one built into the
  camera...The differences in rendering the finest detail (and grain
  texture) that come from the choice of demosaicing algorithm are among
  the main differences between various raw developers; often
  photographers will prefer a particular program for aesthetic reasons
  related to this effect.

The following link seems to be a good overview of various reconstruction algorithms:

A Study of Spatial Color Interpolation Algorithms for Single-Detector Digital Cameras

In higher-end cameras, however, the RAW color channels might not even be spatially interleaved.  Two of the different higher-end schemes (beam splitting, temporal interleaving) are reviewed here.

One type of sensor that I think works more like what you had in mind is the Foveon X3 Sensor. The following paragraph from the link gives a good overview of the resolution issues:

  For example, the dimensions of the photosite array in the sensor in
  the Sigma SD10 camera are 2268 × 1512, and the camera produces a
  native file size of those dimensions (times three color layers). This
  amounts to approximately 3.4 million three-color pixels. However, it
  has been advertised as a 10.2 MP camera by taking account of the fact
  that each photosite contains stacked red, green, and blue color
  sensing photodiodes, or pixel sensors (2268 × 1512 × 3). By
  comparison, the dimensions of the photosite array in the 10.2 MP Bayer
  sensor in the Nikon D200 camera are 3872 × 2592, but there is only one
  photodiode, or one pixel sensor, at each site. The cameras have equal
  numbers of photodiodes, and produce similar RAW data file sizes, but
  the Bayer filter camera produces a larger native file size via
  demosaicing.

I think this is what you're looking for, assuming I understood the question.


Minimum SNR Requirements for Maximum Likelihood Frequency Estimation In certain applications, you have enough SNR available to, for example, perform an FFT and identify peak location and hence the signal frequency. If my understanding is correct, parameter estimation techniques, such as maximum likelihood, will only be useful when your SNR is so low that you can't run a basic peak search to identify frequency.

However, I have seen literature where people plot performance of various estimation algorithms vs SNRs upto 40~50dB. It seems like a waste of processing power to perform ML on a high SNR (40dB is quite a high value) signal and extract the already obvious information?

Is there a 'rule of thumb' to say, use statistical algorithms below such SNR levels?

I doubt that there is such a rule because your assumption- the negative of computational work more than offsets the benefits of lower error probabilities- assumes that computations are expensive and errors are not. The "costs" of computational work and errors depends entirely on the application.


Train Cascade Object Detector example in MATLAB I am trying to use the example from here to detect a stop sign.

But when I change the image to search for stop signs in there are a huge number of false positives.

Am I doing something wrong?

You are not doing anything wrong.  This particular example is there for a quick illustration, and it does not produce a very good detector. It trains a 5-stage detector from a very small training set.  A decent detector may have 20 stages or more, and you would need thousands of positive samples and negative images to train it.

Take a look at this tutorial for more information.


Clean noise in neural data - MATLAB I'm using Matlab to process neurophysiological data and i'm currently dealing with this sort of noise in my signal:

Close up

EDIT:

For clarification, the clean signal should look like this:

EDIT 2:

Clean

Noise contaminated

These spectra were taken from 10s segments as the whole block is considerably bigger. Both spectra are from segments that include the data in the first and third figures (with and without noise bits of data). The noise events are sparse within each data block.

This happens somewhat frequently across datasets and due to it, i'm having some trouble detecting events of interest (i.e. neural spike detection). As i'm trying to detect small amplitude spikes in the signal, it's critical to preserve as much of the original signal as possible.

Any suggestion on how to better approach this problem? What would be te best way to correct these shifts/spikes in the signal, keeping as much of the signal as possible?

Any suggestion/help will be immensely appreciated!

Cheers, M

Following, a way to remove the 3 exponential pulses.

1.- Rebuilding signal + interferer from graph provided

clear all;clc;close all

A=imread('001.jpg');
A1=A(:,:,1);
figure(1);imshow(A1)

A2=~imbinarize(A1);
figure(2);imshow(A2)
% 001

problem: bwmorph breaks causality and splits signal

A3=bwmorph(A2,'skel',Inf);
figure(3);imshow(A3);
% 002

[sz1 sz2]=size(A2)

for k=1:1:sz2
    L1=A2(:,k);
    n1=find(L1);

    A2(:,k)=zeros(1,sz1);
    if ~isempty(n1)
        A2(floor(mean(n1)),k)=1;
    end
end
figure(4);imshow(A2)
% 003

Now signal is correctly mapped.

Since the contaminated signal is available to the question originator, the original data should be inserted here:

s=0
for k=1:1:sz2
    L1=A2(:,k)
    s=[s find(L1)];
end
s(1)=[]

2.- where is mean(s)?

don't do just s=s-mean(s)

Left behind the y=0 axis tick on purpose

% 004

mean_s=287

s=-(s-mean_s)

3.- Rebuilding time reference

t=linspace(0,10,sz2)
dt=mean(diff(t))

4.- Plotting signal + interference

figure(5);plot(t,s);grid on
xlabel('t[seconds]');ylabel('s+i')  % signal+interferer

Comments:

The VLF bending within the initial 5 seconds may come from conducted interference; contact, applying an ECG sensor on a mobile phone?

The base station or mobile seeking answer, starts with certain power and
then as the packets are exchanged Base Station Mobile Terminal tell one another to increase
or reduce power, usually the Base Station starts higher than needed
to make sure reaching cell edge, and then reduces power until MS says
level ok.

The negative sharp spikes show there may be something metallic reflecting the interfering signal.

You have provided a reference signal, what it should be like, but there's
no time reference, just the amount of samples, without further
information, it's not sure the reference signal duration is also 10 seconds, is it?

[pks,locs]=findpeaks(s,'MinPeakHeight',50,'MinPeakDistance',50)
s_peaks=s(locs)

5.- When do the 3 exponential decay interfering pulses start?

t_peaks=t(locs)

6.- What is the mean delay between those 3 pulses measured in time?

diff(t_peaks)

7.- What is the mean delay between those 3 pulses measured in amount samples?

diff(locs)

8.- What is the real amplitude of the pulses? approximately the 1st sample

ds_peaks=s(locs)-s(locs-10)

9.- what is the amount of samples per interfering pulse?

nT1=floor(mean(diff(locs)))

10.- What are the respective attenuation constant of the 3 exponential pulses?

a=-1./(dt*nT1)*log(1./ds_peaks)

11.- 1st interfering pulse start time measured directly on the graph

t0=.688

12.- What is the sample numeral of the start moment of the 1st interfering pulse?

nt0=[1:1:floor(mean(intersect(find(t>t0-.01),find(t<t0+.01))))]

13.- Values of each interfering pulse:

e1=ds_peaks(1)*exp(-a(1)*dt*([1:1:nT1]+t0))
e2=ds_peaks(2)*exp(-a(2)*dt*([1:1:nT1]+t0))
e3=ds_peaks(3)*exp(-a(3)*dt*([1:1:nT1]+t0))
e_interf=[zeros(1,numel(nt0)) e1 e2 e3 zeros(1,sz2-3*nT1-numel(nt0))]
hold all
plot(t,e_interf)

s2=s-e_interf
figure(6);plot(t,s2);grid on
xlabel('t[seconds]');ylabel('s+i-3pulses')  % signal+interferrer

If interested I can show how to use findpeaks to remove the sharp negative pulses.


Blending Artifacts in Photo Stitching I am working on a photo stitching application that uses multi-band blending. I need to get rid of unpleasant edges appearing at some places:

Here is the area of overlap (left - new image added to the mosaic, right - current mosaic contaning pixels of new image on background pixels to improve blending, middle - blending mask):

If I just compute weighted average between left and right image according to mask, the result is of course OK. However, this would leave a visible seam as two images have usually slightly different exposure.

So all three images need to be successively blurred to build a Gaussian pyramid - here is how one level of the pyramid looks like:

You can see that the top part of blending mask "touches" the border. The Gaussian blurring filter reflects on image borders and this causes inaccuracy in lower band.

I colored the images to make the problem more visible:

I am not sure about how to blur the mask so that it would deal with overlap area edges nicely.

Few suggestions:

change behavior of the blurring filter (how?)
extend area of blending mask so that one part never "touches" the border
update weighting masks so that resulting blend mask is more "edge-aware"

Any other suggestions/hints?

I have solved this problem by adding certain padding to overlap area.

There are pixels belonging to image1, image2, background and overlap. A pixels in overlap are successively relabeled to either image1 or image2 depending on neighborhood.

This padding will make some space for the blurring so that sharp edges of the overlap have no change of appearing.

Another treat is photometric calibration, i.e. gain compensation and vignetting removal. This minimizes or even eliminates the differences in low frequency band.

Finally, the blurring have to avoid background pixels (usually black) from bleeding into the image. This can be done by alpha blending or using binary mask as described here.

It is also possible to reduce amount of blurring near edges and corners of the overlap. This, unfortunately, reduces performance of the blending algorithm and the results are somewhat uncontrolled (as the blurring is non-uniform).

Finally, a gradient domain blending improves the quality near edges as it does not rely on Gaussian blurs of fixed sizes.


Calculating gain coefficients for comb filters in Moorer's algorithm Could someone explain how Moorer calculated the gain coefficients for the comb filters in his reverb? In his book, Udo Zölzer's layed out this equation:

$$g = 10 ^ {-3 \left(\mathrm{Reverb \ Time} \cdot \frac{\mathrm{Sampling \ Rate}}{\mathrm{Delay \ length \ in \ samples}}\right)}$$

In his paper for the first comb filter he gives $50 \ \mathrm{ms}$ for the delay and the gain coefficient sampled at $25 \ \mathrm{kHz}$ to be $0.24$ given that reverb time is about $2$ seconds.

So if I plugged in those values I don't get near $0.24$,
$10 ^ {-3 \frac{2 \cdot 25000} {50 \cdot 25}}$.

The values in Moorer's paper were not determined using a formula but they were computed by numerical optimization:

  There is a question, however, of how to set this parameter $g_1$ which controls the roll-off of the filter. To this end, a series of optimizations was done using the Marquardt algorithm (Marquardt, 1963) to match the frequency response of this low-pass filter in a least-magnitude-squared manner to the actual data reported in the literature.

from: J.A. Moorer. About This Reverberation Business. Computer Music Journal, vol.3, no.2, pp. 13-18, 1979.

The formula in Zölzer's book tries to approximate those values from the given parameters. I'm not sure who came up with that formula, but in any case, it doesn't seem to be accurate. E.g., note that for a fixed delay length (not in samples, but in seconds), the formula becomes independent of the sampling frequency, but in Moorer's article the values for $g_1$ do depend on the sampling frequency. Also note that Zölzer's book is not known for its accurate formulas, as shown in this answer to a question about another formula in the book.

In sum, I do not think that the formula is useful, and you better use the values in Moorer's paper as a starting point for your own experiments.


How to recognise partial text on mechanical indicators I would like to know is there any possible approach to recognise digits of the first (and second as well, as you can see in the example figure) order of magnitude. I've thought of neural network with full, half and possibly quarter training samples of the digits I want to recognise, but I don't think this is the best approach. You see, I am new to computer vision and I'd like to see what path should I take to study if there are any for this problem.

If you are using an out-of-the box OCR engine, then there's not much you can do with significantly occluded letters.

However, if you are rolling your own OCR, say with some template-matching, then, there is nothing preventing you from learning/teaching partial letters and their valid combinations (like the 8-9 and 9-0 above). Just treat it as another letter in the set to be recognized.


Trying to convert Raw image(Bayer format) to jpg/png image in Python I am using the Sony Low-light image dataset from the paper Learning-to-see-in-the-Dark. I am trying to convert Raw image (bayer pattern in arw format) into a jpg/png. I used an online tool that gives a very dark jpg for a short exposure image and I have also used Python rawpy library to read and process the raw image. The resulting image from rawpy is shown below:

Rawpy gives a brightened image for a low-light image, which is not desired. I was expecting a dark rgb image. How can I get a dark rgb image in jpg format using rawpy?

Did you solve your problem?

Otherwise, I can introduce you an utility software tool (free to download) to convert any raw images to jpgs. I converted all long/short .RAF and .ARW images (in Learn to see in dark) to jpgs in just 30 sec. Unlike rawpy, it doesn't change (brighten) the images.


What are physical scenarious which can be modeled as multiplicative noise? Consider $s(t)$ as the signal, $w(t)$ as the noise, and $y(t)$ as the captured signal. Usage of additive noise model, that is $y(t)=s(t)+w(t)$, is quite wide spread, example in speech and audio signal enhancement, and is also physically convincing to hypothesize. The convolutive noise model,that is $y(t)=s(t)\star w(t)$, where, delayed copies of the signal overlap together is another widespread model. A physically relevant scenario for this is reverberation modelling.

My question is - Is there a physical motivation to model noisy signal using multiplicative noise model, that is $y(t)=s(t)w(t)$? I am unable to think about any convincing physical domain examples where this happens in reality.

In communications, a very important multiplicative model is a flat fading channel where the received signal can be modeled as

$$y(t)=h(t)\cdot s(t)+n(t)\tag{1}$$

where $h(t)$ is the time-varying attenuation caused by the channel, and $n(t)$ is additive noise.


Two voice pronunciation comparison similarity MFCC + DTW Calculate each MFCC to compare wave file A and wave file B, and then use FastDTW to measure the distance after two sets of MFCCs.

We compared the four wave files and obtained the Euclidean distance value.

The values below are the Euclidean distance values.

  675.0095954620155 A.wav vs. A2.wav

  998.7554375714773 A.wav vs B.wav

  976.1293903229977 A.wav vs. B2.wav

  856.4364672719398 A.wav vs C.wav

  645.8353052519245 A.wav vs C2.wav

I do not know how the Euclidean distance values look for similarity.
Is the similarity rate high when the distance value approaches zero?
How do you know there is a high match with this distance value?
Is not DTW a library to see how similar the voices are?

I will answer your questions in reverse order.

4: DTW (Dynamic Time Warping) is not a library but an algorithm. It allows aligning two sequences by warping them in time. You can use it for pretty much any kind of sequential data, for which a metric (distance) can be defined. Generally, you calculate a distance between each and every point of both sequences and produce a big matrix of accumulated distances. In the end, a path that produces the smallest final cost is selected as the optimal.

Here is an example of two 1D signals: $\color{blue} {x_1}$ is an exponential sweep sine and $\color{orange} {x_2}$ is a sinusoid with a constant frequency. The goal is to warp the $\color{orange} {x_2}$ in time so it would match the $\color{blue} {x_1}$ as close as possible. The figure below shows the outcome of such alignment. The signal $\color{orange} {x_2}$ had been stretched so that it matches $\color{blue} {x_1}$. However, $\color{orange} {x_2}$ is shorter than $\color{blue} {x_1}$, hence it had to be padded towards the end. The subplot at the bottom shows the evolution of the distance for each frame, where the final value is the accumulated cost - the very thing that you are using in your comparison. From this example, you can immediately see some limitations of the DTW.

Just to wrap up, here is the full DTW matrix together with the optimal path.

3 & 2: It's a distance, so the lower the better. The lowest possible distance is $0$, which can be achieved when:

Both signals are exactly the same
One signal is the stretched version of the other one, such as [0, -1, 0, 1, 0] and [0, -1, 0, 0, 1, 1, 0].

1: As I've mentioned before, DTW compares the distance between each and every point and then accumulates the cost of all operations that lead up to that point - the path that will yield the minimal cost corresponds to the best alignment.

In terms of similarity, signals will be considered close to each other if stretching or squeezing certain parts will minimize the difference between the MFCC vectors for each frame.

These answers might be relevant for you: click, click, click.


Phone conversation diarization with LIUM I am trying split call-center recording by speakers. LIUM seems to fit great. But results look random to me (green is etalon, gray is what I get):

So I think I am doing something wrong. Maybe features are not extracted correctly? Command I am running is

java -cp lium_spkdiarization-8.4.1.jar fr.lium.spkDiarization.system.Telephone \
     --saveAllStep --doCEClustering \
     --fInputMask sound.wav --sOutputMask=%s.xml \
     --sOutputFormat=seg.xml,UTF8 sound

sound.wav is RIFF (little-endian) data, WAVE audio, Microsoft PCM, 16 bit, mono 48000 Hz. One can download it here https://yadi.sk/d/w0doqi7QhcjJj . Please, help!

system.Telephone is just for detection of speech, that's why you get only init clusters.  You can use default fr.lium.spkDiarization.system.Diarization it will do proper steps as described in LIUM docs, result will be like this:

;; cluster S0
file 1 0 1135 M S U S0
file 1 5722 300 M S U S0
file 1 9266 359 M S U S0
file 1 12670 370 M S U S0
;; cluster S1
file 1 1135 832 F S U S1
file 1 5372 350 F S U S1
file 1 6022 572 F S U S1
file 1 10167 1550 F S U S1
file 1 13040 466 F S U S1
file 1 13623 352 F S U S1
;; cluster S10 [ merge HCLR 0 = S10 in S11 with 1.208107727685253 ]
file 1 1967 693 F S U S10
file 1 2660 1644 F S U S10
file 1 4304 1068 F S U S10
file 1 6594 937 F S U S10
;; cluster S12
file 1 7531 429 F S U S12
file 1 8692 574 F S U S12
;; cluster S13
file 1 7960 732 M T U S13
;; cluster S16
file 1 9625 542 F S U S16
file 1 11717 953 F S U S16

LIUM models included into jar are trained for 16khz, so the input should be better 16khz, you need to resample your 48khz file to 16khz before you feed into LIUM.

If you want to run on telephone speech, you might need to retrain mail/female/ubm models for that. It is also worth to note that LIUM models are trained for French, for other languages they should be better retrained.

Overall, LIUM is not expected to work out of box, you might need to spend some time to make it work properly.


Correct way to implement Discrete and Stationary Wavelet Transform? Can someone please confirm if this is indeed the correct way to compute DWT:

load image
load Haar wavelet

for i = 1: Level
  periodize edges of image
  convolve image with Low pass Haar wavelet filter and store result in G
  convolve image with Low pass Haar wavelet filter and store result in H

  convolve G with low pass Haar to calculate approximation co-eff
  convolve G with high pass Haar to calculate horizontal co-eff

  convolve H with low pass Haar to calculate vertical co-eff
  convolve H with high pass Haar to calculate diagonal co-eff

  i = i + 1
  downsample approximation by 2
end

In case of SWT, instead of down-sampling the approximation, Haar filter co-efficients are up-sampled twice at each level and the same process is repeated. Am I correct here?

In the Haar wavelet case, I just used the following NumPy:

# assume x is passed in as a variable
haar_transform = zeros(N)
i = arange(N/2, dtype=int)
haar_tranform[i]      = x[2*i] + x[2*i + 1]
haar_transform[i+N/2] = x[2*i] - x[2*i + 1]

I believe the mathematical description includes some scaling by $1/\sqrt{2}$.

I've learned that while academic papers may describe something very complicated there's a simple way to compute it.


How to prove that the peak of the autocorrelation function is at zero lag? Show that for a signal $f(\tau)$ with finite energy and energy autocorrelation function $\phi^e_{ff} (\tau),$$$|\phi_{ff}^e (\tau)| \leq \phi_{ff}^e (0), \ \ \forall \tau.$$

According to my textbook the total energy is given by $\phi_{ff}^e (0),$ but it does not provide any proof. The two energy auto-correlation functions here would become:

$$\phi_{ff}^e (\tau) = \int^\infty_{- \infty} f^*(t) f(t+ \tau) \ dt,$$

$$\phi_{ff}^e (0) =\int^\infty_{- \infty} f^*(t) f(t) \ dt.$$

I was told that we must use the Cauchy-Schwarz inequality for this proof. But I am not exactly where to apply it. Any help would be appreciated.

The Cauchy Schwarz inequality states that:
$$
\left|\int_{-\infty}^{\infty}g_1(t)g_2(t) dt\right|^2 \leq \int_{-\infty}^{\infty}|g_1(t)|^2 dt \int_{-\infty}^{\infty}|g_2(t)|^2 dt
$$

I'm going to assume that $f(t)$ is real, just to make the math a little easier. From the above we can write:
$$
\left|\int_{-\infty}^{\infty}f(t)f(t-\tau) dt\right|^2 \leq \int_{-\infty}^{\infty}|f(t)|^2 dt \int_{-\infty}^{\infty}|f(t-\tau)|^2 dt
$$
For the second integral on the right, just use a variable substitution - let $u=t-\tau$ and then rewrite it with the dummy variable $t$ again - we have
$$
\left|\int_{-\infty}^{\infty}f(t)f(t-\tau) dt\right|^2 \leq \int_{-\infty}^{\infty}|f(t)|^2 dt \int_{-\infty}^{\infty}|f(t)|^2 dt
$$
but the integral on the right is just the autocorrelation function i.e. $ \int_{-\infty}^{\infty}|f(t)|^2 dt= R_{f}(0)$. So we can write
$$
\left|\int_{-\infty}^{\infty}f(t)f(t-\tau) dt\right|^2 \leq R^2_{f}(0)
$$
Note that from the definition of $R_f(0)$ it must be a positive quantity, so taking the square root of both sides and paying attention to the positive and negative quantities on the left hand side, we can write
$$
\int_{-\infty}^{\infty}f(t)f(t-\tau) dt \leq R_{f}(0)
$$

Note - you don't have to use the Cauchy Schwarz inequality - another proof:
$$[f(t)-f(t-\tau)]^2 = f^2(t)+f^2(t-\tau) -2f(t)f(t-\tau)
$$
Integrating both sides from $-\infty$ to $\infty$ and using the same trick (variable substitution for the integral over $f^2(t-\tau)$) we can write
$$
\int_{-\infty}^{\infty} [f(t)-f(t-\tau)]^2dt= 2\int_{-\infty}^{\infty}f^2(t)dt-2\int_{-\infty}^{\infty}f(t)f(t-\tau)dt
$$

Rearranging we can write:
$$
\int_{-\infty}^{\infty}f(t)f(t-\tau)dt = \int_{-\infty}^{\infty}f^2(t)dt -\frac{1}{2}\int_{-\infty}^{\infty} [f(t)-f(t-\tau)]^2dt
$$
Now - the second integral on the right has to be a positive quantity and it is subtracted so we have
$$
\int_{-\infty}^{\infty}f(t)f(t-\tau)dt \leq \int_{-\infty}^{\infty}f^2(t)dt
$$


Bit Error Probability for 8PSK I'm trying to understand how to derive the BEP for 8PSK by geometrically looking at the decision regions.

Say I have this 8PSK signal set and the corresponding decision boundaries.

I assume that '000' is sent, and I want to calculate the BEP. So I have

$$P_b(e) = \frac{1}{3}(1 \times R_1 + 1 \times R_2 + 2 \times R_3 + 1 \times R_4 + 2 \times R_5 + 2 \times R_6 + 3 \times R_7)$$

Since $R_1 = R_4$, $R_3 = R_5$ and $R_2 = R_7$ by symmetry, the above can be simplified to

$$P_b(e) = \frac{2}{3}(R_1 + 2R_2 + 2R_3 + R_6)$$

Then we have

$$R_1 + R_2 + R_3 + R_6 = P\Big(n > d_1\sin\frac{\pi}{8}\Big)\tag{1}$$

and

$$R_2 + R_3 = P\Big(n_0 > d_1\sin\frac{3\pi}{8}\Big) \times P\Big(n_1 < d_1\sin\frac{\pi}{8}\Big)\tag{2}$$

I'm having a hard time visualizing and understanding how $(1)$ and $(2)$ are obtained.

I understand that $n$ is the noise vector, and $n_0$ and $n_1$ are the horizontal and vertical components of the noise vector, and we are trying to find out how far $n$ or $n_0$ and $n_1$ will go to "push" the signal into the wrong decision region.

I've seen the following from this website, and it seems to explain $(1)$, that is, the noise vector $n$ "pushes" the signal past the decision line (dl) into the wrong decision regions $R_1 + R_2 + R_3 + R_6$.

But what about $(2)$? I don't know why we have $\sin\frac{3\pi}{8}$ and $\sin\frac\pi8$.

Also, is the noise vector always perpendicular to the decision boundaries?

The way I learned this is there's two ways to find the probability of error. You can find the exact or as estimated probability. Exact is typically calculated for signal space where the points are in a rectangular pattern, because you're decision regions are nice and easy. You would see this rectangular pattern in a 4-QAM, 16-QAM, etc.

For estimated probability, we use the Union Bound. This states that the symbol probability of error = $P_{se|1} = \sum P_{se}($between point 1 and each other individual point$)$. Out total symbol probability of error

$P_{se} = $$\frac{\sum_{i=1}^M P_{es|i} }{M} $ where M is the number of points on our signal space ($M= 8$ in this case)

This means we're take a single point, and then with each other individual point we're going to create a new decision boundaries exactly half way between to the two and find the probability of error for just those two points. We do that for all the points, add them up, and now we have an estimated probability of error for the symbol. So let's say our noise is Gaussian, the variance of our noise $\sigma^2 = \frac{N_0}{2}$. Let's say D is the distance between two points. The probability of error between those two points $= Q(\frac{D}{2\sigma}) = Q(\frac{\sqrt2 \cdot D}{\sqrt N_0})$ So, let's look at point 000 in your examples. We can use cosine theory to find the distance between two points.

$D^2 = D_1^2 + D_2^2 -2\cdot D_1 \cdot D_2 \cdot cos(\theta)$  where $\theta$ is the angle between two points. Note $D_1$ and $D_2$ are the distances from the origin of our points, and in this case $D_1 = D_2$ because all points are
 the same distance from the origin.

$P_{se} =  Q(\frac{d_1\cdot \sqrt {2+\sqrt2}}{2 \sigma})  +  Q(\frac{d_1\cdot \sqrt {2+\sqrt2}}{2 \sigma})  +Q(\frac{\sqrt2\cdot d_1}{2 \sigma}) + Q(\frac{\sqrt2\cdot d_1}{2 \sigma}) +...$

$P_{se} = 2\cdot Q(\frac{d_1\cdot \sqrt {2+\sqrt2}}{2 \sigma}) + 2 \cdot Q(\frac{\sqrt2\cdot d_1}{2 \sigma}) +...$

We could go on, but since the $Q$ function falls off so quickly, we really only care about the point(s) closes to what we're comparing too. All the other terms are insignificant cared to that. So we end up with $P_se = 2\cdot Q(\frac{d_1\cdot \sqrt {2+\sqrt2}}{2 \sigma})$.

That's the symbol probability of error if we look at the point 000. We can calculate this for the other points, but it's just going to be the exact same due to symmetry. Since the Probability of error is the same for each symbol, we know the probability of error for all symbols is what we found above.

This is for symbol probability of error. To get to bit probability of error we have to look at our bit mapping. Since each point will only have two points that are closest to it, we can implement Gray coding. That means we can divide our symbol probability by $log_2(M)$ where $M =$ the number of different points we have on our signal space, which is 8 in this case. Therefore our $P_{be}$ is

$P_{be} = \frac{P_{es}}{log_2(8)} = \frac{2\cdot Q(\frac{d_1\cdot sqrt(2+sqrt(2))}{(2 \sigma)})}{3}$