Simple and efficient algorithm to detect frequency and phase of a sine signal I need an algorithm to detect frequency and phase of a pure sine signal. The frequency of the input signal changes between 0 and 100 Hz.

The value of the signal gets captured with a frequency of 20 kHz (so I get 20000 values per second) - this is given and cannot be changed. I need to detect frequency and phase of this input signal and using PWM generate MCU interrupts with the same frequency as the input signal.

Can anyone suggest what algorithm to use to do this simple and efficient? Maybe Goertzel algorithm?

Note: I originally posted this answer for the Stack Overflow copy of this question, before realizing that it had also been asked here.  It somewhat duplicates pichenettes' answer, but I felt it still worth (re)posting here, since it includes some extra details.  (Whether those details are useful or not, I'll leave for you and the OP to judge.)

If you know your signal is a pure sine wave, you can just use zero crossing detection.  Each cycle of the sine wave will have two zero crossings: one from negative to positive, and one from positive to negative.  This will be a lot simpler and more efficient than trying to do something fancy like Fourier transforms.

There are a few details to keep in mind, though:

It's OK for the signal to be slightly biased, but if the bias might exceed the amplitude of the signal, you'll need to correct it somehow.

You can either do this before sampling with an analog high-pass filter, or you can track a moving average of your sampled signal and use it as the "zero level" to compare to.  Or you can, instead, look at zero crossings of the difference between successive samples (which correspond to the maximum and minimum of each cycle in the signal) to avoid any bias issues.
If your input (or ADC) is noisy, your samples might randomly fluctuate around the zero level when the signal is close to it.  In such cases, naïvely comparing successive samples might detect multiple zero crossings where there's really just one.

One way to fix this issue is to smooth your signal before processing it, but it may be easier and more efficient to implement hysteresis in your zero-crossing detector.  That is, you'd only detect a positive-to-negative crossing when the signal dips below some pre-set threshold level −ε, and a negative-to-positive crossing only when it rises above +ε, like this pseudo-C code:

boolean isPositive = (firstSample > 0);
while (running) {
    int signalLevel = getNextSample();
    if (isPositive && signalLevel < -threshold) {
        isPositive = false;
        handleZeroCrossing();
    }
    else if (!isPositive && signalLevel > +threshold) {
        isPositive = true;
        handleZeroCrossing();
    }
}

In fact, you may not even need an algorithm for this at all, since, as noted in the comments below, you can implement zero-crossing detection with hysteresis in hardware with a simple Schmitt trigger, which basically converts the sine wave input into a square wave signal with the same frequency and (almost) the same phase, which you can then read as a simple digital input.  You might even be able to use the output of the Schmitt trigger to drive the MCU interrupt pin directly.


Restricting elemets in mfcc of a signal I am getting 20-length mfcc of an audio signal. Each element(vector) is of size 5. I am using Librosa.
How to give custom size to the each element of mfcc (say 6)?

If you are looking to get more mfcc coefficients for the given signal, you can do that by using librosa.feature.mfcc(y=y, sr=sr, n_mfcc=40), here you will get 40 coefficients for each time component.


Re-mixing high and low pass filtered audio has interference in the shared transition band? Disclaimer: this feels like something that should be obvious from an intro class, but for some reason I can't find anything around this.

I'm working on an algorithm that takes a signal, let's call it S_in. I'm running both a high pass and low pass filter (identical parameters otherwise) on it so that I can work on different frequency bands of the signal with slightly different parameters in my algorithm. So, lpf(S_in) and hpf(S_in). When I mix them back together (which I'm doing in the time domain via just simple addition, S_out = lpf(S_in) + hpf(S_in)), I get significant "interference" in the transition band in the resulting signal, usually in the form of something that resembles a band stop/notch filter at my original cutoff frequency. This happens even if I don't run those signals through my algorithm (i.e., S_in != S_out; preferably they would be approximately equal to each other, instead of significantly different).

As a "workaround" I noticed if I offset the addition by as few as 8 or 16 samples on one of the signals (that is, a time delay on one of the streams; also 44100 samples/sec audio fwiw), then the frequency response of the mixed signal looks much more flat relative to my original signal (which is what I want -- the band stop is gone, though there's still clearly some interference pattern in the shared transition band).

What am I doing wrong here? Is this a simple signal processing thing I've forgotten about and can't figure out how to search for it? Is it an artifact of discretization/rounding errors of my implementation? (Though since I can reproduce that same behavior in common audio editing software, I don't think it's just my implementation). I'm using a low order biquad filter for the filters (is it that? obviously a higher order filter would shrink the transition band, but is that the "right" fix, or just a workaround?). Is there a "right" way to mix them back together, or a "right" way to do my overall algorithm where I need to process frequency bands separately and then mix them back together?

You are assuming that a parallel high pass and low pass add up to "flat". In most cases, they do not. You have to carefully design both the phase and the amplitude responses in the transition band to match each other.

There are a few classes of filters that are useful here

1st order lowpass/highpass to add up to unity
Odd order butterworth filters add up to an allpass filter (flat frequency response, but some phase distortion). The two paths are either 90 or 270 degrees phase shifted (depending on order)
Linkwitz Riley filter (two butterworth cascaded) also add up to an allpass (when choosing the appropriate sign for recombination). The two paths are in phase or out of phase (depending on order).
linear phase FIR filters that are designed complementary.
Specific types of filter banks (reconstructing)

All of this only helps if your processing doesn't induce significant differential phase shifts. Unless this phase shift is constant and can be compensated, you will always get some amount of interference.


How to implement bandpass filter on complex valued signal? I am using Scipy to implement bandpass filter but it assumes that positive normalized frequency is passed but I & Q samples range from [Fc-Fs/2,Fc+Fs/2] where Fs is sampling frequency & Fc is centre frequency but desired bandpass filter should filter from [f1,f2] such that they can negative too. How to pass this band information ?

My approach is

Generated a window using Scipy firwin
Convoluted it with signal using Scipy fftconvolve

You can do different things: For example,

use the (frequency-symmetric) real-tapped bandpass that firwin gives you, and after applying that, apply a complex high pass (a Hilbert filter, essentially) to kill all negative frequencies. That sadly leaves you with the original problem (finding a complex-tapped filter using scipy)
Do the same as above, but already convolve the real-valued band pass filter taps with the Hilbert filter taps, so that you only need to apply one filter
Do the "usual" lowpass-to-bandpass transform trick: You design a real-valued low pass filter (e.g. using firwin) with impulse response (==taps) $h_{LP}[n]$ that has the same passband and transition widths as your desired complex band pass. Afterwards, you multiply that $h_{LP}$ with the shift you need to move the 0-centric low pass up to your $f_\text{center}$-centric band pass: $$h_{BP}[n] = h_{LP} \cdot e^{j2\pi \frac{f_\text{center}}{f_\text{sample}}n}\text.$$ The idea is that you shift the low pass in frequency by convolving with a dirac impulse in frequency domain, which is equivalent to multiplying with the complex sinusoid of said frequency in time domain. These answers (1,2) might be of interest to you.


Calculate algorithm's MIPS from cycles? My audio DSP algorithm processes input sound file (45 seconds, 16 kHz, mono). After simulation is done, gprof profiler tells me:

Cycles: total = 6846310893.

My target DSP processor is 600 MHz and executes 1 instruction per cycle, meaning it is 600 "MIPS" (millions of instructions per second).

How much will this algorithm (runnig real-time) take out of those 600 MIPS?
Can this at all be calculated using above data?

The CPU's datasheet will provide information about how many cycles each command requires (e.g. in here look for "cycles", and with something like this you can also work out what else might a command require to determine if additional cycles will be required). Sometimes they might differ (for the same command), depending on the operands.

What you can do is disassemble the processing part, look at the set (i.e. unique) of commands that are used, look up how many cycles they take on average and from that, work-out a best-case / worst-case scenario for your algorithm (taking also into account how many times each command is used, potentially even, in a loop).

If you are already at a low-enough level, so that you are sure about what commands have been used, then your estimate will be closer to the real value.

Hope this helps.


On the time-frequency relationships of median filters I'm applying Matlab's median filter to a data sequence and I'm trying to figure out what filter co-efficients of a linear filter this corresponds to.

Would the normal time-frequency relationship that exists for linear filters maintain itself for the median filter?

The median filter is actually an example of a non-linear filtering operation.

This stands in contrast to linear filtering operations (the 'classical' way that involve the convolution of a filter's co-efficients with a signal). Since convolution in time is the same as multiplication in the frequency domain, you can go back and forth between the domains to figure out co-efficients of your filter in either basis.

However, since medial filtering is a non-linear operation, you cannot do this. There is nothing you can convolve your signal with (ie, no linear operation you can do on your signal) such that you would get the same result as you see by using the non-linear medial filter.


What is the frequency of $\cos x -\sqrt 2\cos\sqrt 2 x$? Question

Considering that $1\over{2\pi}$ is the frequency of $\cos x$ and also of $\cos x - 2\cos 2x$, what is the frequency of $\cos x - \sqrt 2\cos\sqrt 2 x$?

Thoughts

Perhaps "frequency" isn't the interesting quantity. I think the period is infinite. That would make the frequency zero. Maybe what I want to know is bandwidth?

I think the additive combination (superposition) of the two waves produces a component at the higher frequency $1/(2\pi(\sqrt 2-1))$. I would hope that sampling at a frequency of $1/(\sqrt 2\pi -\pi)$ would be sufficient to reveal as many zero-crossings as we may wish. Is that correct? Is it the best answer?

I understood Fourier series and partial differential equations when I took that course and I still get the general idea, but it has been a long time.

Motivation

Mistakenly, I devised and posed this question for my colleagues as a weekly math challenge before I realized that I don't know how to get the answer myself.

I'll provide additionally just some background about my motive for asking this particular question. Maybe it can help you give an answer that is relevant to me. However, this background information in itself is not part of my question. So, I came up with the question about frequency hoping to move us nearer to another question that Mark Kac used in order to elucidate stochastic independence. That is, given any particular value $-2 <\alpha < 2$, in what fraction of the interval $0\leq x < 100$ (just to choose some arbitrary bounds) is $f(x) <\alpha$ where $f(x) = \sin x - \sin\sqrt 2 x$? The interesting destination that I want to reach is that the answer to this further question is (somehow related to) the normal distribution function, yet there is nothing "random" involved.

Of course $f'(x) =\cos x -\sqrt 2\cos\sqrt 2 x$. If I know the "frequency" (or something) of $f'(x)$ then I hope I can find, at least numerically, its zero-crossings. Those include all of the local extrema of $f(x)$. Then, knowing that the sum-to-product formula gives $f(x) = 2\cos(x/2+x/\sqrt 2)\sin(x/2-x/\sqrt 2)$, I think that my favorite computer algebra system and I can find, at least numerically, where $f(x) <\alpha$.

Not sure, whether it's really your question. The signal for which you asked for the frequency actually has two frequencies. There are two components, the two cosines. And each component has a frequency. Assume your frequency is called $f$. So, the signal looks like $\cos(2\pi f x)$. So, getting the frequency of a component is quite easy. Just take everything from it between the brackets and divide it by $2\pi$.

Probably, you already got the frequencies of your signal. But just to be sure, it's: $\frac{1}{2\pi}$ and $\frac{\sqrt{2}}{2\pi}$.

Hope, it was the answer, you're looking for...


$ {L}_{0} $ Pseudo Norm Minimization in Compressive Sensing I have recently taken up studying compressive sensing related papers. Some things are not very clear to me or may be I am not able to visualize the scenario as is said. Like how $l_0$ norm minimization is NP hard? Can somebody explain with an example (step by step example may be?) of a matrix and vector for an equation $y=Ax$, how $l_0$ minimization of the same can be NP hard?

the solution for a sparse recovery problem is given by:
$$\text{min} ||x||_0$$
                    $$\text{s.t} \hspace{2mm}  y = Ax$$

The definition of $||x||_0$ is no. of non-zero entries in $x$. This is also called the sparsity of the vector.

i.e., we are asking for the sparsest solution $x$, that satisfies $y = Ax$.

Consider the simplest case where $x$ is 1-sparse but you don't know the location of that non-zero entry. In such a case, we have $\binom{N} {1}$ possibilities for a 1-sparse vector and to find the solution we have to examine the values of all the $\binom{N} {1}$ possibilities for finding the unique minimizer.
Similarly, if you are told that $x$ is $k$-sparse, you need to search $\binom{N} {k}$ possibilities.
i.e., the algorithm grows as $\binom{N}{k}$ with increase in $k$. Since $k$ is not known a priori, you have to check for all the $N$ possible values of $k$.

Hence, the complexity of the algorithm is $\displaystyle\sum_{i=1}^N \binom{N}{i}$. This is called as NP hard problem which is an acronym for Non-polynomial time complexity.


Software to extract frequency spectrum data from an audio file I've been trying to learn how computers recognize musical notes, and I found apps like spectrogram which shows the frequency spectrum of a sound, but what I need is the data which is plotted or the raw data used. I'm a beginner in this field, so it'd be really great if you could help me with the software or the process used. Thanks!

An important thing to consider is that the pitch of musical notes does not directly correspond to frequency peaks in a spectrogram.  This is because the energy of typical musical sounds are partially to mostly in the overtones and harmonics, at frequencies well above the note pitch (depending on the source of the music).  So it would help to investigate algorithms for pitch detection/estimation, and the difference between those and the frequency spectrum of music.


Looking for a way to monitor raw microphone digital output? This is sort of a cross between a data processing and Physics question, but I would be grateful if anyone here could help me.

My objective is to calculate the sound pressure level (in dB) of a sound generated 1 meter away. I would like to use my condenser microphone to calculate this. The microphone I have is the Blue Yeti, which has a sensitivity of 4.5mV/Pa (reference value 20 microPascal).

So my question is is there a way that the digital output of the microphone can be directly monitored and recorded in units of Volts. Hence, I would use these values to calculate Pascals. Can anyone point me to a (preferably free) program that I can use to do this direct monitoring?

So my question is is there a way that the digital output of the microphone can be directly monitored and recorded in units of Volts. Hence, I would use these values to calculate Pascals. Can anyone point me to a (preferably free) program that I can use to do this direct monitoring?

So, in fact, your sound card is a Analog-to-Digital Converter, which does exactly that: convert the voltage coming from a microphone into a stream of numbers. Congratulations, you already have your data acquisition output :)

That stream of numbers has no physical unit (it is "best-case proportional to input Voltage", but it is a number, not a physical entity anymore)– you will need to calibrate it with a couple of known signals.

To make matters more complicated, the amplitude response of a soundcard can't, technically, be flat. Due to the [Shannon-Nyquist theorem], the microphone voltage signal must first go through a low-pass filter that limits the maximum frequency contained in the signal (you're a physicist: apply Fourier transform to time signal $s(t)$, giving you $S(f)$, find $f_{max}$ for which $S(f>f_{max}) \approx 0$).

Since perfect analog filters don't exist, the suppression "edge" can't be infinitely steep, so you'll see changing effect of voltage on digital numbers with changing frequency.

If you work with only one (or few) frequency(ies), and you have a multiple reference pressure wave source for this(those) at different powers, then go ahead and use the sound card, calibrate using the known powers, and calculate correction factors.

If you need something that actually promises linearity, sound cards with their analog filters optimized for audio perception might not be the tool of choice. Verify your approach with an oscilloscope, then use a suitable data acquisition unit – some digital oscilloscope can record signals for you, there's even digital oscilloscopes that have no display and are only meant for usage with a PC, other people I know use Software Defined Radio (SDR) devices to sample audio/sonar frequency signals instead of RF, yet others have just built their own microcontroller-based data acquisition devices based on 10€ evaluation boards.


Understanding 2D FFT result of an image having a pattern I have created Fast Fourier Transform (FFT 2D) of the following image (without pattern angle markers) using ejectamenta online tool.

Input (Original image is here)

Output I got is the following

As well as input pattern 66 deg gives output with 48 deg maximum magnitude line. (angles are measured approximately). But I believe Fourier pattern should be perpendicular to the input image pattern. But it has a difference. How this happens and what is the problem here. Could someone please explain? I'm trying to derive pattern direction using it's FFT result. Thank you.

I found that FFT output depends on aspect ratio of the input image. Fourier pattern is perpendicular to the input image pattern only if the input image has same width and height.


Sine Wave Aliasing during IFFT I’’m synthesising a simple sine wave using the IFFT.   I’m controlling the pitch of the sine wave by controlling the bin number and the phase for that bin number.

I’ve synthesised an F whose frequency is very close to a bin centre frequency and everything looks good.

however when I detune that same note away from the bin centre frequency I start to get additional harmonics showing up.   Like this:  (Note the -30db peak just after the significant peak, and the appearance of signal all the way to the left.)

Why is this occurring???

This is likely due to a truncation of your sine wave in the time interval you have chosen, such that there is not an integer multiple number of sine wave cycles chosen. I am guessing that your frequency view that you provide above is a subsequent measurement (DFT) from the resulting time domain waveform you have created, and not a view of your original frequency function you are using to create the time domain waveform (as the original function in the frequency domain would not have those artifacts!).

That said, the specific answer will depend on how exactly you create and repeat/extend the final time domain waveform from the fixed length IFFT result that you get, but below are some insights on the effects of waveform truncation in the time domain that may help you to understand how you may be getting such artifacts.

A direct and concise mathematical answer of frequency effects from waveform truncation can be derived by describing your fixed length (sampled) waveform in the time domain as an infinite length waveform multiplied that is not truncated by a rectangular (boxcar) window. Multiplication in the time domain is convolution in the frequency domain, so the result of a boxcar window truncation would be convolving the frequency you would get without truncation with a sinc function (or when sampled, the approximation of a sinc function).

I like to view the resulting effect by picturing the time domain waveform extending from minus infinity to positive infinity in time, made up of your specific waveform snippet (IFFT result) continuously repeating in time. For a direct DFT this analogy this is mathematically valid (meaning the DFT over the finite time interval would match the DTFT of the same sequence repeating it time (just that the DTFT would be a continuous function and the DFT would be samples of this function- see below). Repetition in the time domain over an interval T results in a FT with discrete impulses in the frequency domain spaced at 1/T.).

From this view it may be clearer to see the effects of truncation: if we had 1.5 cycles of a sine wave in your time interval, the sine wave would suddenly snap every 1.5 cycles in the repeating waveform, requiring higher frequencies to achieve this "snap" that takes place in time. In contrast if you had an exact integer number of cycles (for a samples system, samples 0 to N-1, where sample 0 would be identical to the sample N for a pure undistorted sine wave), then the FT result (and DFT and DTFT) would only exist at one frequency location with 0 elsewhere; as expected for a pure tone.

FT, DTFT, DFT Background

I reference FT, DTFT, and DFT above, so I included this graphic and primer below to help clear up any confusion for those less familiar.

What is shown is the FT (or CTFT), DTFT, and DFT for a pure tone (that in "truth" extends from minus infinity to infinity in time), along with a modulated waveform of finite bandwidth.  In the Fourier Transforms, the tone is shown in blue while the modulated waveform is shown in red, and $\omega_s$ is the sampling frequency.

The CTFT is the Continuous Time Fourier Transform, or Fourier Transform (FT) for short.

$$ X(\omega) = \int_{t=-\infty}^\infty x(t)e^{-j\omega t}dt $$

Observe:

continuous in time -> aperiodic in frequency

aperiodic in time -> continuous in frequency

The DTFT is the Discrete Time Fourier Transform; the time interval still goes from minus infinity to infinity, but the system is now sampled. It is a sampled FT, which create repetition in frequency.

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

Observe:

discrete in time -> periodic in frequency

aperiodic in time -> continuous in frequency

The DFT is the Discrete Fourier Transform; the time interval is defined over a finite interval. There is an implied periodicity in time, meaning if you made a waveform by repeating this one in time from minus infinity to infinity, the DTFT of this periodic waveform would match the DFT at the DFT sample locations (see plot above).

$$ X[k] = \sum_{n=0}^{N-1}x[n]  e^{-j\omega_o k/n} $$

Observe:

discrete in time -> periodic in frequency (implied)

periodic in time (implied) -> discrete in frequency

*note the time domain and frequency domain waveforms are not really periodic, as they only exist over a finite number of samples- but I like to make this interpretation, as if you did repeat the waveform in both domains you would get the same valued discrete samples (but they would be continuous waveforms with zero values in between). This interpretation has been extremely helpful for me in working with mixed signal (analog/digital) systems.


Difference between "Discrete Wavelet Transform" and "Discrete Wavelet Decomposition" I have a rough overview on Discrete Wavelet Transform (DWT). However, I am confused about Discrete Wavelet Decomposition and did not find a good reference yet which explains this well. What is it actually about? Is it somehow part of DWT or an inverse operation to it?

The discrete wavelet transform should denote "the operations" that, applied to some data, yield a discrete wavelet decomposition. The first one can be seen as a matrix operator, while the second relates to the actual wavelet coefficients, or the structure of thereof, that you would obtain after the application of the first one.

In everyday language, they are often used interchangeably, by a form of metonymy:

  a figure of speech in which a thing or concept is referred to by the
  name of something closely associated with that thing or concept


Extract parameters of a single harmonic Context

I need to extract the parameters of a single harmonic in a signal.

For instance if

$$
x(t) = A \sin (\omega_0t + d_0) + B \sin (\omega_2t + d_2)
$$
then $f(x, \omega_0) = (A, d_0)$. This can be done offline.

Please consider the fact that $t$ does not cover a full number of periods!

Of course I used a fft() and then extracted the coefficient corresponding to the interesting frequency. The problem is that I compute the whole spectrum when I'm interested in only one value.

What I tried and read

I first tried to use the definition of the DFT (or to project my signal on a sine/cosine plane) by using:

$$
A e^{jd_0} = \frac{2}{N} \sum_t e^{j\omega_0 t} x(t)  ,
$$
with $N$ the number of samples of $t$.

This is what is proposed here:
FFT for a single frequency

This works however only for a full number of periods of the interesting frequency: I can cut my signal to get reasonable results but I'd rather not (and it's not always possible)

What can I do else?

I tried something that I didn't thought of before:

force $\omega_0$ to be a multiple of $2 \pi \frac{fs}{N}$, with $fs$ my sampling rate.

And of course, this is what happens in the fft.

I only ask myself if I can get the same results for $\omega_0\neq2 \pi \frac{k}{N} fs$


Is It Possible to Reconstruct the Original Image from Laplacian Pyramid? In our computer vision course, my professor said we can reconstruct the original image from the Laplacian pyramid. It just feels to me that information is lost after doing subsampling and upsampling.

Is it possible? If so, why? Are the reconstructed image exactly the same as the original pixel by pixel? And is this answer correct?

Generate Laplacian pyramid $(f_2, h_1, h_0)$ from image $f_0$:

Reconstruct $f_0$ from Laplacian pyramid:

source

This depends on the order of upsampling and downsampling.
If the order is correct, then you won't throw away anything and thus you should in principle be able to reconstruct the image.

In general:

$$ \left(\uparrow_n\downarrow_n f\right) \neq \left(\downarrow_n\uparrow_n f\right) $$

Similar things are used when using the Wavelet decomposition on a signal, which again uses up- and downsampling combined with the Wavelet filters.


BER vs Eb/N0 (db) Curve So, I was trying to plot BER (Bit Error Rate) vs. Eb/No (in dB) curve using Python. Here is what I get,

Eb/N0 is in dB on the x-axis

But I know that this is not the correct graph.

Here is my code.

import matplotlib.pyplot as plt
import numpy as np
from scipy.special import erfc

def Q(x):
    return 0.5*erfc(x*(2**(-0.5)))

ebnodb = np.arange(0.1, 12, 0.01)
ebno = [10**(i*0.1) for i in ebnodb]
ebno = np.array(ebno)
y = Q(np.sqrt(2*ebno))

plt.plot(ebnodb, y)
plt.xlabel("Eb/N0")
plt.ylabel("BER")
plt.show()

Can someone figure out, what am I doing wrong?

Usually, such figures show the error rate on a log scale, and $E_b/N_0$ in dB. So you would specify $E_b/N_0$ in dB, then convert it to a ratio, and use that value to evaluate the formula for the bit error rate (as you did). Then you plot the logarithm of the error rate as a function of $E_b/N_0$ in dB. So the main problem is that you used a linear scale for the BER instead of a logarithmic scale.

If you do things right, you'll end up with something like this:


How to decode a Barker code Let's suppose that I send the following sequence of bits:

$x = [1, 0]$

I use a 7 length Barker code ($b = [−1,−1,−1,1,1,−1,1]$) in a BPSK modulation scheme, resulting in the following sequence:

$x_B =\\
[−1,−1,−1,\:\:1,\:\:1,−1,\:\:1\\
\:\:\:1,\:\:\:1,\:\:1,-1,-1,\:\:\:1,-1]$

We send it, through a AWGN channel, and after decoding there are some errors in the sequence:
$y_B =\\
[\:\:1,−1,\:\:1,\:\:1,\:\:1,−1,\:\:1\\
\:\:\:1,-1,\:\:1,-1,-1,\:\:\:1,-1]$

How would I decode this sequence?

Under the white Gaussian noise assumption, the optimal decoder would be to minimize each window of $7$ bits, that is take you first 7 symbols
$$s_1 = [\:\:1,−1,\:\:1,\:\:1,\:\:1,−1,\:\:1]$$
and compute two metrics
$$c_0 = \Vert s_1 - b_0 \Vert_2^2 = 2^2 + 2^2 + 2^2 + 2^2 + 2^2 = 20$$
$$c_1 = \Vert s_1 - b_1 \Vert_2^2  = 2^2 + 2^2 = 4$$
where $b_1 = [−1,−1,−1,1,1,−1,1]$ is the Barker code corresponding to $1$ and $b_0 = [1,1,1,-1,-1,1,-1]$ is the Barker code corresponding to $0$. So if $c_0 < c_1$ then the first received bits is $0$, else $1$.

Next, grab the other $7$ symbols
$$s_2 = [1,-1,1,-1,-1,1,-1]$$
and do the same thing.

PS: You could have also computed a Hamming distance as the criterion chooses the minimum number of errors.


Bessel filters: linear passband phase characteristic only for low pass filters? I've found a test (internet) from the contest to signal processing teacher and it has a question (true) or (false).

I cant find any book or paper to argue this statement. Is it really true?

Filter types other than low pass (high pass, band pass, band stop) are obtained from a low pass prototype filter by frequency transformations (such as $s\rightarrow 1/s$ for a low pass / high pass transformation). These transformation destroy the optimality of the phase response. Nevertheless, filters obtained from a Bessel low pass prototype filter usually exhibit a passband phase response that is closer to a linear response than the phase response of the other standard filter types (Butterworth, Chebyshev, Cauer).


How to locally quantify the 'sharpness' of an image? I am trying to quantify how much sharpness (or acutance) is in a picture which has some bokeh (out of focus background).

I am using the Python scikit image for that. Here is my naive approach:

import matplotlib.pyplot as plt

from skimage import data
from skimage.color import rgb2gray
from skimage.morphology import disk
from skimage.filters.rank import gradient

cat = data.chelsea() # cat is a 300-by-451 pixel RGB image
cat_gray = rgb2gray(cat)

selection_element = disk(5) # matrix of n pixels with a disk shape

cat_sharpness = gradient(cat_gray, selection_element)

plt.imshow(cat_sharpness, cmap="viridis")
plt.axis('off')
plt.colorbar()

plt.show()

So here is the picture (notice the background not in focus)

And here is the gradient, which actually measures the difference in contrast:

This difference in contrast problem is more obvious if one uses the Lena image:

Here, the background is caught as well due to the difference in contrast (there is a black frame on the top right)

Any ideas about how to give an scalar value where only the focused areas are highlighted?

The recent works I am aware of make use of tools that go beyond mere gradients. Here are a few references that could be starting points:

S3: A Spectral and Spatial Measure of Local Perceived Sharpness in Natural Images, 2012, with examples of sharpness maps and Matlab code (that could  be converted to Python)

  This paper presents an algorithm designed to measure the local perceived sharpness in an image. Our method utilizes both spectral and spatial properties
  of the image: For each block, we measure the slope of the magnitude
  spectrum and the total spatial variation. These measures are then
  adjusted to account for visual perception, and then, the adjusted
  measures are combined via a weighted geometric mean. The resulting
  measure, i.e., S3 (spectral and spatial sharpness), yields a perceived
  sharpness map in which greater values denote perceptually sharper
  regions

Image sharpness assessment based on local phase coherence, 2013, with examples  and (BROKEN yet) code (local copy of the code provided by the authors)

  Sharpness is an important determinant in visual assessment of image
  quality. The human visual system is able to effortlessly detect blur
  and evaluate sharpness of visual images, but the underlying mechanism
  is not fully understood. Existing blur/sharpness evaluation algorithms
  are mostly based on edge width, local gradient, or energy reduction of
  global/local high frequency content. Here we understand the subject
  from a different perspective, where sharpness is identified as strong
  local phase coherence (LPC) near distinctive image features evaluated
  in the complex wavelet transform domain. Previous LPC computation is
  restricted to be applied to complex coefficients spread in three
  consecutive dyadic scales in the scale-space. Here we propose a
  flexible framework that allows for LPC computation in arbitrary
  fractional scales.

The given examples and comparisons across different  could provide you with  some hints toward your goal.


Recover signal from its spectral amplitude and envelope I am given an amplitude spectrum (no phase) and the envelope of a signal. What I know about the original signal is, that it can be thought of as wavepacket like structure (sinusoids under a bell curve) and the spectrum is narrow banded.
To be able to further analyse the original signal I would like to recover the input signal from those two given functions.

My questions are: Can I recover the signal from its amplitude spectrum and envelope or would you advice me to do otherwise.

To clarify what is given:

The signal is thought of wavepacket structure, which is

$$x(t) \approx \sum_n\frac{1}{\sigma_n\sqrt{2\pi}}e^{-\dfrac{1}{2}\left(\dfrac{x-\mu_n}{\sigma_n}\right)^2} \cdot \sin(2\pi f_n x)$$

with $n \gt 1$, $\lvert\mu_k-\mu_l\rvert \lt \sigma_k, \sigma_l$ and if $\mu_k \gt \mu_l$ then $f_k \lt f_l$.

The amplitude spectrum is the absolute of the Fourier transform of the signal, which is

$$a(\omega) = \lvert\mathcal{F}\{x(t)\}\rvert$$

The envelope can be thought of as absolute value of the analytic signal, which is

$$e(t) = \lvert x(t) + j \cdot \mathcal{H}\{x(t)\}\rvert$$

I'm not sure if it is possible or even that beneficial to recover the signal and I am pleased for any answer or every hint.

One counterexample is two wave packets of identical spread and amplitude, well separated in time so that their envelopes do not overlap, and of high enough (differing) well-separated frequencies such that their analytic signals can be considered Gaussian-enveloped complex exponentials. The Hilbert envelope will thus consist of two identical well-separated Gaussians. The amplitude spectrum will consist of well-separated Gaussians of identical magnitude. There is not enough information to assign the frequencies to the wave packets.

So in general you cannot resolve the original signal, not even up to an ambiguity of a constant phase shift over all frequencies.

With sufficiently large time spreads of the two wave packets their frequencies need not be far apart to still be well separated in the spectrum, satisfying your constraint for a narrow-banded signal.


short time energy matlab code for audio feature extraction Is there code in matlab for short time energy feature extraction in matlab??

I want to extract audio feature (Short time energy) using matlab code

Short time energy (STE) is simple feature and so is the code to compute it.

ste = sum(buffer(signal.^2, winLen));

where  and winLen is duration of rectangular time window (in samples).

If you want to use diffrent shape of window, then you have to multiply each column returned by buffer by this window, before you sum the columns.


How to apply convolution on a buffer I'm trying to do some simple filtering for an audio signal using a window-sinc low pass filter. Supposing that my buffer has M values and the filter kernel size is N, after convolving these two arrays I would obtain an output of m+n-1 values. Which part should I take into consideration? I'm asking this because, if I read a M values array I should write a M values array. I have tried to take the first M or the last M values for output, but the result is not too good. I was wondering if I have to use a circular buffer in which I put every time M+N-1 values, but I read the first M values.

Output the first M samples of your convolution result
Keep the remaining samples and ADD those to the result next buffer

Google "overlap add" for more information. While overlap add is a frequency domain method, it explains the framing and buffer handling well.


What is going wrong with my IDFT program? I learn IDFT and have created two programs.

Input x(n). It calculates input x(n) values from 2 sinusoids and outputs them.
Output x(n). It performs IDFT on same signal and outputs results.

Results should agree. Unfortunately while first program's output is okay, second's is not. Can you please find the error and say me?

Input x(n) program

// Input x(n) program.
// Grid is created by variables lowTestFreq and highTestFreq
// both of double type, and integer type M.

#include <iostream>
#include "math.h"

int main()
{
//Unchanged variables.
float pi = 3.141592653589793; //pure number.
int n = 0; // sample.
float input_xn = 0;

//Sampling frequency and samples of segment.
int samplFreq = 100; // sample/second.
int N = 40; //sample.

// 1st existing sinusoidal.
float ampl_1 = 70000; //pure number.
float existFreq_1 = 36; // Hz.
float phase_1 = 0.67 * pi; // rad.
float unitAngle_1 = (existFreq_1 / samplFreq) * (2 * pi); // rad/sample.

// 2nd existing sinusoidal.
float ampl_2 = 600000; // pure number.
float existFreq_2 = 25; // Hz.
float phase_2 = -0.37 * pi; // rad.
float unitAngle_2 = (existFreq_2 / samplFreq) * (2 * pi); // rad/sample.

for (n = 1; n <= N; ++n)
{
input_xn =  ampl_1 * sin(n * unitAngle_1 + phase_1)
      + ampl_2 * sin(n * unitAngle_2 + phase_2);
std::cout << "input-x(n)" << input_xn << std::endl;
}
return 0;
}

Output x(n) (IDFT) program

// Grid is created by variables lowTestFreq and highTestFreq
// both of double type, and integer type M.

#include <iostream>
#include "math.h"

int main()
{
//Unchanged variables.
float pi = 3.141592653589793; //pure number.
int n = 0; // sample.
float input_xn = 0; //pure number.
float Re = 0; //pure number.
float Im = 0; //pure number.
float xn_cos = 0; //pure number.
float xn_sin = 0; //pure number.
float co = 0;
float si = 0;
float output_xn_Re = 0;
float sumReal = 0;

//Sampling frequency and samples of segment.
int samplFreq = 100; // sample/second.
int N = 40; //sample.

// 1st existing sinusoidal.
float ampl_1 = 70000; //pure number.
float existFreq_1 = 36; // Hz.
float phase_1 = 0.67 * pi; // rad.
float unitAngle_1 = (existFreq_1 / samplFreq) * (2 * pi); // rad/sample.

// 2nd existing sinusoidal.
float ampl_2 = 600000; // pure number.
float existFreq_2 = 25; // Hz.
float phase_2 = -0.37 * pi; // rad.
float unitAngle_2 = (existFreq_2 / samplFreq) * (2 * pi); // rad/sample.

// Testing grid.
float testFreq = 0; // Hz.
float testUnitAngle = 0; // rad/sample.
float testAngle = 0; // rad.
float lowTestFreq = 0; //Hz.
float highTestFreq = samplFreq; //Hz.
int M = N; // pure number.
int m = 0; // pure number.
float testFreqStep = (highTestFreq - lowTestFreq) / M; // Hz.

for (m = 0; m <= M-1; ++m)
{
testFreq = lowTestFreq + m * testFreqStep; // Hz.
testUnitAngle = (testFreq / samplFreq) * (2 * pi); // rad/sample.

Re = 0; Im = 0;

for (n = 1; n <= N; ++n)
{
input_xn =  ampl_1 * sin(n * unitAngle_1 + phase_1)
      + ampl_2 * sin(n * unitAngle_2 + phase_2);

testAngle = n * testUnitAngle;
xn_cos = input_xn *  cos(testAngle);
xn_sin = input_xn * -sin(testAngle);

Re += xn_cos;
Im += xn_sin;
}
sumReal = 0;
for (n = 0; n <= N-1; ++n)
{
co = cos(2 * pi * m * n / N);
si = sin(2 * pi * m * n / N);
output_xn_Re = Re * co - Im * si;
sumReal += output_xn_Re;
}
std::cout << "output-x(n): " << sumReal / N << std::endl;
}
return 0;
}

I suspect reason is that my IDFT program takes into account every leakage frequency.

Also in complex numbers multiplication I do not take into account imaginary unit.

I apologize for the errors.

My answer to my question is the following C++ program. It performs DFT and IDFT on 20 input values in array, and outputs the same values. Sample rate for calculate input signal from existing sinusoids is 71114 Hz, and N = 20. Existing sinusoids are two:

1. Amplitude = 995,574, frequency = 20247 Hz, phase = 0.34 π.
2.    "      = 492,356      "     = 3818 Hz,    "   = 0.15 π.

Sample rate for DFT and IDFT is 45940 Hz, and N = 20. I point out that existing frequencies are not multiples of sample rate by N. So, there is leakage but is not problem. Any one can change values. These in program are random in limits. Any one can use first program in my question, for calculate input x(n) from any sinusoids.

     #include <iostream>
     #include "math.h"

     int main()
     {
     //Unchanged variables.
     float pi = 3.141592653589793; //pure number.
     int n = 0; // sample.
     float Re = 0; //pure number.
     float Im = 0; //pure number.
     float xn_cos = 0; //pure number.
     float xn_sin = 0; //pure number.
     float cosine = 0;
     float sine = 0;
     float sumReal = 0;
     float real = 0;
    int m = 0;
    int no = 0;

    //Sampling frequency and samples of segment.
    int samplFreq = 45940; // sample/second.
    int N = 20; //sample. Number of input samples.

   // Testing grid.
    float testFreq = 0; // Hz.
    float testUnitAngle = 0; // rad/sample.
    float testAngle = 0; // rad.
    float lowTestFreq = 0; //Hz.
    float highTestFreq = samplFreq; //Hz.
    int M = N; // pure number. Number of DFT output samples.
    float testFreqStep = (highTestFreq - lowTestFreq) / M; // Hz.
    float xn[] =
{ 1.09595e+06,
 635609,
-544799,
 640870,
 1.40537e+06,
-142213,
-392826,
 1.00018e+06,
 307386,
-1.16555e+06,
-211430,
 518555,
-998272,
-1.20759e+06,
 390434,
 26244.3,
-1.18203e+06,
 25149,
 1.08298e+06,
-200068,};
   int No = N; // pure number. Number of IDFT output samples.

  // IDFT including DFT.
  for (no = 0; no <= No - 1; no++)
    {
            sumReal = 0;
            for (m = 0; m <= M - 1; ++m)
            {
                    testFreq = lowTestFreq + m * testFreqStep; // Hz.
                    testUnitAngle = (testFreq / samplFreq) * (2 * pi); // rad/sample.

                    Re = 0; Im = 0;

                    // Start DFT.
                   for (n = 0; n <= N - 1; ++n)
                    {
                            testAngle = n * testUnitAngle;

                            xn_cos = xn[n] * cos(testAngle);
                            xn_sin = xn[n] * -sin(testAngle);

                            Re += xn_cos;
                            Im += xn_sin;
                    }
                   // End DFT.

                   cosine = cos(((no * 2 * pi) / N) * m);
                    sine = sin(((no * 2 * pi) / N) * m);
                    real = Re * cosine - Im * sine; // Product of complex
//numbers (Re, Im i) * (cosine, sine i) = (Re * cosine + Re * sine i + Im i *
//cosine + Im i * sine i) = (Re * cosine + Re * sine i + Im i * cosine - Im *
//sine), for, by definition, i * i = -1. Then, first and last product in
//parenthesis are real numbers. Other two are imaginary and their sum is
//practically zero. For the values in array are between -26 and 25.
//More over because  input values are real, output should be real.
                    sumReal += real;
            }
            std::cout << sumReal / N << "\n";
    }
    return 0;

}


IFFT and sample count issue in OFDM I wrote a python program to generate OFDM packets (i.e. the IQ samples) with a custom payload to, for example, be transmitted via an SDR.
After mapping the chunk of bits corresponding to one ofdm symbol to a constellation mapping (e.g. 64QAM), I perform the IFFT (using numpy.fft.ifft) and add the 16-sample cyclic prefix.
In the IEEE standard, the resulting 80-sample time snippet is multiplied by a windowing function (x0.5 for first and last value). However, when I test the whole thing and perform FFT on the 64-sample snipped again the whole mapping has a rather large error.

If I omit the windowing its perfect. I can't really understand the heavy influence, since its only one single value that has been altered by the window function in the end (concerning the fft).

Now I've read in the standard that symbols are appended one after another with one sample overlap, which I am currently not doing, so I guess there lies my mistake. If I append them with an overlap of one though, I would end up with a wrong overal sample count, wouldn't i? ( i.e. $80 * nSymbols - nSymbols$ instead of $80 * nSymbols$)

Thanks for every hint!
Cheers Lucas

First value is from cyclic prefix and last value is the last sub-carrier which you are messing up by windowing of 1 sample. As a result of this, the frequency orthogonality of the sub-carriers is lost. And this inter-carrier interference is causing your QAM sliced symbols to randomly change phase and amplitude around true value.

The way actually windowing is done is by adding cyclic suffix of one sample in the end also like cyclic prefix and then overlap the next symbol with the first one by 1 sample. Then ICI won't happen. For a detailed code of how windowed overlap and add (wola) is performed you can check my code at : 5GNR Wola Waveform Simulation


How do they produce an impedance/phase plot for a speaker with a chirp that lasts 1 sec? Commercial software is available that sweeps a speaker in about 1 second and produces impedance/phase plot from 1 Hz to 20 kHz.
How do they do it? Obviously not by a single tone sweep over the full range.
So it must be a clever waveform that is frequency-stepped over some range while doing a DFT for each step. But what is this waveform? Or is there another way?
I can hear the sweep in the speaker and it is not an impulse 'click' but more of a chirp that only lasts 1 sec.

I need to build one for another purpose.

Thanks for any insights

Quoting Wikipedia article on Electrical impedance:

  Impedance is defined as the frequency domain ratio of the voltage to the current.

You can measure instantaneous current indirectly by connecting a small resistor (say 0.1 ohm) in series with the speaker and by measuring the instantaneous voltage across it, as function of time. Because the resistor is purely resistive, by Ohm's law $I = V/R$. To know the instantaneous voltage seen by the speaker, subtract the instantaneous voltage drop across the resistor from the instantaneous amplifier output voltage, which you can also measure as function of time.

Figure 1. Impedance curve test circuit. Voltage is to be measured at the test points TP2 and TP3 referenced to TP1.

In order to do the frequency domain division, play back the test signal, measure both instantaneous voltage and instantaneous current across the speaker as function of time, digitize those two signals, do a discrete Fourier transform on both, and complex divide bin by bin in the frequency domain. As long as there are no bins that are (too close to) zero-valued, the result should be meaningful. A short chirp, even a click, would work as long as it contains all frequencies of interest. The discrete Fourier transform should be long enough to fit the portion of the current signal where "things are happening". That can be significantly longer than the chirp, if it is short.

This kind of testing relies on the assumption that the system is linear. Some causes of nonlinearities (mainly at the speaker if the amplifier has good specs) can be avoided by limiting the peak amplitude of the test signal. A chirp is better test signal than an impulse, because it carries more power for the same peak amplitude, improving the signal-to-noise ratio.


what is the mel scale? I am not sure I understand what the Mel Scale is.

Googling doesn't give me various answers. I seem to be getting the same response again and again.

Which would be something like:

"The mel scale reflects how people hear musical tone"

First of all, does the concert pitch has a specific frequency? So why do we even have a scale stating how do people hear it?

I am currently trying to understand how MFCC uses them, and it seems like it maps the frequency to the Mel scale, or there are some filter banks which are applied on the given frames

To stand on @hotpaw2's answers, think of Mel as one kind of pyscho-acoustic scale, derived from a set of experiments on human subjects, others are Bark & ERB

Why have such a scale?

Imagine a 100 Hz sine wave playing in your head, .... wait for it to stick ...okay... Good.
Imagine now, a 200 Hz sine wave playing ... marinate in it. Good work.
Now compare the two, they have some "perceptual distance". Perceptual only because we're using your soft-grey matter as a measuring instrument.
Finally, if you were you repeat this experiment with 1100 Hz & 1200 Hz, and again with 10,100Hz and 10,200Hz, your opinion of this percpetual distance might diminish, i.e. 100 & 200 Hz can sound "farther apart" than do 10,100 & 10,200Hz, even though in all cases the differences are equivalent to 100Hz.

These pyschoacoustic scales try to capture these distances from low (20Hz) to high frequency (20kHz).


How to differentiate an empty aerosol from a full one using only sound measurements? The challenge I have is to detect an empty spray container using sound measurements only, since I can't place any sensor directly into the container because of the mechanism that tightens it.
My main idea was to calculate the FFT of a 2-second recording just when the mechanism squeezes the trigger of the sprayer and then get the maximum peak frequency.
However, I don't know if it is appropriate to use peak detection because when it is empty, it does not emit sound.

Any ideas for solving this ?

What I would do is :

1) Acquire the audio signal (or in fact, several, like ten) corresponding to a full container. I'd then quantify the energy of the signal, by squaring it and dividing by the number of samples. The obtained quantity is my reference energy. If several signals were measured, I take the average of calculated energies. Let's call it $E_{full}$

2) Do the same with an empty container (which is basically measuring ambient noise if it does not emit any sound like you say). Let's call it $E_{empty}$

From this point, you can define any energy threshold, with which you will compare the energy of the sound emitted by your experimental aerosol.

You can also do it simpler, by detecting for example the maximum peak amplitude of the full and empty aerosol and defining peak-amplitude threshold accordingly.

The threshold definitions can also be fine-tuned by measuring other aerosols with varying fillings, like half-empty, 25 % etc...

And finally, all the measurements I suggest implies that the sound is measured in steady state.


Kalman Filter - Velocity [Matlab] I have a quite typical Kalman filter to design. I really read a lot of articles about the design of this filter but the performances of my filter are still quite bad.

Here is my situation. I have a quite good measurement signal of my position (let's say a very small white noise) and a pretty noisy measurement signal of my velocity (big white noise). I want to estimate a new position (which will be I guess not very different from my measurements) and my velocity as well (which will be more different)

Here is my Matlab code: (I don't want to use the Matlab Kalman function ;) )

%% Kalman Filter Design
dt = 0.01;

%% Define coefficient matrices
% x = [position;velocity;acceleration]

% True State equation
% x = F * x(t-1) + B * u + w
F = [ 1 dt dt^2/2
      0  1 dt
      0  0 1];

B = [0
     0
     0];
% Measurement of the true state equation
% z = H * x(t) + v
H = [1 0 0
     0 1 0];

%% Define noise

%Process noise(white noise)
% w = Gw * a
a = [0 ; 0 ; wgn(1,1,100)];
Gw = [0 0 0
     0 0 0
     0 0 1];

%Process noise covariance matrix
Q = Gw * Gw' * cov(a);

%Measurement noise (white noise)
v = [wgn(1,1,0.001) ; wgn(1,1,20)];

%Measurement noise covariance matrix
R = cov(v)* [0.5 0
              0  1];

%% Kalman Filter
x_estimate = [0;0;0];
P = Q;

position_estimate = [];
velocity_estimate = [];
acc_estimate = [];
P_mag_estimate = [];

predic_state =[];
predic_var = [];

z = [position_meas,velocity_meas];

for t = 1:length (z)

%Predicted state estimate
x_estimate = F * x_estimate;
predic_state = [predic_state; x_estimate(1)];

%Predicted estimate covariance
P = F * P * F' + Q;
predic_var = [predic_var; P];

%Innovation covariance
S = H * P * H' + R ;

%Kalman Gain Predict measurement covariance
K = P * H' * inv(S);

% Updated  state estimate
x_estimate = x_estimate + K * (z(t) - H * x_estimate);

%Update covariance estimation
P = (eye(size(P,2)) - K * H) * P;

%Store for plotting
position_estimate = [position_estimate; x_estimate(1)];
velocity_estimate = [velocity_estimate; x_estimate(2)];
P_mag_estimate = [P_mag_estimate; P(1)];

end

My question is: how can I find Q and R? Does it depend of what I want to do?

Thanks a lot for your help

I think I've managed to implement your system, but in R rather than matlab (sorry, I have no access to matlab and scilab stopped working for me on OSX).

The implementation gives the following estimates of the position and velocity, and seems to work even with the large velocity noise.

The top plot is the true (red), noisy (black), and estimated (blue) position. The bottom plot is the true (red), noisy (black), and estimated (blue) velocity.  I've zoomed in on the bottom plot to see the true and estimated values because the noise is so large.

To your questions:

  My question is: how can I find Q and R?

Well, from my code Q is

  Q <- matrix(c(0,0,0,0,0,0,0,0,1),3,3)

which would be

  Q = [0,0,0; 0,0,0 ; 0,0,1];

in matlab, I think.

And R is

  R <- matrix(c(0.001,0,0,20),2,2)

which is

  R = [0.001 0; 0 20];

in matlab.

I've generated my acceleration as a vector that is T long with a variance of 1; you seem to have generated just one value for it with a variance of 100.

As for the second part of your question:

  Does it depend of what I want to do?

Generally, the better the values the Kalman filter use match the "true" values, the better the Kalman filter estimates the state.  I say "true" instead of true because sometimes we can't know what the truth is, so we have to guess it.

The only leeway I see in what you've presented is what the value if Q is. It's not clear to me what variance you assumed the acceleration to have.  If it is a constant, then the variance will be zero; if it's a Gaussian random variable, then you'll need to set the Q matrix using the variance of that random variable.

Also, this all assumes that the process noise and the measurement noise are uncorrelated (or independent). If there is correlation, then you'll somehow need to take that into account.

R Code

  # 26818
  T <- 128;
  dt <- 0.01

  t <- seq(0,T-1)*dt

  xkm1km1 <- matrix(c(0, 0, 0),3,1)
  Pkm1km1 <- matrix(c(1000,0,0 ,0,1000,0, 0,0,1000),3,3)
  #H = [1 0 0
  #     0 1 0];
  H <- matrix(c(1,0,0,1,0,0),c(2,3))

  #F = [ 1 dt dt^2/2
  #      0  1 dt
  #      0  0 1];
  A <- t(matrix(c(1,dt,dt^2/2,0,1,dt,0,0,0),3,3))
  B <- matrix(c(0,0,1),c(3,1))

  true_acceleration <- rnorm(T,0,1)
  measurement <- t(matrix(rep(0,T*2),c(2,T)))
  x <- t(matrix(rep(0,T*3),c(3,T)))

  position_noise <- rnorm(T,0,0.001)
  velocity_noise <- rnorm(T,0,20)

  v <- t(array(c(position_noise, velocity_noise),c(T,2)))

  Q <- matrix(c(0,0,0,0,0,0,0,0,1),3,3)
  R <- matrix(c(0.001,0,0,20),2,2)

  # Generate the data using the signal model
  for (k in 2:T)
  {
    x[k,] = A %*% x[k-1,] + B * true_acceleration[k-1]
    measurement[k,] = H %*% x[k,] + v[,k]
  }

  library("MASS") # For pseudo inverse ginv()

  zhat <- matrix(rep(0,2*T),c(T,2))

  for (k in 1:T)
  {
    xkkm1 <- A %*% xkm1km1
    Pkkm1 <- A %*% Pkm1km1 %*% t(A) + Q
    K <- Pkkm1 %*% t(H) %*% ginv( H %*% Pkkm1 %*% t(H) + R)
    z <- matrix(c(measurement[k,]), 2, 1)
    xkm1km1 <- xkkm1 + K %*% (z - H %*% xkkm1)
    Pkm1km1 <- (matrix(c(1,0,0,0,1,0,0,0,1),3,3) - K %*% H) %*% Pkkm1
    zhat[k,] <- as.numeric(H %*% xkkm1)
  }

  par(mfrow=c(2,1))
  plot(t, measurement[,1], type="b", pch=19)
  lines(t, x[,1], col="red",  lwd=10)
  lines(t, zhat[,1], col="blue", lwd=5)

  plot(t, measurement[,2], type="b", pch=19, ylim=c(-1,1))
  lines(t, x[,2], col="red",  lwd=10)
  lines(t, zhat[,2], col="blue", lwd=5)

  errs <- c(sum((signal - zhat)^2), sum((signal - lpf)^2))
  print(errs)


Implementation of a Variable Fractional Delay with Lagrange Interpolation using Farrow Structure I am writing a C++ simulation software working in time domain. I generate regularly sampled data, and need them to be delayed, in "real-time", by a variable fractional delay.

This is a pretty common problem, and very close to all the "resampling" issues. After a short bibliographic search, it seems that the most commonly used method is Lagrange Interpolation with a Farrow structure. Based on J. O. Smith, "Interpolated Delay Lines, Ideal Bandlimited Interpolation, and Fractional Delay Filter Design", I used the following implementation:

which includes cascading filters with fixed coefficients $1+z^{-1}$, modulated by cascading factors that depend on the variable delay. This can be seen in the equation

$$
\hat{H}_\Delta(q^{-1}) = \sum_{n=0}^N \hat{h}_\Delta[n](q-1)^n.
$$

My implementation in C++ works well for interpolation orders up to 51.
Unfortunately, for larger interpolation orders the numerical errors increase rapidly and become dominant. After a quick inspection, it seems that the buffered samples, in the $1+z^{-1}$ cells of increasing index, have an increasing order of magnitude. It helps when I use long double, but this becomes insufficient when the order is of the order of 60: we loose a lot of numerical precision when we combine floating-point values of around $10^{10}$ with values of the order of $1$.

Do you know a variation of the Farrow structure or another implementation for Lagrange interpolation that works well for large orders, i.e. that remains computationally efficient while requiring only reasonable numerical precision in the memory buffers?

Thank you in advance for your time :)

The task at hand is a pretty common one; most symbol timing synchronizers need to do exactly what you describe.

Same goes for many types of resamplers – they need to reconstruct varying fractional inter-sample "snapshots" from the input sample stream to construct a non-rationally related rate output.

So that's where I'd start looking.

In fact, the farrow structure idea takes you pretty far!

Imagine you could have a large bank of linear-phase filters. A linear-phase filter has a constant group delay.

By modulating said filters, you can change the phase slope and hence the group delay. Now, you just pick the filter that fits your needs.

In a polyphase filterbank, such as implemented by the GNU Radio "PFB Arbitrary Resampler" block, exactly that is done: You construct a whole bunch of filters, and for each sample, pick the one that's closest to what you need (or the two closest, and linearly interpolate between these).

Because it's used as a resampler, for every non-1.0 resampling ratio, the next sample will be choosen from a different filter branch, with the appropriate delay.

Now, since you, most of the time, just want to delay, you seldom need the "resampling" capability of that resampler. (Resampling is exactly what happens when you change the delay linearly)

So, instead, you'd have a single filter (if possible, use a filter that you need to have, anyway), and you'd modulate it.

Think about this in frequency domain: assume you had a linear, continuous Fourier transform that you apply to your signal and you multiply that with both the FT of your filter impulse response and $e^{j2\pi\tau f}$ before transforming back to time domain. To adjust delay, you adjust $\tau$. Tadaaah! Variable fractional delay.

Problem certainly is that you're in DSP, so you don't get that linear Fourier transform to implement linear convolution; you get the DFT, but multiplication in discrete frequency domain is equivalent to circular convolution in time domain.

But that's been solved already – fast convolution allows you to apply long filters with little effort (and low numerical problems – the FFT is exceptionally stable, usually even more stable than the direct convolution!). Now, if you go ahead and implement that e.g. via the overlap-save method, all you need to do is make sure that the phase of the complex exponential doesn't "jump" between FFT segments. You can then change (even gracefully!) the $\tau$ between individual convolutions.


How to Map CWT to Synchrosqueezed wavelet transform? I don't understand the mapping time-scale plane to the time-frequency plane in synchrosqueezed wavelet transform, i.e. $(3)$. You can find the paper here.

For the given signal of $x(t)$ and mother wavelet of $\psi$, the continous wavalet transform is:

$$W(u,s)= \int_{-\infty}^{+\infty}  x(t)\psi_{u,s}^*(t)dt\tag{1}\\$$

From $(1)$, a preliminary frequency $\omega (u, s)$ is obtained from the oscillatory behavior of $W (u, s)$ in $u$:

$$\omega(u, s) = − i\left(W(u, s)\right)^{−1} \frac{\partial }{\partial u} W(u, s)\tag{2}$$

$W(u,s)$ is then transformed from the time-scale plane $(1)$ to the time-frequency plane $(3)$. Each value of $W(u, s)$ is reassigned to $(u, \omega_{l} )$, where $\omega_{l}$ is the frequency that is the closest to the preliminary frequency of the original (discrete) point $\omega(u, s)$.

$$T(u, \omega_{l}) = \left(\Delta\omega\right)^{-1} \sum_{sk:|\omega(u, s_{k} )−\omega_{l}|\leq\Delta\omega/2}^{} {W(u, s_{k})s_{k}^{−3/2}\Delta s}\tag{3}$$

Let me explain the intuition briefly.  The authors of the paper you've cited assume that the signal $x(t)$ can be written in the form
\begin{align*}
x(t) &= \sum_{k=1}^K a_k(t) \exp(2\pi\mathrm{i} \phi_k(t)),
\end{align*}
where the $a_k$ denote instantaneous amplitudes, the $\phi_k$ denote instantaneous phases and the $\phi'_k$ denote instantaneous frequencies. The continuous wavelet transform enables the visualization of the $\phi'_k$ in the time-scale plane, but the visualization is necessarily blurry due to the Fourier uncertainty principle.  The idea is that the peaks of the continuous wavelet transform can give you an idea of what the $\phi'_k$ values are, but the energy is spread out in the time-scale plane, so the exact values of $\phi'_k$ may not be easy to obtain.  The goal of the synchrosqueezing transform (SST) is to partially undo the blurring by a frequency reassignment, which can be explained as follows.

First of all, let's note that each scale $s$ corresponds to a natural frequency $\xi$ satisfying the relation $s = c/\xi$, where $c$ is the center frequency of the mother wavelet $\psi$.  Now, suppose the time $u$ is fixed.  If $\xi = c/s$ is close but not exactly equal to an instantaneous frequency $\phi'_k(u)$, then the coefficient $W(u,s)$ will have some nonzero energy (i.e., $|W(u,s)|^2 > 0$).  The idea of synchrosqueezing is to move all this energy away from the frequency $\xi$, and reassign its frequency location closer to the instantaneous frequency $\phi'_k(u)$.  Then, we can arrive at a time-frequency representation where the energy is more closely concentrated around the instantaneous frequency curves.

So, for each scale $s$, we compute the frequency reassignment  $\omega(u,s)$ of the wavelet transform coefficient $W(u,s)$ by the formula you cited above.  As Daubechies, Wu, and Lu discovered in the paper you've linked, it turns out that $\omega(u,s)$ is often a good approximation to the instantaneous frequency curve $\phi'_k(u)$ for values of $(u,s)$ where $\xi = c/s$ is sufficiently close to $\phi'_k(u)$. (In fact, when $\phi'_k(u)$ is a constant function, $\omega(u,s)$ exactly equals $\phi'_k(u)$ for values of $s$ where $\xi = c/s$ is sufficiently close to $\phi'_k(u)$!)

Then, the computation of the SST in the continuous setting is as follows. First, fix a time of interest $u$. Next, compute the frequency reassignment $\omega(u,s)$ for all scale values $s$. Then, for each frequency of interest $\eta$, we compute the SST $T(u,\eta)$ by adding up all values $W(u,s)$ where the reassigned frequency $\omega(u,s)$ equals $\eta$. This means

\begin{align*}
T(u,\eta) &= \int_\mathbb{R} W(u,s) \delta(\eta - \omega(u,s)) s^{-3/2} ds,
\end{align*}

where $\delta$ is the Dirac delta. (To ensure the convergence result for SST, the authors use a formulation for $T$ which includes an approximation to $\delta$ rather than $\delta$ itself, but in computational practice I've never found this approximation to be necessary.) The multiplication by $s^{-3/2}$ is necessary for reconstruction purposes (which you can read more about in the paper).

In practice, one is limited to the discrete setting where we have only finitely many possible frequency bins $\eta_\ell$. In this case, for a fixed frequency $\eta_{\ell_0}$, one finds all the values of $\omega(u, s)$ which are closer to $\eta_{\ell_0}$ than to any other frequency bin $\eta_\ell$. We can more explicitly formulate the SST in the discrete setting as follows.

For simplicity's sake let's assume the possible frequency values $\eta_\ell$ are uniformly spaced by a distance $\Delta \omega$. So, in the discrete setting, the formula above becomes
\begin{align*}
T(u, \eta_\ell) &= \sum_{s: |\omega(u,s) - \eta_\ell| < \Delta\omega/2} W(u,s) s^{-3/2},
\end{align*}
where I haven't bothered to give an explicit notation that shows that the values $u$ and $s$ are discrete.

By doing this reassignment, you end up with a representation $T$ which is sparser than $W$, and hopefully very sharply concentrated about the curves $\phi'_k$. (In practice, for visualization purposes, it's better to reassign not the original coefficient $W(u,s)$ but the magnitude-squared $|W(u,s)|^2$, or even just the constant number $1$. This is because the coefficients $W(u,s)$ are complex-valued, and summing these coefficients may not actually leave you with a very large energy. However, if you additionally want to use the reconstruction formula given in the paper to compute the $x_k$ from the SST, then you need to sum the $W(u,s)$ in order for that formula to work.)

By the way, I am leaving out the fact that we generally throw out coefficients $W(u,s)$ that do not fall above a certain threshold.  This is because the computation of $\omega(u,s)$ is not necessarily accurate for small $W(u,s)$.


Comparison of Linear Convolution and N Point DFT I am studying Discrete Fourier Transform currently and I have a doubt in that.

Consider two sequences $$\left\{1,2,3,4\right\}$$ and $$\left\{0,1,0,0\right\}$$ When I convolve them linearly, I get this $$\left\{0,1,2,3,4,0,0\right\}$$

However, if I take the 4 point DFT of $\{1,2,3,4\}$ and $\{0,1,0,0\}$ and then multiply them, and after that if I take IDFT of the result, I get this $$\{4, 1, 2, 3\}$$

$$\mathrm{DFT}({1,2,3,4}) = \{10, -2+2j, -2, -2-2j\}$$
$$\mathrm{DFT}({0,1,0,0}) = \{1, -j, -1, j\}$$

Element wise product $$\{10, 2+2j, 2, 2-2j\}$$
$$\mathrm{IDFT}({10, 2+2j, 2, 2-2j}) = \{4, 1, 2, 3\}$$

I know that the output of linear convolution here would need at least 7 elements and I am calculating the 4 point DFT/IDFT which of course would not be enough.

But what inference should I take from this. The result of IDFT looks somewhat like circular convolution but it is not exactly the result that you would get even after circular convolution.

The result {4,1,2,3} is the circular convolution of {1,2,3,4} and {0,1,0,0} which you correctly get by taking the inverse DFT of the product of the DFTs of the two sequences.

We can check this by doing the circular convolution the long way via matrix multiplication as follows:

$x(n) = [1,2,3,4]$

$y(n) = [0,1,0,0]$

To solve for the circular convolution $x(n) (*) y(n)$, (where I use $(*)$ to denote circular connvolution) we can do that with multiplication as $Xy$, where $X$ is a matrix formed by repeating $x(n)$ in each column with a circular rotation of shift n with n = 0,1,2,3 and $y$ is $y(n)$ as a column vector:

$$ Xy =\begin{bmatrix}
1 & 4 & 3 &2 \\
2 & 1 & 4 & 3 \\
3 & 2 & 1 & 4 \\
4 & 3 & 2 & 1 \\
\end{bmatrix} \begin{bmatrix}0 \\ 1 \\ 0 \\ 0 \end{bmatrix} = [4,1,2,3]$$

Confirming the relatinship:

$$x[n](*)y[n] = IFFT(FFT(a)FFT(b))$$

Where the ($*$) operator here is specifically circular convolution. And the multiplication of the two FFT results is the element by element product.

Also to add to this: In comparison, if you instead conjugate one of the FFT results prior to multiplication as follows, the result is a circular cross-correlation of the two sequences!

$$XCORR = IFFT(FFT(a)FFT(b)^*)$$


Autocorrelation of Addition of Two Independent Signals Given a random signal $ Z \left( t \right) $ which is addition of two independent signals $ X \left( t \right) $ and $ Y \left( t \right) $ with constant parameters $ a $ and $ b $:
$$ Z (t) = aX(t) + bY (t) $$
If the auto correlation function of $ R_{XX} \left( \tau \right) $ and $ R_{YY} \left( \tau \right) $ is known (Assume both are Stationary Signals).
What would be the Auto correlation function of $ Z \left( t \right) $?
Is it given by (I skipped the steps):
$$ {R}_{ZZ} \left( \tau \right) = a^2 R_{XX} \left( \tau \right) + b^2 R_{YY}\left( \tau \right) $$
I neglected the $ {R}_{XY} \left( \tau \right) $ contribution because it is zero as $ X \left( t \right) \perp Y \left( t \right) $ (Independent Random Signals).

You're correct as the Cross Correlation function vanishes.
This has the implicit assumption the process has zero mean (Actually, at least one of them).

Namely, in order to have $ {R}_{XY} \left( \tau \right) = 0 $ having $ X \left( t \right) \perp Y \left( t \right) $ isn't enough but at least of them has zero mean (Namely, $ \mathbb{E} \left[ X \left( t \right) \right] = 0 $, $ \mathbb{E} \left[ Y \left( t \right) \right] = 0 $ or both).

Pay attention that in practice for a single realization the empirical cross correlation won't vanish.


number of possible component in sinusoidal model Suppose that we have the following model

$y(t) = A_1\sin(\omega_1 t+\phi_1) + A_2\sin(\omega_2 t+\phi_2) + ... + A_p\sin(\omega_p t+\phi_p) + z(t)$

My question is not related to how to determine the number of components, but I am interested in the maximum. For example, let us assume that we have sampled this model and get sample with size let's say $294$, the question arises what can be the maximum number of components? For instance, can it be more then $147$? or more then half the size of the signal sample? Or can it be in general something $20-50$, for instance? I am not expecting an exact number, but a possible maximum number which can be in this case.

I guess that your problem can be modeled as a harmonic retrieval problem.
You can check the Kung's method.


The fastest way to "fill holes" in 3D? I have a 3D binary image that I must fill all of the holes in the image.

"Fill holes 3D" in this case, is basically a 3D version of

BW2 = imfill(BW,'holes') from MatLab,

but I am not focusing on Matlab here. I am focusing on the algorithm in general.

Currently, I simply perform a "flood fill 3D" algorithm on all voxels at the border area of the 3D matrix. Anything not filled are holes so take this as a mask and remove all voxels on the real image with this mask.

Sure, the algorithm works, but somehow, I feel this is inefficient.

What is the fastest way to "fill holes" in 3D in this case?

I don't have the privilege to comment, so I'm going write my comment as an answer.

I'm guessing your function, as well as imfill in Matlab, performs the "filling" iteratively.

Two basic improvements to speed these kind of iterative functions up would be to 1. increase the seed points, i.e. using other possible locations in addition to the border area, maybe base on image intensity before the binerization process, or whatever characteristics the images have, and 2. try to parallelize your function (see Dirk-Jan Kroon's Region Growing function for an example).  If this doesn't speed things up, or you just have too many seed points, then I would suggest you start implementing a function that uses your GPU for parallel processing.

Now there may be "non-iterative" functions that fill 3d holes... so I'll leave that to the experts here on dsp.se.  It would also be nice if you can show us a slice of your 3D image for a better understanding of what you are dealing with.


Kalman filter - understanding the noise covariance matrix What is the significance of the noise covariance matrices in the Kalman Filter framework?

I am referring to:

process noise covariance matrix Q, and
measurement noise covariance matrix R

at any time step t.

How do I interpret these matrices? What do they represent? Do they talk about how one observation's noise varies with respect to another observation's noise in the state vector?

Roughly speaking, they are the amount of noise in your system. Process noise is the noise in the process - if the system is a moving car on the interstate on cruise control, there will be slight variations in the speed due to bumps, hills, winds, and so on. Q tells how much variance and covariance there is. The diagonal of Q contains the variance of each state variable, and off diagonal contain the covariances between the different state variables (e.g. velocity in x vs position in y).

R contains the variance of your measurement. In the above example, our measurement might just be speed from the speedometer. Suppose it's reading has a standard deviation of 0.2 mph. Then R=[0.2^2]=[0.04]. Squared because variance is the square of the standard deviation.

Q is in state space, and R is in measurement space. In the example above, our state might be position only $[x, y]^T$, and measurement space is velocity $[v]$. That is problematic because that is not velocity in terms of x and y - you need the heading to convert. The Kalman filter matrix H is used to do that conversion, and in nonlinear systems you tend to have to linearize that in some manner.

Shameless plug: my free book on the Kalman filter goes into this in great detail: https://github.com/rlabbe/Kalman-and-Bayesian-Filters-in-Python


Is $y[k] = y[k-1] + x[k]$ an integrator? It looks exactly like an integrator to me.
Since $$y[k] = y[k-1]+x[k] = y[k-2] + x[k-1] +x[k] = \sum{x}$$

Applying the Z-transform gives
\begin{align}
Y(z) &= Y(z)\cdot z^{-1} + X(z)\\
\Rightarrow\frac{Y(z)}{X(z)} &= \frac{1}{1-z^{-1}}
\end{align}

When I convert $\frac{1}{1-z^{-1}}$ into Lapace transform in Matlab using d2c, it does not return $\frac{k}{s}$. Instead, it returns $\frac{1+s}{s}$. It seems there is a proportional term inside, when thinking in PID control point of view.

Could anyone explain on that?

The system

$$y[n]=y[n-1]+x[n]\tag{1}$$

is an ideal accumulator, i.e., it computes the cumulative sum of the input samples:

$$y[n]=\sum_{k=-\infty}^nx[k]\tag{2}$$

It is in a way analogous to a continuous-time integrator, but this doesn't mean that you will necessarily obtain an ideal integrator by transforming the discrete-time system to a continuous-time system. There are several methods to do the conversion, and I'm not sure which one you used when calling d2c.

In any case, the properties of the transformed continuous-time system not only reflect the properties of the original discrete-time system, but also the properties of the transformation. So the properties of a discrete-time system should be investigated directly in the discrete-time domain, not by transforming it to a continuous-time system.


Convolutions with changes to the argument I think I understand what happens when I shift the argument, but I'm not sure what should happen when the signal is compressed or expanded. In particular I'm trying to figure out what happens when the convolution where $y(t)=x(t)*h(t)$ is changed to $x(t/k)*h(t/k)$ or $x(t/k)*h(t)$. Supposedly the first should be $y(kt)/k$ and the second should be unable to be done but I don't know how to get there. When I work the first one out, I get $\int{x(\frac{\tau}{k})h(\frac{t-\tau}{k})}$. This doesn't look like something I can get by changing the argument of $y$, so I must be missing something. Any help on how to go about doing these convolutions with arguments other than t would be helpful.

Let me show you how to manipulate the convolution integral.
Let
$$y(t)  = x(t) \star h(t) = \int_{-\infty}^{\infty} x(\tau) h(t-\tau) d\tau \tag{1} $$
Let $a > 1$ that we have a compressed signal $y(at)$ ;
$$y(at)  = \int_{\tau=-\infty}^{\infty} x(\tau) h(at-\tau) d\tau \tag{2} $$
Substitude $\tau' = \tau/a$ so that $\tau = a \tau'$ and $d\tau = a \cdot d\tau'$ into Eq.2 :
\begin{align}
y(at)  &= \int_{\tau'=-\infty}^{\infty} x(a \tau')~ h(at-a\tau') ~a ~d\tau' \tag{3}\\ \\
&= a \int_{\tau'=-\infty}^{\infty} x(a \tau') ~ h(a(t-\tau')) ~d\tau' \tag{4}\\ \\
&= a \int_{\tau=-\infty}^{\infty} x(a \tau) ~ h(a(t-\tau)) ~d\tau  = a \cdot \big( x(at) \star h(at) \big) \tag{5}
\end{align}
Eq.5 shows that the compressed output $y(at)$ is obtained by the convolution of the compressed input $x(at)$ and the compressed impulse response $h(at)$, and (also the result scaled by $a$)....
Note that convolution operator assumes an LTI (linear time-invariant) framework. However, the compressor is not an LTI system and it does not have an impulse response $h_c(t)$, so that $y(at) = y(t) \star h_c(t)$ can not be true.
Assume for a moment that the compressor is an LTI system with impulse response $h_c(t)$. Then we can consider the generation of the compressed output in two (serial) LTI stages:
$$ x(t) \longrightarrow \boxed{ h(t) } \longrightarrow \boxed{ h_c(t)} \longrightarrow y(at) \tag{6}$$
From LTI theory we know that the cascade of two LTI blocks is also LTI with an equivalent impulse response given by $h_e(t) = h(t) \star h_c(t)$ hence Eq.6 can be written as :
$$ y(at) = x(t) \star h_e(t)  = x(t) \star ( h(t) \star h_c(t)) = x(t) \star h(at)\tag{7}$$
Which is contradicting the result we have obtained in Eq.5, hence the assumption that the compressor could have been represented with an impulse reponse was wrong.
Another implication of this is also observed in the time-reversal of the output $y(t)$ :
$$y(-t) = x(-t) \star h(-t) \neq x(-t) \star h(t) \neq x(t) \star h(-t) $$
The reason of the unequal signs is again about reversal operation not being an LTI system.
On the other hand, if you want to apply an LTI operation on the output such as $y_d(t) = y(t-d)$ , then it can be simply be reflected in the convolution as:
$$ y(t-d) = y(t) \star \delta(t-d) = (x(t) \star h(t)) \star \delta(t-d) = x(t) \star h(t-d)$$...


Is Phase Locked Loop is essential for BPSK signal reception? I'm noob in DSP studying BPSK communication between speaker and microphone using acoustic signal. I read article of this BPSK communication which says it needs Phase Locked Loop(PLL) sequence to match both ends communicate precisely. But in my textbook, there is no mention of PLL in BPSK demodulation.

I want to communicate Tx/Rx within close range to get channel tap corresponding to its multipath effects in my room. In this case, does PLL is essential to synchronize the signal in both ends? Or I can communicate well without PLL and can get channel tap?

There are three different kinds of synchronization in a passband digital communications system:

Carrier synchronization: the receiver needs to know the exact frequency and phase of the carrier used by the transmitter.
Symbol synchronization: the receiver needs to know the optimum instants to sample the matched filter's output.
Frame synchronization: the receiver needs to know where each transmitted "character" (or word, or byte) starts.

So, if you're doing passband BPSK, you'll need to perform all three. Carrier synchronization is typically done with a PLL (or a Costas loop) and you can find the details in several answers on this website.

Many textbooks ignore these subjects. A textbook that does a great job of teaching these three aspects of a receiver is "Telecommunications Breakdown" by Johnson and Sethares, an early version of which is available for free from the authors' website.

Having said that, there is a way to avoid doing carrier synchronization, at the cost of energy efficiency, by upconverting the baseband BPSK signal using AM DSB-LC instead of the more common AM DSB-SC. So, you first design your baseband signal: $$s_{BB}(t) = \sum_k a_k p(t-kT_P),$$ where $a_k \in \lbrace -1, 1 \rbrace$, $T_p$ is the pulse interval, and $p(t)$ is the pulse shape. Then, the DSB-LC passband signal is $$s_{PB}(t) = (A_m + s_{BB}(t)) A_c \cos(2\pi f_c t),$$ where $A_m > -\min(s_{BB}(t))$ (so that $(A_m + s_{BB}(t))>0$) and $A_c$ is the carrier amplitude. Then, the baseband signal can easily be recovered by a simple envelope detector (identical to conventional analog modulation).


IFFT and OFDM upconversion I am trying to understand some basic concepts with IFFT as represented here OFDM Modem page 50-53, and 108-114.

Some basics that are represented -

The base sampling rate is 8khz.
The symbol period is 32ms, which would calculate to 250 total samples.
The fft size is 64.
The bin separations are 32hz.
The bandwidth used is from 1000-3000hz(1024-3072?)

I cannot get my mind around how to get from a ifft size of 64 to an output bandwidth of 1000-3000 and 32khz spacing.  Is there some basic reason why I wouldn't use a much larger ifft size of something like 256, and then use only 32-96? That would calculate to 31.25hz per bin.

In short, this thesis is proposing that you use a fft size of 64 to cover a bandwidth of 2khz beginning at 1khz, 1k-3k.  This isn't coming together for me at all, how this can get shifted.

If I have to be math'ed here, be gentle.

First attempt at upconversion using a spreadsheet:

=(C12*ROW(C12)*EXP(-1*2*PI()*10*ROW(C12)))/ROW(C12)

Where 10 is the carrier freq, C(12) is a row containing a sine lookup value.

To explain the supposed contradiction between the number of samples calculated by symbol length and sampling rate on the one hand and the IFFT size on the other hand, you have to take into account that the IFFT output is upsampled by a factor of 4. This means that for every output sample of the IFFT block there are four samples after upsampling. There's a cyclic prefix of 16 samples so that we arrive at a total number of samples of $N=4\cdot(64+16)=320$ as stated on page 53.

On the other hand we can calculate the number of samples by $N=Tf_\mathrm{S}=32\mathrm{ms}\cdot8\mathrm{kHz}=256$. The symbol duration $T$ has been taken from page 108. It think the guard interval has not been taken into account there and thus $N=64\cdot4=256$.

The process that is described as "upconversion" is exactly the frequency shift that you are observing (cf. figure 3-8). Citing from the report:

  Once the baseband signal has been interpolated it is multiplied with
  the audio band carrier frequency to obtain a audioband representation.

To make this more clear using a formula, let $x(t)$ be the complex, analog OFDM signal (after digital-to-analog conversion (DAC)) and let $f_0$ be the desired carrier frequency. Then the radio frequency, upconverted signal $y(t)$ is
$$
y(t) = x(t)\operatorname{exp}(j2\pi f_0 t)
$$
This is a basic concept required in almost every radio system. For more information you might want to do some reserach on radio frequency upconversion and I/Q modulation.

The baseband bandwidth of an OFDM signal is approximately $f_\mathrm{S}/2$. But again, here we have to take the upsampling into account. If it is done correctly, it won't change the bandwidth of the signal. I.e. the sampling frequency before upsampling is actually $f_\mathrm{S}/4 = 2\mathrm{kHz}$ and we have to use this frequency for calcuating the OFDM signal bandwidth wich here is $B=1\mathrm{kHz}$. Finally, the band occupied by the radio frequency signal is $f_0 - B\ldots f_0+B$ which corresponds to 1kHz to 3kHz.

Concerning your suggestion to use a greater IFFT size and to left some subcarriers unmodulated to achieve upconversion - this is a valid approach in this specific example. Note however, that in general it requires a larger sampling frequency which is often undesirable.


Kalman for 3D position and 1D orientation Is the following state correct for a Constant Acceleration (CA) model KF applied for tracking an object which moves in (x, y, z) and have the freedom of rotation $\phi$ around the $z$ axis?

\begin{equation}
 x=
 \begin{bmatrix}
 x, \dot x, \ddot x, y, \dot y, \ddot y, z, \dot z, \ddot z, \phi, \dot \phi, \ddot \phi
 \end{bmatrix}
\end{equation}

UPDATE:

I can only measure $(x, y, z)$, directly. After, I computed $\phi$, rotation around the $z$ axis, using $\arctan(\frac{\delta y}{\delta x})$, indirectly. The measurements that I input to KF are: $x, y, z, \phi$.

Assuming:
\begin{equation}
\pmb{\alpha}=
\begin{bmatrix}
1  & dt & 0.5 \times dt^2 \\
0  & 1  & dt \\
0  & 0  & 1
\end{bmatrix}_{3 \times 3}
\end{equation}

\begin{equation}
\pmb{\beta}=
\begin{bmatrix}
1  & 0 & 0
\end{bmatrix}_{1 \times 3}
\end{equation}

$A$, the state transition matrix which applies the effect of each state parameter at time $t-1$ on the state at time $t$, is built using:
\begin{equation}
x_{k+1}=
\begin{bmatrix}
\pmb{\alpha}_{3 \times 3}  & 0_{3 \times 3} & 0_{3 \times 3} & 0_{3 \times 3} \\
0_{3 \times 3}  & \pmb{\alpha}_{3 \times 3} & 0_{3 \times 3} & 0_{3 \times 3} \\
0_{3 \times 3}  & 0_{3 \times 3} & \pmb{\alpha}_{3 \times 3} & 0_{3 \times 3} \\
0_{3 \times 3}  & 0_{3 \times 3} & 0_{3 \times 3} & \pmb{\alpha}_{3 \times 3} \\
\end{bmatrix}_{12 \times 12}
\cdot
\begin{bmatrix}
x\\ \dot x\\ \ddot x\\ y\\ \dot y\\ \ddot y\\ z\\ \dot z\\ \ddot z\\ \phi\\ \dot \phi\\ \ddot \phi
\end{bmatrix}_{k}
\end{equation}

In the same manner, the transformation matrix $H$, which map state vector parameters into the measurements domain, is built using:
\begin{equation}
H=
\begin{bmatrix}
\pmb{\beta}_{1 \times 3}  & 0_{1 \times 3} & 0_{1 \times 3} & 0_{1 \times 3} \\
0_{1 \times 3}  & \pmb{\beta}_{1 \times 3} & 0_{1 \times 3} & 0_{1 \times 3} \\
0_{1 \times 3}  & 0_{1 \times 3} & \pmb{\beta}_{1 \times 3} & 0_{1 \times 3} \\
0_{1 \times 3}  & 0_{1 \times 3} & 0_{1 \times 3} & \pmb{\beta}_{1 \times 3} \\
\end{bmatrix}_{4 \times 12}
\end{equation}

Reading your link in more detail, it says:

  Acceleration ($\ddot{x}$ & $\ddot{y}$) as well as position ($x$ & $y$) is measured.

which, to me, is not a constant acceleration.  What it is doing is a discretized changing acceleration model, assuming that the acceleration is constant between each sample period.

Now to your question: is inclusion of $\phi, \dot \phi,$ and $\ddot \phi$ in the state vector valid?

The answer is yes, but you need to include the dynamics of how those states interact with the other states.

OK, so now you've added the extra information, I can take a stab at it.

It should work with the following provisos:

The noise distribution on the angle $\phi$ will be uniform not Gaussian.  I can't find a good reference, but this discussion on comp.dsp mentions it.  That means the KF equations for $\phi$ won't be quite right.  I don't think this is a big deal.
There will be some coloration between the noise on $\phi$ and the noises on $x$ and $y$.  I'm not sure what this should be, but it may mean you need to introduce some off-diagonal terms in the $Q$ matrix in the KF formulation.


DSP signal filtering and number formats I'm working on a dsPIC33, using Audio Codec Board - Proto to read audio samples and for reproduction.
In order to implement a relatively simple signal processing algorithm, I'm reading audio samples and I'm supposed to process them and pass them to DAC for reproduction.
The dsPIC33 I'm using filters the signal using Q15 format, but the samples I'm getting from Codec board are signed 2's complement, 16-bit. Thus far, I managed to get the reading and send it to DAC, without any data manipulation, but now I need to filter the signal.
Is there a way to overcome this format barrier? I suppose that conversions from signed 2's complement to Q15 and vice versa are necessary.

I suppose that conversions from signed 2's complement to Q15 and vice versa are necessary.

Yes & No.  Q15 (signed) uses 2's complement, the bits are the same so no conversion is necessary. The difference between a Q15 and Q0 aren't the bit them selves but your interpretation of them. If you multiply two 16-bit numbers you get a 31-bit long result (one less than 32 because of a redundant sign bit). If you want to convert this to Q15 you take the top 16 bits, if you want Q0 you take the bottom 15 bits plus the sign bit.
The problem is that fixed point audio processing is difficult and 16-bit audio processing is extremely difficult. You need to constantly manage head room and optimize signal to noise ratio for each processing step. Fixed point filters require careful management of section ordering, section gain, section topology, noise shaping, rounding strategies, limit-cycle prevention, etc.
If you can afford it, I strongly recommend to do this in 32-bit floating point, which is much, much easier. If you can't, consider 32-bit or at least 24-bit fixed point. If that's not possible either, at the very least use double precision (fixed point) accumulators and filter state variables. If that's also not an option, you have a very rough project ahead of you.


Estimate gap between two adjacent images I have two images of a panoramic view but there is a small gap (<24 pixels) in the horizontal dimension between them. I would like to do interpolation/inpainting to fill the gap, but the exact size of the gap is unknown.

Is there a way to estimate the best value for the gap? My first thought is to detect lines/curves in the images (Canny edge detection et al) and somehow try to connect the lines at the adjacent edges, and see which value for the gap gives the smoothest continuation. However I don't known where to begin, or if there is a better/simpler solution.

There are multiple papers dealing with this issue, using similar methods (extrapolation beyond the image edges to create a virtual overlap).

Alignment and Mosaicing of Non-Overlapping Images by Yair Poleg & Shmuel Peleg

Building a Mosaic from Non-Overlapping Images by Benjamin Choi

A Novel Technique for Non-overlapping Image
Mosaicing based on Pyramid Method by Achala Pandey & Umesh C. Pati


White noise generation for different bandwidth signals in MATLAB I'm working on simulating an uplink scenario in satellite communication. The goal is to send a preamble + data to the satellite from ground station. Data symbols have a bandwidth of $15\ \rm kHz$ and the preamble has a bandwidth of $180\ \rm kHz$.
I'm simulating the transmitter, channel and receiver in MATLAB. I've a question regarding how to introduce white noise.
I read that variance of white noise depends on the bandwidth of the signal.
$$SNR = \frac{P}{NB}$$
Where $P$ is the power of signal, $N$ is PSD of white noise, $B$ is the bandwidth of the signal. In above case, the signal contains 2 parts, preamble + Data. Preamble has $180\ \rm kHz$ of bandwidth, whereas data has only $15\ \rm kHz$.
How to generate the white noise in this scenario? What should be the variance of the noise?

One technicality to be aware of is that white Gaussian noise has infinite bandwidth. What you are asking for is filtered white noise, which has a flat PSD in a given bandwidth, and zero power outside of it. It is common to call this filtered noise "white noise", but you should be aware of the difference.
With that out of the way, consider the channel between the ground station and the satellite. The noise in this channel is modeled as white Gaussian, with $\text{PSD} = N/2$. Assuming a filter with bandwidth $B$, the power of the filtered noise is, as you say, $NB$ (because the power is the integral of the PSD from minus infinity to infinity, which includes negative frequencies).
Now you need to define the receiver filter. It is very unusual, in my experience, to transmit a signal with two different bandwidths. You have two choices for the receiver design: a single filter with bandwidth 180 kHz, or two filters, one of 180 kHz and another of 15 kHz.
The single-filter design is obviously cheaper and much simpler to implement. The drawback is that it is quite suboptimal, as the data signal will be subject to much more noise than it would with a 15 kHz filter.
The two-filter design will be much more complex (I have never seen a real system implemented like this). The benefit is that it will be optimal.
In Matlab, you can generate filtered white noise of variance 1 with randn. After that, it's just a matter of scaling the noise appropriately to get the variance you want.