Phase synchronization in BPSK In BPSK for demodulation at the receiver end we do phase recovery between received wave and carrier wave. We can find the phase difference between the carrier and received waves using :

PLL
And  Hilbert transform.

Phase locked loop recovers the phase difference after some iterations and becomes constant (when it gets a lock on exact phase difference).

This is the code for finding phase difference using Hilbert, the function take two waves as parameters and returns the phase difference between them.

def phase_shift(carr_wave, rec_wave):
    assert len(carr_wave)==len(rec_wave)
    carr_comp = hilbert(carr_wave)
    rec_comp = hilbert(rec_wave)
    c = np.inner( carr_comp, np.conjugate(rec_comp) ) / math.sqrt( np.inner(carr_comp,np.conjugate(carr_comp)) * np.inner(rec_comp,np.conjugate(rec_comp)) )
    phase_diff = np.angle(c)
    phase_diff = abs(phase_diff)
    return phase_diff

We can find the phase difference using Hilbert function using whole waves(carrier and received) or phase difference at every bit( bit by bit )

My question is that:

In practice phase recovery/synchronization is done just at starting of two waves(carrier and received) or at every bit?

If i do the bit by bit phase synchronization:My demodulation result gives all zeros or ones here is why.

say my carrier wave: $S(t) = \cos( 2\pi f_c t )$ for bit 1, and $S'(t) = \cos( 2\pi f_c t + \pi )$ for bit 0

Lets say at the receiver end there is a phase difference ($\varphi$) at every bit, between received and carrier waves.

So recovered carrier looks like this $\cos(2\pi f_c t + \varphi)$ for 1 and $\cos(2\pi f_c t + \varphi + pi)$ for 0

Hence demodulation will be  $A_m(t)\cos(2\pi f_c t + \varphi)*A\cos(2\pi f_c t + \varphi) = A^2$. for bit 1

and $A_m(t)\cos(2\pi f_c t + \varphi + \pi)*A\cos(2\pi f_c t + \varphi + \pi) = A^2$ for bit 0.

Hence both bit zero and one give a positive value so how do we decide the decision boundary.

Should phase recovery be done only at the beginning of the received
wave?
If so how do we incorporate if there is a phase error with time?

Phase (or carrier) Recovery for BPSK can be done over the entire sequence using the information from every sample.  Here are common approaches to doing Carrier Recovery:

Frequency Doubling (squaring): If you square a BPSK signal (multiply it by itself) a strong tone will be created at twice your carrier frequency. The Squaring operation strips the data modulation; since the data is modulating the phase back and forth 0 to 180°, and since squaring doubles the frequency, and doubling the frequency doubles the phase, then the phase after doubling will modulate back and forth 0 to 360°, hence the modulation is removed. A PLL is used to track the new unmodulated signal at twice the carrier frequency, which is effectively a self tuning narrow band filter, and creates a copy of the signal with better noise performance. This is followed by a frequency divide by 2 and the carrier is recovered.

Cross product detectors: Regarding other approaches to carrier recovery, that are applicable to BPSK, QPSK and QAM in general, see the block diagram below showing the typical components of a carrier recovery loop:

This loop works in conjuction with a timing recovery loop, and most timing recovery loops can establish timing with a carrier offset, so assume that the timing recovery function has established proper timing and the carrier recovery is working off of the pre-decision best estimate of each symbol. For more info on timing Loops see  Location of Matched Filter and Gardner Timing Recovery for Repeated Sybmols.

Without going into deeper details of proper digital control loop design (which is necessary to understand), I will provide here details on implementations of the phase detector for use in carrier recovery.

First, importantly, the phase between two complex samples is directly proportional to the imaginary portion of the complex conjugate multiplication of those two samples. (Specifically it is proportional to the sine of the phase, but for small angles $\textrm{sin}(\theta ) \approx \theta$ ).  That is a mouthful! Luckily there is a simple and intuitive explanation that I have explained in detail in this post: How to correct the phase offset for QPSK I-Q data, but here is the direct mathematical explanation with reference to the upper plot below:

$Ae^{j(\omega t + \phi)} [Ae^{(j\omega t)}]^*= Ae^{j(\omega t + \phi)} Ae^{(-j\omega t)} = A^2e^{j\phi}= A^2cos\phi + jA^2sin\phi$

With the imaginary part equal to our output of interest: $A^2 sin\phi$

As described in the link, for two complex samples $I_1+jQ_1$ and $I_2+jQ_2$, and treating the $A^2$ as a proportionality constant, this becomes the cross product:

$A^2 sin\phi = I_2Q_1-I_1Q_2$

For small angles, phase is proportional to $I_2Q_1-I_1Q_2$, otherwise known as the cross product. Wow, it doesn't get much simpler than that!  (And for larger angles, all the way from $+\pi$ to $-\pi$ we can still use this phase detector, just be aware of the changing gain in the discriminator when designing the loop; mostly affecting acquisition time.)

See the upper plot below showing this, and first make note that we get a fairly linear phase discriminator output for small angles centered about 0° phase; further this discriminator will provide non-ambiguous phase information (although non-linear due to sine function) for angles from $+\pi$ to $-\pi$.

The lower plot shows the similar relationship in using a multiplier as a phase detector for real signals. In this case the output is proportional to the cosine of an angle, and if used in a carrier recovery loop, the loop will lock to a 90° phase offset since the 0 crossing of the discriminator occurs at 90°. In this post I am focusing on the complex implementation that would be more prevalent in an all digital design.

The question then in the specifics is what are the two signals that you use to compare phase? Here there are several possible approaches, but a common theme is you need to remove the data that is causing the phase to switch back and forth 180°- this can be done with the frequency doubling approach, comparing the doubled (squared) modulated signal with a local reference at twice the carrier, or a decision directed approach where you compare the received signal at the proper timing location to the closest decision of where that symbol should be. Since in this case it is comparing the signal to itself, the modulation effect is naturally removed.

There are also implementations that end up being a Frequency Lock Loop (so will lock in phase but to an arbitrary phase offset) that work by comparing the phase of the current symbol to the phase of the previous symbol (with data removed), or even the phase between samples within a symbol. This measures a change in phase versus a change in time, which is by definition frequency. This approach uses the imaginary of the complex conjugate method, and is often referred to as a cross product frequency discriminator.

Please refer to How to correct the phase offset for QPSK I-Q data for more details on using this concept for phase detector implementations.

This decision directed approach is basically the Costas Loop implementation; see below the implementation that works for BPSK and QPSK (adding more decision boundaries extends this to QAM). Notice that this is simply the implementation of $I_2Q_1-I_1Q_2$ described above, with one of the I, Q samples being the pre-decision sample, and the other being the post-decision sample:

And for context I show below a complete all-digital carrier recovery loop implementation for BPSK, QPSK and QAM, without the design details; the complexity is more in the Loop Filter and how to set the coefficients, other than that the block diagram shown is rather straight-forward: the accumulator followed by the two look-up tables is simply a complex NCO (numerically controlled oscillator), and the four multipliers and adders is a full complex multiplier implementation. The assumption is that earlier in the system the real signal is down-converted from RF and sampled with I and Q components (in either the analog or digital domain).  It is also feasible that the input be real (all on I and Q set to 0, and the digital multiplier in its down-conversion process by multiplying with a complex NCO will create the quadrature I and Q channels needed).


The Fourier transform of a triangular window function Windowing functions are explained on Wikipedia and they include plots of the Fourier transform for several windowing functions. For the triangle window  they provide

The above Fourier transform plot has a horizontal axis that goes from -40 bins to +40 bins, so this seems to be a discrete Fourier transform. However, I don't get a plot like that when I apply an FFT to 80 samples of a triangular window. I am further confused in that the graphic of the Fourier transform looks like a continuous plot rather than 80 discrete samples.

My questions are:
How does one create the plot for the Fourier transform of a triangular window?
Is the plot for a continuous Fourier transform or discrete Fourier transform?

Not the perfect fit, but a quick-and-dirty base to play on:

%nBin = 81;
nBin = 2049;

y=(20*log(fftshift(abs(fft(triang(81),nBin)))))';
x = linspace(-40,40,length(y));
basey = min(0,min(y));
ShadingColor = [0 0 1];
h = fill([x x(end) x(1)], [y basey basey], ShadingColor);
set(h, 'EdgeColor','none', 'FaceAlpha', 0.5);


Frequency Representation of Downsampled Signal I'm trying to follow the steps from Oppenheim for the derivation of the frequency representation of a signal which is to be down sampled by a factor of M.

It makes sense to me that:

But then the book says to make a change of variables from $r$ to $i$ and $k$ where $i$ = $0, 1, ... M-1$, and $k$ = (-$\infty,\infty$). However, I don't really understand how this change of variables comes logically. After making this change of variables it is then easy to see how the downsampled frequency representation is just a sum of scaled and shifted representations of the original sampled sequence. But I would like to be able to make the bridge better between steps and how such a change of variable was concluded.

In the final result, you want to express the spectrum $X_d(e^{j\omega})$ in terms of $X(e^{j\omega})$, the spectrum of $x[n]=x_c(nT)$. Since $X(e^{j\omega})$ is already periodic, it must be possible to represent $X_d(e^{j\omega})$ as a sum of a finite number ($M$) of shifted versions of $X(e^{j\omega})$. This is why the original infinite sum is split up into a finite sum of infinite sums, the latter being shifted versions of $X(e^{j\omega})$.

I think you shouldn't worry if you think you wouldn't have come up with that change of variables yourself. What is important is that you understand what's going on.

Furthermore, I think it's instructive to understand the derivation of the expression for $X_d(e^{j\omega})$ without introducing an auxiliary continuous-time signal $x_c(t)$:

$$\begin{align}X_d(e^{j\omega})&=\sum_{n=-\infty}^{\infty}x[Mn]e^{-jn\omega}\\&=\sum_{n=kM}x[n]e^{-jn\omega /M}\end{align}\tag{1}$$

where in the second sum we only sum over indices $n$ that are integer multiples of $M$. If we introduce a sequence $d[n]$ which equals $1$ for $n=kM$ ($k\in\mathbb{Z}$), and zero otherwise, we can rewrite $(1)$ as

$$X_d(e^{j\omega})=\sum_{n=-\infty}^{\infty}x[n]d[n]e^{-jn\omega /M}\tag{2}$$

An expression for the sequence $d[n]$ satisfying our requirements is

$$d[n]=\frac{1}{M}\sum_{l=0}^{M-1}e^{j2\pi ln/M}\tag{3}$$

Inserting $(3)$ into $(2)$ gives the final result:

$$\begin{align}X_d(e^{j\omega})&=\frac{1}{M}\sum_{n=-\infty}^{\infty}x[n]\sum_{l=0}^{M-1}e^{j2\pi ln/M}e^{-jn\omega /M}\\
&=\frac{1}{M}\sum_{l=0}^{M-1}\sum_{n=-\infty}^{\infty}x[n]e^{-jn(\omega-2\pi l)/M}\\
&=\frac{1}{M}\sum_{l=0}^{M-1}X\left(e^{j(\omega-2\pi l)/M}\right)\tag{4}\end{align}$$


Technique of sampling a high-frequency signal (Mhz) with a lower sampling rate (khz)? I currently have a .wav signal that was recorded at a 48 kHz sample rate, with a central frequency of 5.260 MHz and bandwidth of 4 kHz. I'm trying to get some information from this signal, specifically the location of the frequency peaks using matlab, which I would expect to be around 5.260 MHz. The code I'm using to generate the power spectral density is as follows:

path = '5.260MHz.wav';
f0 = 5.260e6; % frequency 5.260 MHz

% sample properties
[x, Fs] = wavread(path); % Fs = 48 kHz
alias = 48e3*floor(f0/48e3);      % starting frequency of target alias

%estimate spectrum
[psd, f] = pwelch(x, 512, [], [], Fs);

% If the frequency lies in one of the mirror frequency bands,
% we have to rotate about half the sample rate.
if (mod(f0, 48e3) > 24e3)
    f = Fs - f;
end

plot(f + alias, psd)

This script produces the following plot

The bandwidth does seem to be 4 kHz as expected, but the frequency peaks appear to be in the wrong place. I would expect them to occur between 5.256 and 5.264 MHz, or maybe 5.258 and 5.262 MHz, not between 5.276 and 5.280 MHz.

However, when I generated my own files with the same central frequency and bandwidth, I got the following when I ran the above script:

% signal properties
f0 = 5.260e6;           % frequency 5.260 MHz
sf = 4e3;               % bandwidth 4 kHz

% sample properties
fs = 24e6;              % sampling frequency 24 MHz
N = 10e6;               % number of samples

% generate random frequency modulated sinusoidal signal
fi = smooth(randn(N, 1), 11);
fi = fi / std(fi) * sf + f0;
x = cos(2 * pi * cumsum(fi) / fs);

% downsample to 48 kHz (factor 500)
x = x(1 :500: end);
fs = fs / 500;

wavwrite(x, fs, 'test_5.260MHz.wav');

This plot is exactly as expected, using the same script as the first plot. I'm wondering if I should be treating the first input different in some way to have it graphed correctly.

So far I've tried testing some other generated cosine signals using the same code as above, and they all graph correctly. However, every data signal I've tried has been off to some degree (I can post different sample plots if that would be useful). It could be an error with the data collection instrument, but it's far more likely to be a bug in my code.

You cannot capture a 5 MHz signal by sampling at 48 kHz. That doesn't mean there won't be anything in your sample, it means it will represent the signal incorrectly.

What happens if you use a too small sampling frequency is called aliasing: the original frequency shows up as another frequency altogether. The translation from the real to the sampled frequency is like a process of wrapping or folding the infinite frequency axis into your limited frequency interval $[0 ; \frac12 f_s]$, or $[-\frac12 f_s ; \frac12 f_s]$ for complex-valued signals. $\frac12 f_s$ is also called Nyquist frequency.

In your example, not just the main frequency but the total bandwidth of the signal of 800 kHz is much larger than the sampling frequency. Because of this, I wouldn't expect any clearly defined peak to be present. Since you seem to find peaks, my guess is that that is due to your very primitive approach to spectral estimation. If you use a better method, e.g. Welch's modified periodogram, which is implemented in Matlab as pwelch, I'd expect that you see a very flat (white) spectrum.

Addition based on nimrodm's comment: My initial statement has to be relativized if the bandwidth of the signal is smaller than $\frac12 f_s$. In that case, if there is additional knowledge about the true signal frequency, this knowledge can be used to resolve the ambiguity introduced by aliasing, and perfect reconstruction of the true signal is possible. See undersampling. However, it is not possible to determine the true frequency from the undersampled signal and nothing else.

A simulated example: Matlab code to generate a frequency-modulated signal with peak at 5 MHz and bandwidth 400 kHz, initially sampled at 24 MHz.

% signal properties
f0 = 5e6;               % frequency 5 MHz
sf = 0.4e6;             % bandwidth 400 kHz

% sample properties
fs = 24e6;              % sampling frequency 24 MHz
N = 10e6;               % number of samples

% generate random frequency modulated sinusoidal signal
fi = smooth(randn(N, 1), 11);
fi = fi / std(fi) * sf + f0;
x = cos(2 * pi * cumsum(fi) / fs);

Estimating the spectrum using pwelch

% estimate spectrum
[psd, f] = pwelch(x, 512, [], [], fs);
figure
plot(f, psd)
xlabel('f [Hz]')
ylabel('normalized power spectral density [1/Hz]')

gives this result:

If the signal is downsampled to 48 kHz

% downsample to 48 kHz (factor 500)
fs = fs / 500;
x = x(1 :500: end);

the estimated spectrum looks like this:

This is, as expected, basically a white (flat) spectrum, and the many small "peaks" are random fluctuations due to the fact that this spectrum was estimated based on only a 500th of the original 10 million samples. The estimate could be improved by repeating the procedure and averaging over the different results.

If the bandwidth is changed to 4 kHz

sf = 4e3;               % bandwidth 4 kHz

the results of spectal estimation on original and downsampled signal look like this:

The narrower peak fits nicely into the range 0 to 24 kHz of the downsampled spectral estimate, it just necessarily appears at the wrong place, in this case 8 kHz.

8 kHz = 5 MHz – 104 · 48 KHz is the smallest alias of 5 MHz with respect to a sampling rate of 48 kHz, the next largest is 48 kHz – 8 kHz = 40 kHz (mirroring around the Nyquist frequency of 24 kHz). All the larger aliases derive from these two by adding integer multiples of 48 kHz: 56, 88, 104, 136 kHz up to 4952, 4984, 5000, 5032 kHz and so on.

Since the spectral estimation cannot distinguish these aliases, the frequency axis of the spectrum plot becomes ambiguous: instead of representing frequency components of the signal from 0 to 24 kHz, it could as well represent frequencies from 48 to 24 kHz (backwards because of the mirroring!),

or from 48 to 72 kHz, and so on.

If we know the frequency content of the signal lies around 5 MHz, we can select the corresponding aliased frequency interval, 4992 to 5016 kHz,

and recover the correct location of the peak at 5 GHz.

This trick does not work, however, if the bandwidth of the signal is larger than $\frac12 f_s$, because the true frequency content of the signal has been aliased across several of these frequency intervals.


How is a "good" or "bad" singer defined in signal processing terms? I was watching the American and Indian idol auditions on youtube earlier in the day. I watched both successful and failed contestants and I was able to distinguish the really bad from the good singers even if I did not recognize the song or had ever heard the melody. Among the good singers, however, my limited knowledge of music prevented me from further discerning how good or bad they were. That set me thinking and I have the following questions to ask:

When we say that a singer is "good", what do we mean by it in signal processing terms?
What factors/characteristics should we consider when we say (with some authority) that the person is a "good" singer? Is it pitch, the "tune", "melody"? How much weight should these (and other parameters) have when taking such a decision?
How important a role would pitch play in this? I mean, if a man was to sing My heart will go on in perfect tune (with low pitch), would he get a score higher than a woman who sang the same song (but much less perfectly)?
When the contestants are shown singing, there is no background score playing. Yet, we are able to figure out who is singing well and who is not. So does that mean that the background score is pretty much irrelevant when it comes to deciding good/bad singing ability?
Can the above points be extended to all types of songs and genres? For instance, in most classical music (vocal), the singer is pretty much primary. The instruments etc. are present simply to aid the singer. These songs are tougher to sing (any Hindustani classical or Largo al Factotum). Songs from rock/heavy metal (We will rock you/Suicide blonde) are a little easier to sing (IMO), because the instruments are primary and vocals are secondary.
How do we (if at all) compare two songs?

I suppose that the above problem would eventually lead me to look at karaoke scoring algorithms, of which there are good open-source options (Performous, UltraStar Deluxe etc.). But before I look at these in detail, I thought it would be great to put this out and see what ideas others in the community (with a far greater knowledge of music and ability) have on this.

P.S: Reference to some papers/material on this would be great.

P.S: Some of the above examples/points are my opinions. Given what I know about music, they may be way off the mark. Please feel free to add/suggest more questions or correct opinions I may have expressed. I will be only too happy to make the necessary edits.

First order, some DSP analysis might be able to detect central pitch, pitch variation (vibrato), timbre (overtone richness), and note onset, and say whether these are grossly out-of-tune, thin, or off-time, by how much and how often.

However, second-order, good musicians play with subtile variations in pitch (equal vs. just intonation, barbershop quartet harmonization as the key changes, etc.) and in timing (jazz, blues, "swing" rhythms, etc.).  How to "score" the latter is probably more a topic for computer AI than DSP.


How to implement compressed sensing reconstruction? I am new to the field of Compressive Sensing. I'm trying to implement an example in this link. This example have described and implemented a sample tone reconstruction carefully, but unfortunately, there is no use of l1-magic toolbox to reconstruct the signal using compressive sensing minimization.

I know that l1eq_pd function in l1-magic calculates x in Ax=b for compressed sensing but when I use this function it returns an error Error using linsolve: Matrix must be positive definite.

Have anyone solved this minimization before?

Does anyone know any substitute for this toolbox?

How can I format my code in stack exchange (below code) for Matlab using <!-- language: lang-or-tag-here --> ?

Appendix 1:

My Matlab code:

%% signal initialization

n = 1/40000:1/40000:1/8;
f = sin(1394*pi*n) + sin(3266*pi*n);

%% Random Sampling

% b = Phi * f    (random sampling)
% c = Psi * f    where Psi = IDCT (sparsifying)

m = floor(rand(1,500)*length(f));
b = f(m);       % random samples of (f)
c = idct(f);    % sparsed (f)
plot(f,'b');hold
plot(m,b,'k.');title('Original Signal(f) and random sampless(b)');
axis([1,1000,-3,3])
legend('Original Signal','Randomly Sampled Signal')
figure, plot(c), axis([0,650,-10,10]);title('Sparse samples (c = IDCT(f))');

% f = Psi * c
% Psi = DCT
% f = DCT(c)
% c = IDCT(f)

%% solution 1 (x)

% Ax = b
% x = A\b

D = dct(eye(length(n),length(n)));
A = D(m,:);
% sound(f)
x = (A\b')';
b_hat = dct(x);
% sound(b_hat)

%% solution 2 (y)

% Ay = b
% y = pinv(A) * b

y = (pinv(A)*b')';
figure,plot(y),axis([0,650,-10,10])
figure,plot(dct(y)),axis([1,1000,-1.5,1.5])

%% solution 3 (s1)

s1 = l1eq_pd(y',A,A',b',5e-3,20); % L1-magic toolbox

This piece of code is a good example to understand compressed sensing and works correctly except the last line I am questioning.

I have had encountered this problem before. In using l1-magic toolbox, the first thing you need to do is to modify the file l1eq_pd.m.

To do so find all linsolve instructions and remove the third parameter, for example on line 141 change the following line from

[dv,hcond] = linsolve(H11p, w1p,opts);

to,

[dv,hcond] = linsolve(H11p, w1p);

  This is because the structure opts specify the property of the matrix
  A in solving Ax=b linear system of equations and omitting the
  structure opts in the linsolve function will allow Matlab to choose
  the appropriate solver automatically [Reference]

If the errors persists, your own code should be scrutinized.

I suggest to have a look on this examples codes.


Is it possible to estimate what the white noise floor should be, and after cross-correlation? I am simply sending a signal from a function generator, split it into two signals, and then send it into two inputs of an oscilloscope which can interface with my PC. I take the FFT of these two signal inputs and see a noise floor around 1e-7 V, but how can I know if this is the right noise floor value?

If I was to send this same signal, duplicate it again, and input it into this oscilloscope. Cross-correlate these two signals in the frequency domain (multiply with a conjugate), then I should observe a lower noise floor as the noise of the systems (function generator, oscilloscope, wires) should cancel out. But I observe a noise floor of around 1e-14 V which is really low... so I am not sure how to validate this, if it works or not.

edit: reworded my question better with right terminologies

I assume you notation "^-7V" refers to a voltage noise level of $1E-7$ volts, and then when you say it is in the frequency domain, then I assume this is the average of the FFT magnitude for all the bins that are not part of your signal (the noise floor). This average noise level is completely dependent on the number of bins in your FFT.

I do not understand clearly what your processing would be to "cross-correlate in the frequency domain", but based on your results I assume you are simply multiplying and then taking an FFT (which is not "cross-correlate in the frequency domain). If this was the case, and it was a complex conjugate multiplication then you would get the cancellation you describe, and to the degree the input signals are correlated, and the system is matched then you would likely be seeing the noise floor of your instrument.

Generally, when you perform an actual cross correlation between two noise signals, if it is the same noise at the same rms level, then the correlation result will go up in magnitude by a factor of 2. If the two noise sources are independent and identically distributed, then the correlation result will go up by $\sqrt{2}$. (Independent noise sources sum directly by their variance, correlated noise sources sum directly by their standard deviations).


Sampling fast sinusoidal signal with low sampling rate I have sampling problem that I am trying to solve analytically

I have sampling frequency, let say 1kHz.
I have number of target frequencies, let say 10kHz,100kHz,1MHz

For each  target frequency I need to find nearest frequency so that after 20 cycles of
sampling I will have 20 points that will be evenly distributed within cycle of sampled signal.  I mean that for each sample sampling point will be shifted 1/20 of period of sampled signal.  The target is to get approximation of amplitude of fast signal with slow sampling.

Do anyone have any idea how can i solve it?

In practice it sounds difficult to measure the amplitude of a 1 MHz signal (or some other very high frequency) with a 1 kHz sampler; in general, the 1 kHz analog-to-digital converter will probably rely on the input voltage being stable for some small fraction of the 1 ms sampling period; if the signal is oscillating 1000 times within that period, it's likely that the sample-and-hold will effect some smoothing/low-pass filtering, or that the conversion may be completely wrong because the signal changed so much during the conversion process.

Assuming you really do have some very fast sample-and-hold so you can record the level of your high-frequency signal "instantaneously", then your best strategy might be to sample at random times (preferably without any quantization of the sampling instant).  If you sample at a fixed rate, you run the risk of sampling at an exact multiple of the high-frequency signal's period.  So if you were taking a snapshot of a 1 MHz signal every 1 ms (i.e., every 1000 cycles) -- so your signal periods were exactly 1000:1 -- you'd see the same part of the cycle every time, and you'd have no idea what the actual waveform looked like.

However, this is statistically unlikely.  In general, when undersampling a sinusoid, you see what looks like a sinusoid with frequency $\Omega_{alias} = |\Omega_{true} - r.\Omega_{sampling}|$ rad/sec, where $r$ is the integer that makes the result of the expression less than or equal to the nyquist rate, $\Omega_{sampling}/2$ .  Except in the case where some value of $r$ makes this evaluate to zero (aliased signal appears as a constant value), you'll sample a sinusoid whose amplitude follows the true amplitude of the underlying high-frequency signal.  But of course you have no idea what the value of $r$ actually is, so you don't know the true frequency.

In general, if you want to know the peak amplitude, you're going to need to observe for a complete cycle of $\Omega_{alias}$, which could be arbitrarily small (or zero, as discussed).  But assuming it's somewhat uniformly distributed, the chance of getting an $\Omega_{alias}$ smaller than 1 Hz with 1 kHz sampling is only 1 in 500 (because we have to include -1..1 Hz, i.e. 1/500th of the 1 kHz interval), so for almost every case we'd be fine just taking the largest value over 1 sec.  You could reduce the chance of error by increasing the observation period.

Or you could have multiple, irrationally-related sampling periods, on the assumption that even if one of them happens to be close to an exact subdivision of the high-frequency period, the others will not.  There will always be a frequency that is the product of each of the individual sampling frequencies that they are all exact submultiples of, but you can make that pretty large if you have more than a few sampling systems.


Relationship between matrix rank and beamforming I always encounter the term matrix rank in papers related to beamforming. I am only familiar with the basics of beamforming (delay sum beamformer, basic capon). Can someone explain the significance of matrix rank in beamforming in a layman's term?

Typically, when doing any sort of adaptive bamforming, one needs to invert a (square) (covariance) matrix and it needs to be full rank in order to do that. Actually there are work arounds if it isn't full rank and it doesn't always require a literal inversion, like using rank one updates of QR or Cholesky decomposition.  So, that leads into what is a rank one matrix? If $\mathbf{x}$ is a column matrix, the outer product , $$ \mathbf{X} = \mathbf{x} \mathbf{x}^H
$$ is a rank one matrix. A matrix of rank $N$ is, assuming each $\mathbf{x}$ span a linear space, is a sum of $N$ rank one matrices.


Principal and Actual Values for Phase Response Is there  a way to retrieve both Principal and Actual values for a Phase Response using the fvtool (Open Filter Visualization Tool) in MatLab?

Thanks

the Principal value, $ -\pi < \Phi \le +\pi $, of phase is always

$$\begin{align}
\Phi & \triangleq \operatorname{Arg}\{H(e^{j \omega})\}   \\
 &= \operatorname{Arg} \bigg\{ \Re\{H(e^{j \omega})\} + j \,\Im\{H(e^{j \omega})\} \bigg\} \\
 &= \operatorname{atan2}\left(\Im\{H(e^{j \omega})\},\,\Re\{H(e^{j \omega})\}\right) \\
\end{align}$$

the four-quadrant arctangent:

$$\begin{align}
\operatorname{atan2}(y,x) &= \begin{cases}
\arctan(\frac y x) &\text{if } x > 0, \\
\arctan(\frac y x) + \pi &\text{if } x < 0 \text{ and } y \ge 0, \\
\arctan(\frac y x) - \pi &\text{if } x < 0 \text{ and } y < 0, \\
+\frac{\pi}{2} &\text{if } x = 0 \text{ and } y > 0, \\
-\frac{\pi}{2} &\text{if } x = 0 \text{ and } y < 0, \\
\text{undefined} &\text{if } x = 0 \text{ and } y = 0 \\
\end{cases} \\
\\
&= \begin{cases}
\arctan\left(\frac{y}{x}\right) &\text{if } x > 0, \\
\frac{\pi}{2} - \arctan\left(\frac{x}{y}\right) &\text{if } y > 0, \\
-\frac{\pi}{2} - \arctan\left(\frac{x}{y}\right) &\text{if } y < 0, \\
\arctan\left(\frac{y}{x}\right) \pm \pi &\text{if } x < 0, \\
\text{undefined} &\text{if } x = 0 \text{ and } y = 0 \\
\end{cases} \\
\end{align}$$

if, by "Actual value", you mean the unwrapped phase, if you have $\frac{N}{2}$ data points for frequency

$$ \omega[k] \triangleq \frac{2 \pi}{N}k $$

then, assuming real coefficients, define your phase at DC to be:

$$ \phi[0] = \Phi[0] = \begin{cases}
0 \quad & H(e^{j 0}) = H(1) \ge 0 \\
\pm \pi & H(e^{j 0}) = H(1) < 0 \\
\end{cases} $$

then for $ 1 \le k < \frac{N}{2} $ unwrapped phase is

$$ \phi[k] = \phi[k-1] + \Phi[k] - \Phi[k-1] + \theta[k]  $$

where

$$ \theta[k] = \begin{cases}
2 \pi \quad & \text{if } \Phi[k] < \Phi[k-1] - \pi  \\
0 & \text{if } \Phi[k-1] - \pi \le \Phi[k] \le \Phi[k-1] + \pi  \\
-2 \pi \quad & \text{if }  \Phi[k-1] +\pi < \Phi[k]   \\
\end{cases}$$

that's pretty much what MATLAB's unwrap() does.

there's a better way to do the phase difference of adjacent discrete frequencies:

the unwrapped phase for $ 1 \le k < \frac{N}{2} $ can also be expressed as

$$\begin{align}
 \phi[k] &= \phi[k-1] + \Phi[k] - \Phi[k-1] + \theta[k] \\
  &= \phi[k-1] + \Delta \phi[k] \\
\end{align}$$

where

$$\begin{align}
 \Delta \phi[k] &= \Phi[k] - \Phi[k-1] + \theta[k] \\
  &= \phi[k] - \phi[k-1] \\
  &= \arg\{H(e^{j \omega[k]})\} - \arg\{H(e^{j \omega[k-1]})\} \\
 \\
 &= \arg\left\{ \frac{H(e^{j \omega[k]})}{H(e^{j \omega[k-1]})} \right\} \\
 \\
 &= \arg\left\{ \frac{\Re\{H(e^{j \omega[k]})\} + j \Im\{H(e^{j \omega[k]})\}}{\Re\{H(e^{j \omega[k-1]})\} + j \Im\{H(e^{j \omega[k-1]})\}} \right\} \\
 \\
 &= \arg\left\{ \frac{(\Re\{H(e^{j \omega[k]})\} + j \Im\{H(e^{j \omega[k]})\})(\Re\{H(e^{j \omega[k-1]})\} - j \Im\{H(e^{j \omega[k-1]})\})}{(\Re\{H(e^{j \omega[k-1]})\} + j \Im\{H(e^{j \omega[k-1]})\})(\Re\{H(e^{j \omega[k-1]})\} - j \Im\{H(e^{j \omega[k-1]})\})} \right\} \\
 \\
 &= \arg\left\{ \frac{\Re\{H(e^{j \omega[k]})\}\Re\{H(e^{j \omega[k-1]})\} - j^2 \Im\{H(e^{j \omega[k]})\}\Im\{H(e^{j \omega[k-1]})\} + j \Im\{H(e^{j \omega[k]})\}\Re\{H(e^{j \omega[k-1]})\} - j\Re\{H(e^{j \omega[k]})\}\Im\{H(e^{j \omega[k-1]})\}}{|H(e^{j \omega[k-1]})|^2} \right\} \\
 \\
 & = \arg\left\{ \begin{matrix} \Re\{H(e^{j \omega[k]})\}\Re\{H(e^{j \omega[k-1]})\} + \Im\{H(e^{j \omega[k]})\}\Im\{H(e^{j \omega[k-1]})\} \\ + \ j  \big(\Im\{H(e^{j \omega[k]})\}\Re\{H(e^{j \omega[k-1]})\} - \Re\{H(e^{j \omega[k]})\}\Im\{H(e^{j \omega[k-1]})\} \big) \end{matrix} \right\} \\
 \\
 & = \arctan \left( \frac{\Im\{H(e^{j \omega[k]})\}\Re\{H(e^{j \omega[k-1]})\} - \Re\{H(e^{j \omega[k]})\}\Im\{H(e^{j \omega[k-1]})\}}{\Re\{H(e^{j \omega[k]})\}\Re\{H(e^{j \omega[k-1]})\} + \Im\{H(e^{j \omega[k]})\}\Im\{H(e^{j \omega[k-1]})\}} \right) \\
 \\
\end{align}$$

because $\Delta \phi[k]$ is expected to be

$$ -\frac{\pi}{2} < \Delta \phi[k] < \frac{\pi}{2} $$

even smaller in most non-pathological cases.

that might be the best way to unwrap phase given the real and imaginary parts of your frequency response function.


How to get to shoulders localization on image photo? I have a project on my university. I have to localize from video where are shoulders of person in movie.

Do have any advise how to get to this?

I thought about corner detection or some kind of shape detection. But I'm still not sure what next. We can treat video like image sequence (I wrote this, but I think is obvious)

Are there just shoulders in the video, or is there a head attached to them?

If it is the latter, then you can use vision.CascadeObjectDetector from the Computer Vision System Toolbox for Matlab.  It has the "upper body" mode, which effectively detects head and shoulders.


Delta-limiting filter Consider a filter which tracks $x_i$, clipping $|x_i - x_{i-1}| >$ some limit, say 1:

in:  0 1 2  10 0 0  10 10 10  0 0 0 0
out: 0 1 2   3 2 1   2  3  4  3 2 1 0  -- |delta| <= 1

(in Python with Numpy:  cumsum( clip( diff( x ), - limit, limit ))

Can anyone point me to a description of such filters,
either in theory or in practice ?
Or is there no theory -- for the task of smoothing concrete $x$,
just vary the limit and see what happens ?

(Added: "This tutorial was written for normal engineers, who do not have nonlinear filters for breakfast." -- F.Daum, Nonlinear filters: beyond the Kalman filter. Is this online anywhere ?)

I guess that would be a slew rate limiter. This is concept is mostly used an amplifier design as a practical constraint of the circuit. I haven't seen it applied as a digital filter. It is certainly very non-linear and it would make a poor "smoothing filter" as it's highly dependent  on the absolute amplitude. Could you shed some light on the specific application and explain what a regular low pass filter can't do the desired smoothing job?


Change in frequency on differentiation Is there any possible periodic  signal can exist(even mathematically) whose period gets change after differentiation?

No, in a conventional sense of a "periodic signal" phrase, but, if you permit me to delve into a math subtlety, differentiation can turn an aperiodic waveform to a periodic one:
$$
\frac{\mathrm{d}}{\mathrm{d}t}(a\cdot t + b\cdot \cos(\omega\cdot t)) = a-b\cdot\omega\cdot\sin(\omega\cdot t)
$$
Notwithstanding a dubious usefulness of this excursion.


Deciding an EQ's Center Frequencies I am currently making a relatively simple EQ of a music app, for which we decide to have 5 bands. Other than that, we want to have a bass boost knob, which, as I view it, is just another Peak EQ filter — so a 6 band EQ.

Are there any standards on what a 5(or 6) band EQ's frequencies should be?
What should the bass boost's frequency be?

Also, we are making only one Q knob to control the Q for all bands, simply because there isn't any space left for more knobs to be added, but we still want to provide this functionality.

Is this a good practice/design?
If not, is it at least acceptable or useful?

The audible range is about 10 octaves, and usually the center frequencies of a graphic equalizer would be distributed equally spaced on a log scale to cover that range. Common equalizers have either $30$ bands (with $1/3$ octave filters) or $10$ bands (with $1$ octave filters).

If you want $5$ bands, you could choose filters that cover approximately $2$ octaves. The five center frequencies would then be something like

\begin{matrix}
30\,\text{Hz} & 125\,\text{Hz} & 500\,\text{Hz} &2\,\text{kHz} &8\,\text{kHz}
\end{matrix}

If you have an equalizer like that I'm not sure if an additional bass boost button makes much sense. If you want $6$ frequencies, just divide the $10$ octaves into $6$ equally spaced bands (on a log scale).


Is there a means to identify/separate multiple out of phase copies of the same frequency? Given a signal that consists of three 1.xxx kHz signals added together (where xxx is never more than 0.01). Each of the three 1kHz signals is 60 degrees out of phase (first is offset by 0 degrees, next is offset by 120 degrees, and the third is offset by 180 degrees), What's the best way to detect the fact that there are three "instances" of very similar frequency present, and the phase offset of each?

I tried using FFT but it doesn't seem to have a notion of phase. I searched the internet a bit for FFT and phase but can't find any way to use FFT to extract multiple instances of a frequency separated only by phase.

I think in order to use FFT in the traditional sense, I would need to have an extreme number of bins (1000*sample_rate?).

I built a gnuradio flow graph to experiment with the concept.

The resulting FFT shows peaks at 0kHz and 1kHz. I can't think of a way to extract the original three signals out of the original.

My use case deals with input signals that differ in frequency by as little as 0.001Hz. Detecting the presence of the component signals and identifying  their precise frequency is what I need to do. Their frequencies vary ever so slightly over time.

It is not possible. Let's consider the case of two cosines, with angles $\theta = 2\pi f_0 t + \alpha$ and $\phi = 2\pi f_0 t + \beta$. Then their sum is equal to $$2\cos\left(\frac{\alpha-\beta}{2}\right)\cos\left(2\pi f_0 t + \frac{\alpha+\beta}{2} \right).$$

Note that, if you measure the amplitude and phase of this signal, you can find plausible candidates for $\alpha$ and $\beta$, but in general there will be multiple possible solutions. Note that the signal could also be zero! As you increase the number of original cosines at $f=f_0$, the ambiguity will also increase.

As a final note, the DFT does provide phase information. The GnuRadio tools, however, only plot the magnitude of the DFT, and there the phase information is lost.


Could you please confirm my method or correct me on computing the phases of my signals? I did some measurements, and I have two received signals in two different distances.

By doing the fft on the signals, I have the spectrum, then by using the angle() function in MATLAB on my fft, I computed the phase for each frequency.

I understand that the angle() function just give me the angle of each point of fft. I'm confused, is this the phase of my signal? Is this the phase relative to what signal?

As far as I know, phase of a signal, a Line-of-sight signal, after traveling a distance of L is computed as this:

-2*pi* f*(L/c)

c: is the speed of light.

The number that the angle(fft()) gives me is very different. Not even close to the above formula.

And I did all of this in MATLAB.

As Marcus Müller says, phase is relative.  Consider the Discrete Fourier Transform equation:

$$X[k] = \sum_{n=0}^{N-1}x[n]e^{-j2 \pi nk/N}$$

The DFT at element $k$ is the inner product of the signal with a complex exponential with frequency $2{\pi}k/N$ cycles/sample and an initial phase of 0.  Thus the phase of element $k$ is the phase difference of a particular frequency component of the signal relative to a particular complex exponential.

The equation you referenced is the phase shift of a single tone that travels a distance $L$.


Ultrasound Signal Processing Methods I am sending 6.6Vpp square wave across a narrow-band, 40 kHz Transmitter to a narrow-band, 40 kHz Receiver. I want to find a very accurate distance between the two (sub millimeter).

My assumption so far is that if I can get very consistent results at a fixed distance, that I can calibrate the speed of sound and an offset error time to get consistent results (and absolute, not just relative) at a bunch of different results.

Unfortunately, I've tried 3 different signal processing methods and none of them are yielding the perfect results I'm looking for. I'm sampling at 500 kHz.

Threshold method: I transmit 10-20 periods, and then look for the received waveform to pass a given threshold, then I take this sample to be the TOF. Obviously if there is noise, there is going to be some error, so its understandable why this doesn't work.
BFSK method: I transmit 50 periods at 41 kHZ, 20 periods to linearly transition to 39 kHz, and then 50 periods at 39 kHz. I receive the waveform, and use a zero detection algorithm, then detect instantaneous frequency. When I detect the frequency crossing at 40 kHz, I use this to determine TOF. This is my most accurate method. But I come across a strange error: sometimes I get all of my results with precision all within 0.5mm (500 trials), but othertimes I run this code, I get a BUNCH of errors, like 50% of the errors are greater than 5 mm.
Cross correlation method: I store a previous received signal that was taken at a known distance in memory, then I correlate my received signal to this result. My correlation only uses 4-5 periods on each side of the peak of the trial signal to reduce the computation time of the correlation.

Any recommendations on real-time, robust, sub-mm signal processing methods for ultrasound rangefinding?

Cross correlation should work.  I think the problem is the waveform that you are using.  A square wave has bad auto-correlation properties.  If it is a periodic square wave it will have multiple peaks.  It sounds like you are just using a single pulse which is better, but it will still have a gradual roll-off which is a problem.  Instead, use a Barker code, which will give you a nice sharp peak and very low side peaks.  I would BPSK modulate the barker code and then cross correlate what you receive with the "ideal" signal.


Welch’s Algorithm of Averaging Modified Periodograms I am second year Undergrad student in computer engineering . I have a homework problem where I needed to compute Power Spectral Density . I'm using MATLAB as I noticed that it has good documentation and also a function to that.
I'm going through this page to understand using dsp.SpectrumAnalyzer System object™
https://www.mathworks.com/help/dsp/ug/spectral-analysis.html
I used this object to get my results , but i failed to understand how it is doing so. I vaguely understand the concept of PSD that it is computed by Fourier transform of AutoCorr of discrete sequence as FT cannot be applied directly to the random processes.
I'm stuck and confused on this point itself

If I have "m" samples and I need to compute PSD , Is this the correct procedure as described in POINT 1

Divide m Samples into N point segments . --> this results in total of m/N segments (true ?)

Each segment which has N points is further divided into L overlapping segment having M points .

And these M points might have D overlapped points
Is my understanding correct ?

Thanks for taking time to read my question.

Understanding the Windowing Method in PSD Calculation

Please go through the link .. Someone has posted their code snippet which matches with MATLAB output. Hope that might give you better understanding


What is the relationship between the discrete time and continuous time variables? While going through Proakis's Digital Signal Processing (page 21) ,he stated that
if a continuous time signal $~x(t)~$ that has been sampled each $~T~$ seconds to produce a discrete time signal $~x(n)~$ then the relationship between the variables $t$ and $n$ is :

$$   t=nT \tag($$

Question : in the LHS we have a continuous variable whereas in the RHS we have a variable that can only take step sizes of $T$ , so clearly $t$ and $nT$ do not span the same range , then how is the formula above justified ?

It's about the relationship between the discrete time signal $x_d[n]$ and the continuous-time signal $x_c(t)$:

$$x_d[n]=x_c(nT)\tag{1}$$

So you formally replace the variable $t$ by $nT$ but this just means that you sample the continuous-time signal at sample instants $t_n=nT$. So $t=nT$ is only true for the values of $t$ that we're interested in, and these are the discrete values $t_n=nT$.


Calculate the SNR of wav file? I am doing a project in which i have a "recorded voice". Now i have applied the noise reduction and get the "filtered voice". Now i want to calculate the SNR of recorded voice and filtered voice in matlab. I do not have pure voice for reference (As in all the previous post, pure signal is taken as reference in calculating SNR). I am confused how to do it. Please help !!

When you dont have the exact voice as a reference you can only go to estimate/measure the noise. Assuming you are facing white noise, the simplest thing is to look for either time-domain or frequency-domain parts of the signal, where no voice is present. In these parts there is only noise. Hence, if you measure the average signal power in these parts, you get the noise power P(noise). Then, you measure the overall signal power, which is P(voice) + P(noise). Finally, you can have
$$ SNR=\frac{P(voice)+P(noise)-P(noise)}{P(noise)} $$


Autocorrelation of white noise Autocorrelation of white noise should have a strong peak at "0" and absolutely zero for all other $\tau$ according to this. Then why is output of this code a cone shape (with the expected of strong peak at "0" instead of being flat elsewhere except at "0")

y = ones(1,1000);
z = awgn(y,1);
[m l] = xcorr(z);
plot(l,m);

Problem 1:

As pointed by learner, awgn adds noise to the sequence y. So you are not plotting the autocorrelation of white gaussian noise, but the autocorrelation of white gaussian noise plus a constant.

Use zeros instead of ones ; or analyze the difference $z - y$ ; or since $y$ is a constant signal, remove the mean of $z$.

Problem 2:

Once point 1 is solved, and since your sequence is of finite length, you will be plotting the autocorrelation of white noise multiplied by a square window. What you will see is thus the autocorrelation of white noise convolved by the autocorrelation of a square window. The autocorrelation of a square window has a triangular shape.

Another way to look at it: the further you move from 0, the less data is present in your input vector. For example, in the sequence:

$-1, 4, 5, 6, -2$

There are 4 pairs of samples distant by a lag of 1, 3 pairs of samples distant by a lag of 2 (-1 and 5 ; 4 and 6 ; 5 and -2), and only two pairs of samples distant by a lag of 3 (-1 and 6, 4 and -2). Thus, there won't be as much data to estimate correctly the autocorrelation function for the larger values of the lag.

You can artificially compensate for the windowing by providing the 'unbiased' argument to xcorr.

Problem 3:

Keep in mind that the theory is about infinite length sequences (or for finite length sequences, about the expected value of the result).


SC-FDMA vs single carrier signal PAPR Why is the PAPR of a single user SC-FDMA signal different than the PAPR of a  conventional single carrier signal if they both use the same modulation? Aren't they both single carrier signals?

Consider some actual use cases for further insight:

Unmodulated Carrier: The peak to average ratio is 3 dB given the relationship between the rms and peak level of a sine wave.

Unfiltered Single-Carrier BPSK or QPSK: This would be the modulation of rectangular shaped symbols, so would have a 3 dB peak to average ratio just as in the case of an unmodulated carrier. Great for peak-average ratio but terrible for bandwidth constraint. Given we have to share bandwidth efficiently (in terrestrial communications) we don't have this luxury. GPS is an example where the signal can be essentially unfiltered and we have the benefit of power efficiency.

Root-Raised Cosine BPSK or QPSK: This is what is traditionally done  but the tighter the filtering (smaller $\alpha$ in the filter design) results in more overshoot of the symbol transitions and slower trajectories that may go through the origin resulting in a higher peak-average ratios.

pi/4 QPSK: This modulation changes phase by $\pi/4$ between symbols thus avoiding having the slower transitions from filtering from ever going through the origin (Meaning minimizing the lowest instantaneous power in the signal envelope, increasing the average relative to the peak). This modulation is therefore considered for lower power (battery operated) applications since we are able to drive the transmit power amplifier further into saturation meaning higher efficiency.

GMSK: Gaussian Minimum Shift Keying and other constant-envelope modulations have a 0 dB peak-average ratio and offer moderate bandwidth efficiency so are an excellent choice for power efficient (battery operated) applications.

OFDM: To the extent each bin in the FFT representation of an OFDM symbol is independent and in consideration of each other randomly generated, the resulting signal as the summation of multiple independent identically distributed random processes will be well approximated by a Gaussian distribution given the central limit theorem. A Gaussian distribution has a relatively high peak to average ratio, and we need to decide on a significance to establish what is the "peak" (usually on a distortion criteria from clipping anything beyond the established peak). In SC-OFDM we do have some control over the association of the signals in each bin, allowing for lower peak-average ratios than MC-OFDM where each bin can represent a different data source.


Conflict with the properties of convolution excuse me in advance for the lenghty question, consider the following signals:

$x_1[n]$, $x_2[n]$, $y[n]$ and

$z_1[n] = x_1[n] * y[n]$

$z_2[n] = x_2[n] * y[n]$

In addition, we know that:

$x_2[n] = x_1[n] - x_1[n-3]$

Thus, we can write:

$z_2[n] = (x_1[n] - x_1[n-3])*y[n]$

$z_2[n] = (x_1[n]*y[n]) - (x_1[n-3]*y[n])$

$z_2[n] = z_1[n] - (x_1[n-3]*y[n])$

Now for the question; the solution to the exercise I'm working on states that:

$z_2[n] = z_1[n] - z_1[n-3]$

Now this would mean that the following:

$z_1[n-3] = (x_1[n-3]*y[n])$ is true.

Isn't $z_1[n-3] = x_1[n-3]*y[n-3]$ ?

What makes me hesitate is that this development seems to support the given answer:

$x_1[n-3] = x_1[n]*\delta[n-3]$

$x_1[n-3] * y[n] = (x_1[n]*\delta[n-3])*y[n]$

$x_1[n-3] * y[n] = (x_1[n]*y[n])*\delta[n-3]$

$x_1[n-3] * y[n] = z_1[n]*\delta[n-3] = z_1[n-3]$

Which answer is correct and why is there a conflict with both reasonings? Thank you!

I believe that the problem you have with this exercise is mainly caused by bad notation. Writing convolution as

$$z[n]=x[n]\star y[n]\tag{1}$$

(which I think is bad notation) it is tempting to believe that

$$z[n-1]\stackrel{?}{=}x[n-1]\star y[n-1]\tag{2}$$

which is not the case.

If you take the definition of (discrete) convolution

$$z[n]=\sum_kx[k]y[n-k]\tag{3}$$

then it's obvious that

$$z[n-1]=\sum_kx[k]y[n-1-k]=\sum_kx[n-1-k]y[k]\tag{4}$$

i.e., only one of the two signals needs to be delayed by one sample, not both.

A better notation is

$$z[n]=(x\star y)[n]\tag{5}$$

Now you can define $x_1[n]=x[n-1]$ and $y_1[n]=y[n-1]$ and you have

$$z[n-1]=(x\star y_1)[n]=(x_1\star y)[n]\tag{5}$$

If you view $z[n]$ as the output of a system with impulse response $y[n]$ and input $x[n]$, then it should also be obvious that delaying $z[n]$ by one sample is equivalently achieved either by delaying the input signal, or by adding a delay to the impulse response, but not by doing both.


time scaling and shifting of cosine in Fourier transform I've met some problems when calculating the Fourier transform of $\cos(at+b)$. I want to use the shifting and scaling properties to solve this problem. First, when I look up in the book and some online resources, I saw a post about the Fourier transform of the signal $x(3t-6)$, where it tells me that $x(3t-6)$ is scaling $x(t)$ by $3$ and shift it by $2$, which I will get the answer that
$$F\{x(3t-6)\} = \frac13 X(j\omega/3)exp(-j2\omega)$$
This indicates that the Fourier transform of $x(at+b)$ is $$\frac{1}{|a|}X(j\omega/a)\cdot \exp(j\omega b/a)$$

Thus, I apply the conclusion to the Fourier transform of $cos(at+b)$ is

$$\frac{1}{|a|}\exp(j\omega b/a)\,\pi\, [\delta(\omega/a−1)+\delta(\omega/a+1)]$$

How do I proceed to the right answer which is

$$\pi [\exp(jb)\delta(\omega-b)+\exp(-jb)\delta(\omega+b)]$$ ? I appreciate your help.

(I'm new here so I haven't figure out how to type the fractions or exponentials, sorry for the bad typing)

Using the shifting and scaling properties of the Fourier transform is a rather complicated way of computing the Fourier transform of $\cos(at+b)$. A more straightforward way is to realize that

$$\cos(at+b)=\frac12\left[e^{jat}e^{jb}+e^{-jat}e^{-jb}\right]\tag{1}$$

From the Fourier transform pair

$$e^{jat}\Longleftrightarrow 2\pi\delta(\omega-a)\tag{2}$$

we immediately obtain

$$\cos(at+b)\Longleftrightarrow \pi\left[e^{jb}\delta(\omega-a)+e^{-jb}\delta(\omega+a)\right]\tag{3}$$

But having said that, your result is correct. You just need to simplify it appropriately. In order to do so, you need to know two things. First,

$$\delta\left(\frac{\omega}{a}\right)=|a|\delta(\omega)\tag{4}$$

and, second,

$$f(\omega)\delta(\omega-\omega_0)=f(\omega_0)\delta(\omega-\omega_0)\tag{5}$$

for any $f(\omega)$ that is continuous at $\omega_0$.

From $(4)$ you get

$$\delta(\omega/a-1)=|a|\delta(\omega-a)\tag{6}$$

And from $(5)$ you get

$$e^{j\omega b/a}\delta(\omega-a)=e^{jb}\delta(\omega-a)\tag{7}$$

Applying similar transformations to the second term in your result leads to the form given in $(3)$.


How to distinguish signal and noise eigenvalues of autocorrelation matrix for MUSIC algorithm? MUSIC (multiple signal classification) algorithm is introduced in Schmidt, Ralph. "Multiple emitter location and signal parameter estimation." IEEE transactions on antennas and propagation 34.3 (1986): 276-280.

An important part of this algorithm is to distinguish eigenvalues froma sample autocorrelation matrix ($R_{N}=\frac{1}{N}\sum{XX^{H}}$) into two groups: noise eigenvalues and signal eigenvalues.

Eigenvalues are plotted as below. A few nontrivial signal eigenvalues and lots of trivial (zero or close to zero) noise eigenvalues:

The figure above is in a paper Quinlan, Angela et al. “Model Order Selection for Short Data: An Exponential Fitting Test (EFT).” EURASIP J. Adv. Sig. Proc. 2007 (2007): n. pag.. The paper is telling how to distinguish signal eigenvalues and noise eigenvalues. The authors use exponential fitting test (EFT), which is not that simple, but seems effective.

On the other hand, I've found two implementations of MUSIC algorithm in Python and MATLAB. Interestingly, they've used very simple approach: calculate a difference (i.e., a slope) between adjacent eigenvalues and find an index of the most largest difference. Then, noise eigenvalues are with indices larger than the index of the largest difference.

Intuitively, EFT can coincide with the largest difference approach. However, it is just an intuition, and does not include an analysis or reasoning.

Is it okay to use the largest difference approach to distinguish signal eigenvalues and noise eigenvalues? If so or if not, what is a good explanation of it?

This problem is defined as estimating the dimension of a model. I would recommend use something like  Bayesian Information criteria (BIC) or Minimum description length (MDL).


How to change volume of a PCM 16 bit signed audio? I know I can multiply samples, then clip but perceived volume is non-linear for humans.

Can you please help with a formula.

Changing the volume of an audio signal must be done by applying a gain (multiplication) - and optionally clipping if your system has a limited dynamic range. This is as simple as that. Applying a non-linear function to an audio signal will cause distortion and add harmonics, and you don't want this to happen - you want to modify the loudness of the signal, not its timbre. [To be fair, there are non-linear processings designed to change the perceived loudness of the signal without affecting the timbre, within a given dynamic range constraint (eg. multiband compression), but it doesn't look like this is what you need.]

Where non-linearity and fancy response curves come to play is when designing a user interface - when deciding on the relationship between the position of the control (knob or slider, whether on a GUI or as physical hardware) and the gain applied to the signal. This is where perception matters, because the users will expect a mapping between the position of the slider and their perception of loudness. Please note that even if the relationship between the position of the volume control and the gain applied to the signal is non-linear, the process of applying the gain to the signal is linear, and non-linearity would be unwanted there!

When it comes to physical volume controls, eg in hifi systems or personal audio players, the relationship between the knob position and the attenuation is closer to an exponential curve, though its shape has been tweaked and is constrained by the manufacturing process - sometimes it's just two or three linear segments. You can find those curves in the datasheets from manufacturers ("A" taper). Mixing console faders usually have their response compressed so that the upper half of their travel covers the useful range of -20 dB..+6dB.

In the software world - at least for music production - it is most common to have volume/gain knobs calibrated in dB. For example, if you have a 100 pixels long volume slider graduated from -48dB to +6dB, the gain applied to the signal would be $10^{\frac{-48 + 54 \frac{x}{100}}{20}}$.


Synchronization Accuracy Required For Communications One of the biggest problems in phase modulated communication is discovering when your signal has arrived. High accuracy in this area is important, as even a quarter wavelength error can cause the estimated data bit to be completely wrong. One solution to this problem is to use a preamble code, which is a waveform whose autocorrelation has a very tight mainlobe and low sidelobes.

My question is, are preamble codes actually used in practice? If so, which ones? How is it that such a high accuracy of quarter wavelength is achieved in the real world using a preamble code? If this bound isn't correct, what bound do we really need to achieve synchronization using say BPSK?

Not only are you dealing with the uncertain delay or phase but also importantly there will be frequency offsets; inevitably your transmitter will not be at the exact frequency of your receiver. Carrier and symbol timing recovery loops are common receiver implementations. A preamble code is not necessarily required as such recovery and synchronization schemes can be done using the data directly, but having known sequences can aid the recovery to be faster and operate in a lower SNR condition. For a example of a carrier recovery implementation that is suitable for BPSK/QPSK/QAM see

Recovering signal for psk

For timing recovery the Gardner Loop and Mueller &Mueller synchronizer are commonly used. I posted some information on the Gardner Loop here:

Gardner Loop


Most appropriate spectral analysis method for a recording of a saxophone note? I'm currently doing an investigation where I am attempting to resonant frequencies of an Alto Saxophone at various different notes. I have taken audio recordings of the notes with a sample rate of $44100\textrm{ Hz}$, and exported them into CSV files. Now, I am attempting to analyze this data using the signal processing plugin for MATLAB, but I'm feeling a little bit overwhelmed with all the different options, as I'm not sure which spectral analysis method is most suitable for the my task.

In MATLAB, I'm seeing a number of spectral analysis methods including FFT, Burg Method, Welch Method, Multitaper Method, and more - but I'm not sure of their specific applications.

Again, my goal is to clearly identify the major spectral peaks (formants) which correspond to the resonant frequencies of the Saxophone bore. The data was taken with a condenser microphone at a sample rate of $44100\textrm{ Hz}$, with the Saxophone bore exactly 1 meter from the microphone. Ultimately, I'm looking to create a plot similar to the image I have attached.

Please excuse my ignorance on this subject, this is my first time using MATLAB and I am not well-versed in digital signal processing. If anyone could help me with this, I would be extremely grateful.

The plot OP shared looks similar to log magnitude spectrum(black line) plotted with its envelope(red line). so from what OP says, OP need to plot log magnitude spectrum of recorded audio/signal. One thing to remember for this analysis is that same note should be played for entire duration of the signal to be analysed. if not, one should use proper windowing methods to ensure that signal remains stationary(in same note) for the entire duration of window under analysis. I am assuming that a note corresponds to a few particular frequencies.

that is the way if you want to plot/analyse a graph similar to the one shared. are you also looking for other methods for determining the peaks corresponding to a particular note at a time?


DSP programming in C/C++ I have exposure to MATLAB and have some basic exposure to signal processing and plotting in MATLAB... but I am more comfortable with C/C++ than MATLAB because I have been working in Image processing with OpenCV... I was just curious: Is there any DSP library or platform (as OpenCV is for image processing) for doing signal processing in C/C++?

Octave is a multi-platform open source math and matrix toolkit. It has a command line interpreter aimed to be very similar to MATLAB, but there is also a C++ API available for use. Since you refer to signal processing in contrast to image processing, I assume you mean audio processing, so you might need to look into the "signal" and "audio" packages in octave (which provide e.g. wavread). Keep in mind however, that the API is primarily designed for internal use in the octave CLI, and that it is considered not stable (i.e. the bindings might change). Also, since the library is not primarily meant to be linked with, the documentation is poor to say the least.

Alternatively, there is also GNU Scientific Library, which is not primarily aimed at signal processing, but provides complex numbers, matrices, fft and so on. On Windows, you would have to compile that yourself using Cygwin.


Frequency spectrum of signal - Matlab Here is the code I use to plot a function in frequency domain in Matlab:

    dt = 1/10000; % sampling rate
    et = 0.1; % end of the interval
    t = 0:dt:et; % sampling range
    y = 2+sin(2.*pi.*50.*t)+18.*sin(2.*pi.*90.*t)+6.*sin(2.*pi.*180.*t); %     sample the signal
    subplot(2,1,1); % first of two plots
    plot(t,y); grid on % plot with grid
    xlabel('Time (s)'); % time expressed in seconds
    ylabel('Amplitude'); % amplitude as function of time
    Y = fft(y); % compute Fourier transform
    n = size(y,2)/2; % 2nd half are complex conjugates
    amp_spec = abs(Y)/n; % absolute value and normalize
    subplot(2,1,2); % second of two plots
    freq = (0:100)/(2*n*dt); % abscissa viewing window
    stem(freq,amp_spec(1:101)); grid on % plot amplitude spectrum
    xlabel('Frequency (Hz)'); % 1 Herz = number of cycles/second
    ylabel('Amplitude'); % amplitude as function of frequency

The problem is, when I zoom in my graph I don't see peaks exactly on 50Hz, 90Hz and 180Hz.

What did I do wrong in my code?

Discrete Fourier transform doesn't actually decompose a signal into sum of sinusoids it was composed from. DFT projects signal onto a specific set of discrete sinusoids known as Fourier basis. Sinusoids used to compose the signal may or may not coincide with the DFT basis. If they are not, they are being decomposed into a number of DFT basis elements.

There's another way of how to look at it. DFT deals with periodic functions. If you plot 3-4 periods of your signal interval, you'll see it doesn't look like a sum of pure sinusoids and has discontinuities - if sinusoids it was composed from don't coincide with a Fourier basis.


Extracting fundamental amplitude/phase from only half a period of a pseudo-sinusoïd Before you read my problem, I've summed up my experiments in a notebook available here.

So, I have a few noisy samples (~50) of a pseudo sinusoidal signal for which I use a traditional discret Fourier transform in order to extract the amplitude and phase of the fundamental (real signal is much noisier).

fft[1]=5.0, 60.0°

Now my problem is that i sometime have only half the data available, ie only half a peridod of that sinusoid.
But I know I almost exactly have sample of half that period, and I know the signal is almost a sinusoid.

How can I proceed to extract the fundamental characteristics anyway?

I've noticed that zero-padding the missing half and doubling the amplitude works quite well.

fft[1]=2.5, 59.1°

But in my case the sinusoid is not DC-centered on zero, and in that case zero-padding doesn't work.

fft[1]=10.6, -82.5°

I would need to pad with the DC of my signal, which is unknown. I've noticed that if you have half a period, the average of first and last samples should give the DC, but this won't be reliable with noisy data.

fft[1]=2.5, 59.1°

Thanks in advance for your insights!

If you know the frequency (and it seems that you do), simply construct a basis set of vectors from

1) DC (all ones)

2) Half a Sine wave

3) Half a Cosine wave

These form an orthogonal basis just like a DFT.  In fact you can think of it as a "one half bin and DC bin" specialized DFT.

Discard the DC reading (unless you care) and use the Sine and Cosine coefficients just like a regular DFT bin from which you can calculate magnitude (which gives you amplitude with a multiplication) and the phase (using atan2).

This will give you a best fit answer.  No window function and no zero padding.

Small correction (sorry):  The basis is not orthogonal.  You will need to solve a system of equations.

$$ x[n] = A \cos( \omega n ) + B \sin( \omega n ) + C $$

As a vector equation:

$$ \vec x = A \vec c + B \vec s + C \vec u $$

With some dot products, you get the system of equations with $A,B,C$ as the unknowns to be solved for.

$$ \vec x \cdot \vec c = A \vec c\cdot \vec c  + B \vec s\cdot \vec c  + C \vec u \cdot \vec c  $$

$$ \vec x \cdot \vec s = A \vec c\cdot \vec s  + B \vec s\cdot \vec s  + C \vec u \cdot \vec s  $$

$$ \vec x \cdot \vec u = A \vec c\cdot \vec u  + B \vec s\cdot \vec u  + C \vec u \cdot \vec u  $$

Note:

$$ \vec c \cdot \vec u = \vec u \cdot \vec c = 0 $$

$$ \vec c \cdot \vec s = \vec s \cdot \vec c = 0 $$

$$ \vec u \cdot \vec u = N $$

Where $N$ is your sample count.  Thus the equations simplify to:

$$ \vec x \cdot \vec c = A \vec c\cdot \vec c  $$

$$ \vec x \cdot \vec s = B \vec s\cdot \vec s  + C \vec u \cdot \vec s  $$

$$ \vec x \cdot \vec u = B \vec s\cdot \vec u  + C N $$

In order to get the best results with noisy data it is best to build the basis vectors over the full range of data you have and then solve the full set of equations.  A little more work, but likely worth it.

Okay, for the OP and maybe Ben.  Here is a sample program for the best fit of a partial cycle.  In this case 0.6 of a cycle.

import numpy as np
import matplotlib.pyplot as plt

#==========================================================
def main():

#---- Set Parameters, Base Frequency and Phase

        f_s = 50.0
        p   = 60.0 * np.pi / 180.0

#---- Set Domain Time and Angle for partial cycle

        t = np.arange( 0.0, 1.0, 1.0 / f_s )
        a = 2.0 * np.pi * t * 0.6

#---- Construct the Signal

        x = 10.0 * np.cos(        a + p ) \
          +  0.2 * np.sin(  2.0 * a )     \
          +  0.2 * np.sin( 10.0 * a )

#---- Build Basis Vectors

        c = np.cos( a )
        s = np.sin( a )
        u = np.ones( a.shape )

#---- Do Basis Dot Products

        cc = c.dot( c )
        ss = s.dot( s )
        uu = u.dot( u )
        cs = c.dot( s )
        cu = c.dot( u )
        su = s.dot( u )

        print cs, ss, uu, cs, cu, su

#---- Do Signal Dot Products

        xc = x.dot( c )
        xs = x.dot( s )
        xu = x.dot( u )

        print xc, xs, xu

#---- Build Coefficient Matrix

        M = np.array( [[ cc, cs, cu ],  \
                       [ cs, ss, su ],  \
                       [ cu, su, uu ] ] )

        print M

#---- Build Values Matrix

        V = np.array( [xc, xs, xu] )

        print V

#---- Solve for the best fit "Bin Value"

        R = np.linalg.solve( M, V )

        print R

#---- Calculate the Best Fit parameters

        A = R[0]
        B = R[1]
        C = R[2]

        D   = np.sqrt( A*A + B*B )
        phi = np.arctan2( -B, A )

        print D, phi, p

#---- Construct the Best Fit signal

        y = D * np.cos( a + phi ) + C

#---- Show the Results

        plt.plot( t, x )
        plt.plot( t, y )

        plt.show()

#==========================================================
main()

And these are the results:

These are the values:

Underlying target:

    10.0
p =  1.0471975512

Found best fit:

D   = 10.1939671316
phi =  1.04122684844

Changing the parameters a little bit for a smaller portion (0.3) of the cycle and moving it up by 20.

        x = 10.0 * np.cos(        a + p ) \
          +  0.2 * np.sin(  2.0 * a )     \
          +  2.0 * np.sin( 10.0 * a )     \
          + 20.0

Results:

Notice this isn't as accurate because of the smaller cycle portion.

D   = 10.4409622651
phi =  0.943824736628

However, if the signal is a pure tone:

        x = 10.0 * np.cos( a + p ) + 20.0

The results are still exact on less than a third of the cycle:

D   = 10.0
phi =  1.0471975512


Wiener Filter Additive Noise Uncorrelated I have been trying to understand for a while now why a non-causal Wiener filter that has a frequency response of the type

H(jω)=S_yx(ω)/S_xx(ω)

where y(n) = x(n) + n(n), x(n) the input signal, n(n) the noise and both uncorrelated,

can be simplified to

H(jw) = S_xx(jw) / ( S_xx (jw) + S_nn (jw) )

I understand the denominator but why is it that when x and n are uncorrelated

S_yx(jw) = S_xx(jw)

i.e. I don't understand Step 3 in this question but I don't have enough reputation to ask for a clarification/add a comment :(

Let me change the notation slightly to make things clear. Let $d[n]$ be the desired signal, $x[n]$ is the noisy input signal, and $v[n]$ is the noise. Furthermore, we assume that

$$x[n]=d[n]+v[n]\tag{1}$$

and that the desired signal $d[n]$ and the noise $v[n]$ are uncorrelated and at least one of them has zero mean, which implies that their cross-correlation equals zero:

$$R_{dv}[n]=E\{d^*[k]v[k+n]\}=0\tag{2}$$

From $(2)$ it follows that the cross-correlation between the desired signal $d[n]$ and the input signal $x[n]$ is given by

$$\begin{align}R_{dx}[n]&=E\{d^*[k]x[k+n]\}\\&=E\{d^*[k]\left(d[k+n]+v[k+n]\right)\}\\&=E\{d^*[k]d[k+n]\}+E\{d^*[k]v[k+n]\}\\&=E\{d^*[k]d[k+n]\}=R_{d}[n]\tag{3}\end{align}$$

From $(3)$ it follows that the cross-power spectral density of $d[n]$ and $x[n]$ equals the power spectral density of $d[n]$:

$$S_{dx}(\omega)=S_{d}(\omega)\tag{4}$$

With the given assumptions, the frequency response of the non-causal Wiener filter can be written as

$$H(\omega)=\frac{S_{dx}(\omega)}{S_x(\omega)}=\frac{S_d(\omega)}{S_d(\omega)+S_v(\omega)}\tag{5}$$


Pitch detection using Cepstrum in matlab What i am trying to do is to detect the pitch of a given sound signal, so for that purpose, i wanted to find cepstrum of a given signal.

This is my code in MATLAB:

[x fs] = audioread('a.wav');
N = length(x);
n = 0:N-1;

figure, stem(n, x);

X = fft(x, N);
length(X)

figure, stem(n, X);

Y = X(1:6400);
N = length(Y);
n = 0:N-1;

figure, stem(n, Y);

Z = log(abs(Y));
figure, stem(n, Z);

z = ifft(Z, N);
figure, stem(n, z);

Now, all i did was that i found fft of audio signal and then i cut off the second half of it because it is basically the same thing i have in the first half, then, i just followed the standard procedure for cepstum, i moved to logarimic representation and then i found inverse fourier of it. Based on what i saw on the internet, i expected evenly spaced peaks since that represents existence of a definite pitch in the signal, however, this is what i get:

So basically i have no idea what can i do with this, since this what i have is speech signal, i am not sure why there's no any evenly spaced peaks. When i zoomed in a bit, this is what i got

:

This also doesn't help at all, i am having trouble with understanding the way that cepstrum enables us to determine pitch frequency so basically, even if i had something what would seem to me as a proper cepstrum, i doubt that i would be able to determine pitch frequency from it, so any help is highly appreciated!

I wrote a very very basic cepstrum pitch track in the past, well we can try using a pure senoidal signal to show how its works... let me see

For Cepstrum I always have used this steps:

Apply hamming windows in the signal
Apply FFT
Get magnitude
Convert to log scale
Apply IFFT

The equation for cepstrum:

IFFT(log(abs(FFT(s))))

But you can use FFT or IFFT, take a look:

IFFT(log(abs(FFT(s)))) == real(FFT(log(abs(FFT(s)))))

The difference is the scale representation, if do you end using FFT you need extract just the real information, for both above equations you will get the same shape:

For IFFT(log(abs(FFT(s)))):

For real(FFT(log(abs(FFT(s))))):

This is a cepstrum shape example from 4096 points sine in 440hz sampled at 44100hz

Now try get the maximum from the half cepstrum data to find the periodicity...

see my extremely basic code of how to get a sine frequency, just to exemplify to you:

Fs=44100;

%creating one simple signal

f=440; %we need find this
frame = sin(2*pi*f/Fs*(1:2048*2))';

win = hamming(length(frame));

%samples multplied by hamming window
windowedSignal = frame.*win;

fftResult=log(abs(fft(windowedSignal)));

ceps=real(ifft(fftResult));

nceps=length(ceps)

%find the peaks in ceps

peaks = zeros(nceps,1);

k=3;

while(k <= nceps/2 - 1)
   y1 = ceps(k - 1);
   y2 = ceps(k);
   y3 = ceps(k + 1);
   if (y2 > y1 && y2 >= y3)
      peaks(k)=ceps(k);
   end
k=k+1;
end

%get the maximum ...
[maxivalue, maxi]=max(peaks);

result = Fs/(maxi+1)

my result was close 441hz for better results maybe do you need apply parabolic interpolation ....


How to measure distance between two 3d objects using cameras I need some help regarding measuring distances between two balls and their position relative to surface(coordinates).
Dimensions of balls and surface are known.
Do i need just one camera to do this(measuring with ray diagrams and trigonometry)or more?

I hope this is the right thread,I tried on physics and photography before and it wasn't correct thread

Since nobody has answered this for you yet, I'll give it a shot.

You have to understand how to do ray tracing math, which means vectors.  You also have to be able to calibrate your camera so you know for each pixel, what is the direction of the the ray coming from the focal point.

It makes a big difference if the objects are in motion or are still.  If they are in motion you will need to deal with blurry edges.

Since the screen is two dimensional, you will only be able to make two dimensional direction measurements.  However, with something like a ball with a known diameter, you can calculate how far it is from the focal point using trigonometry.

Having two cameras on the scene can give you much more accurate depth measurements.  Again, this requires knowing the rays and doing trigonometry.

If your surface is flat, you can define a mapping between your pixel locations and a coordinate system on the surface.

Many years ago, I wrote a program using a webcam, where when I identified a sheet of notebook paper on the floor by clicking the four corners with my mouse, the program calculated the location of the camera relative to the paper within a tenth of an inch.

Once you have the objects located in 3D space, finding the distance between them is merely an application of the distance formula, aka Pythagorean Theorem.

Hope this helps.

Ced

=============================================

Followup (too long for an add on comment above):

@M.M., I've never used OpenCV so I searched on "lens distortion matrix openCV" and found this:

Camera calibration With OpenCV

This approach is different than what I did, though somewhat similar.

Suppose $(m,n)$ is the pixel location and $(x,y,1)$ is the direction vector where the z-axis is straight out from the lens.  I used calibration points, and a standard best fit approach to find the coefficients of:

$$ x = A_x m^2 + B_x n^2 + C_x m n + D_x m + E_x n + F_x $$
$$ y = A_y m^2 + B_y n^2 + C_y m n + D_y m + E_y n + F_y $$

This required finding the inverse of a 6x6 matrix to solve.  I wrote the program in VisualBasic3 (I think) on an i386, to give you an idea of when.

Ced

=============================================

More detail:

Spheres are easy because they look like a solid circle from any angle.

Find the center of the circle in pixel coordinates.  A centroid calculation is a brute force way to do this, but there are more efficient ways.

Translate the pixel coordinates of the center to a ray emanating from your lens.

Then process several points around the edge.

For each point:

1) Translate to a ray

2) Calculate the angle to the center ray.

3) Figure the distance to the center of the sphere from the lens

Average the distances (they should ideally be identical) to get a good read.

Calculate the location in 3D coordinates of the center of the sphere from the center ray direction and the calculated distance.

Calculate the distances between the spheres using the distance formula.

Angle calculation formula for two vectors $\vec a$ and $\vec b$

$$ \theta = cos^{-1} \left( \frac{ \vec a \cdot \vec b }{ \| \vec a \| \| \vec b \|  } \right) $$

Distance to sphere where $r$ is the radius

$$ d = \frac{ r }{ \sin \theta } $$

The center of the sphere:

$$ \vec c_k  = d \frac{ \vec a_k }{ \| \vec a_k \| } = ( x_k, y_k, z_k ) $$

The distance between spheres:

$$ d_{k,l} = \sqrt{ ( x_k - x_l )^2 + ( y_k - y_l )^2 + ( y_k - y_l )^2 } $$

Your surface will require a different approach, but since you haven't provided any details, I can't help you there.

Ced


What is the effect of aliasing on the magnitude of the autocorrelation? I've a question about the effect of aliasing on the magnitude of autocorrelations. From a simulation in MATLAB, I don't see any effect of aliasing or any need to anti-alias filter when I take the magnitude of the autocorrelation. Which means I can undersample the data and then take the autocorrelation. There is a paper "Effects of Aliasing on Spectral Moment Estimates Derived from the Complete Autocorrelation Function" which says something like what I claim. Would anybody please let me know if I've made a mistake?

Decimating before calculating the autocorrelation, in the presence of noise, is inferior to calculating the autocorrelation using the full dataset. Assume that the signal of interest is embedded in white noise. The vector $x[n], n = 0, 1, ..., N-1$ consists of samples from a discrete random process. The autocorrelation function of the vector $x[n]$ is:

$$
A_x[k] = \frac{1}{N-k} \sum_{i=0}^{N-1-k} x[i]x[i+k]
$$

That is, $k$ is the lag used for the autocorrelation calculation. In your proposed scenario, you are decimating the autocorrelation function output by a factor $D$ (i.e. you are only calculating the function for lags $0, D, 2D, ...$) and comparing that result to the autocorrelation function of $x[n]$ decimated by the same factor $D$. Let $x_d[n]$ be the decimated sequence; its autocorrelation function is:

$$
A_{x_d}[k] = \frac{D}{N-k} \sum_{i=0}^{\frac{N-1-k}{D}} x[iD]x[(i+k)D]
$$

(for simplicity here, I have assumed that $D$ is a factor of $N$ in the above equation)

Your inquiry can be written as:

$$
A_x[kD] \stackrel{?}{\approx} A_{x_d}[k]
$$

$$
\frac{1}{N-kD} \sum_{i=0}^{N-1-kD} x[i]x[i+kD] \stackrel{?}{\approx} \frac{D}{N-k} \sum_{i=0}^{\frac{N-1-k}{D}} x[iD]x[(i+k)D]
$$

Looking at this qualitatively, the summation on the left hand side has more terms than its counterpart on the right side. If $x[n]$ is second-order stationary, then the expected value of each term in each sum is the same; the act of averaging multiple samples that have the same expected value increases the signal to noise ratio. Stated a little differently, you can think the terms in each sum as samples from a new random process:

$$
y[n] = x[n]x[n+kD]
$$

Since the noise present in $x[n]$ is white, the expected value of $y[n]$ is the true autocorrelation of the signal of interest at lag $kD$. Therefore, we would like to accurately estimate the expected value of $y[n]$. Our method for doing so is by calculating a sample average; it can easily be shown that the variance in the sample average estimator decreases given a larger sample size, converging to the actual expected value as the number of samples tends to $\infty$.

So, if there is white noise present in the signal (which is often the case), you're going to get a better estimate of the underlying signal's second-order statistics by using a larger sample size in the calculation (this might sound intuitively obvious). In the context of your two approaches, this is accomplished by using the full, non-decimated signal in the autocorrelation calculation and decimating afterward (i.e. only calculating the result for certain lag values).


Binary classification of grayscale image with little texture I want to carry out a binary classification/segmentation on a grayscale image with very little texture. The only prior knowledge that's supposed to be available is that the object of interest is rather dark (i.e. low gray values) and comprises a significant portion of the image (maybe around 30%-70%). The object of interest should basically be a connected component (although this needn't be the case every single time), but can contain "blobs" of brighter clutter objects. Also, there might be patches of dark pixels in the clutter areas as well, but they should not be as homogeneous and as large as the object of interest.

Of course, I have already checked on simple thresholding (e.g. Otsu's method) or 2-means classification followed by some kind of in-painting, which both lead to similar results. However, these results are only satisfactory if the object of interest is very homogeneous and clearly separable from the clutter object's brighter gray value range. For rather heterogeneous sample images, there is often also clutter, which is put into the object class.

Therefore my question: Is there any interesting literature available about single-channel image classification without exploitation of texture features, but maybe incorporating spatial context? I have already done some research on google, but I have a hard time coming up with proper search terms...`

ADDITION:
Active Contours is very close tho the solution I am looking for. However, they generally work best for one distinct, convex, object, which is enclosed by background. In my situation, however, many convex objects are located in front of the background, so either one active contour for every foreground object is needed or a special case of active contours considering the background as object, being able to deal with clutter objects WITHIN the shape of the object, is needed. Any such method or literature available?

Two large classes of techniques that may be of interest to you:

Mumford Shah functional based techniques.
Active Contour techniques.

Some googling on these terms should open some new doors for you.


Finding Minimum SNR or Minimum Detectable Signal Based on Noise Spectral Density Imagine I have an $n$-bit sensor, whose noise spectral density is $V/\sqrt{Hz}$. $V$ is an exemplary measurement unit. The desired measurement bandwidth is also known. Is there a mathematical way to find a minimum detectable signal or minimum SNR or clean (unaffected by noise) sensor bits or ... based on the available noise density?
As far as I understand, minimum detectable signal concept is normally used for sender/receiver systems. However, for my sensor, I would like to know how powerful the measurand quantity has to be, such that it can be distinguished from the noise.
I will very much appreciate it, if you provide references for your answers.

Let's not get too in the weeds of using bits or noise densities to determine your minimum detectable signal (MDS) just yet. What you're asking is a more fundamental question about determining what value (in terms of SNR) you need to declare a detection. The answer to the question "What SNR do I need in order to detect a signal in noise?" is actually up to you. I'll use a simple radar example here since detection schemes tend to be simpler when compared to more complicated communication systems like PSK, FSK, etc. where you're detecting a sequence of bits or whatever.
In detection theory there are the concepts of probabilty of detection and probability of false alarm. In a simple detector you are trying to choose some absolute threshold value such that

You want to detect the signal(s) of interest (detection)
You want to avoid false detections from noise (false alarm)

Immediately you can see that there are two competing requirements when setting the threshold:

Setting a high threshold will help in avoiding detections from noise, but will also make it harder to detect the signal of interest.
A low threshold will make it easier to detect the signal of interest, but now noise can also be detected more easily causing a false alarm.

An example here is a radar return from a target at 700 m. It has been processed using a matched filter.

The dotted lines are two thresholds we are considering. If we use the blue one, we'll easily avoid false alarms but we're also increasing the chance that the target peak will dip below the threshold, in which we'll miss a detection. The red one will make it very likely that our target will peak above it, but also will the noise which leads to more false alarms.
So, decreasing the threshold increases the probability of detection and the probability of false alarms. Increasing it does the opposite. Having said that, it is up to you to decide which probability to prioritize, since some systems care more about false alarms (or vice versa) and then determine the other. Decide what probabilities you're willing to live with and then set your threshold accordingly. There is no magic number.
There is more to this as this was a simple single pulse example. You can start adding more pulses and play other games, but the basic considerations are the same. You can extend these concepts to communication systems where you can start talking about bit and symbol error rates, etc.


How many taps does an FIR filter need? I am looking to design a set of FIR filters to implement a low pass filter. I am also trying to reduce the latency of the signal through the filter so I am wondering what the minimum number of taps I can use might be.

I know that more taps can lead to a sharper cutoff of the frequency and better stop band rejection etc. However what I'm interested in is more fundamental - if I want to implement a low pass filter with cutoff at $\frac{f_s}{100}$ say does that mean that I need at least 100 taps in order to attenuate the lower frequency signals? Or can I get away with less taps and if so is there some theoretical lower limit?

Citing Bellanger's classic Digital Processing of Signals – Theory and Practice, the point is not where your cut-off frequency is, but how much attenuation you need, how much ripple in the signal you want to preserve you can tolerate and, most importantly, how narrow your transition from pass- to stopband (transition width) needs to be.

I assume you want a linear phase filter (though you specify minimum latency, I don't think a minimum phase filter is a good idea, in general, unless you know damn well what you're going to be doing with your signal afterwards). In that case, the filter order (which is the number of taps) is

$$N\approx \frac 23 \log_{10} \left[\frac1{10 \delta_1\delta_2}\right]\,\frac{f_s}{\Delta f}$$

with

$$\begin{align}
f_s &\text{ the sampling rate}\\
\Delta f& \text{ the transition width,}\\
        & \text{ ie. the difference between end of pass band and start of stop band}\\
\delta_1 &\text{ the ripple in passband,}\\
         &\text{ ie. "how much of the original amplitude can you afford to vary"}\\
\delta_2 &\text{ the suppresion in the stop band}.
\end{align}$$

Let's plug in some numbers! You specified a cut-off frequency of $\frac{f_s}{100}$, so I'll just go ahead and claim your transition width will not be more than half of that, so $\Delta f=\frac{f_s}{200}$.

Coming from SDR / RF technology, 60 dB of suppression is typically fully sufficient – hardware, without crazy costs, won't be better at keeping unwanted signals out of your input, so meh, let's not waste CPU on having a fantastic filter that's better than what your hardware can do. Hence, $\delta_2 = -60\text{ dB} = 10^{-3}$.

Let's say you can live with a amplitude variation of 0.1% in the passband (if you can live with more, also consider making the suppression requirement less strict). That's $\delta_1 = 10^{-4}$.

So, plugging this in:

$$\begin{align}
N_\text{Tommy's filter} &\approx \frac 23 \log_{10} \left[\frac1{10 \delta_1\delta_2}\right]\,\frac{f_s}{\Delta f}\\
&= \frac 23 \log_{10} \left[\frac1{10 \cdot 10^{-4}\cdot10^{-3}}\right]\,\frac{f_s}{\frac{f_s}{200}}\\
&= \frac 23 \log_{10} \left[\frac1{10 \cdot 10^{-7}}\right]\,200\\
&= \frac 23 \log_{10} \left[\frac1{10^{-6}}\right]\,200\\
&= \frac 23 \left(\log_{10} 10^6\right) \,200\\
&= \frac 23 \cdot 6 \cdot 200\\
&= 800\text{ .}
\end{align}$$

So with your 200 taps, you're far off, iff you use an extremely narrow pass band in your filter like I assumed you would.

Note that this doesn't have to be a problem – first of all,  a 800-taps filter is scary, but frankly, only at first sight:

As I tested in this answer over at StackOverflow: CPU's nowadays are fast, if you use someone's CPU-optimized FIR implementation. For example, I used GNU Radio's FFT-FIR implementation with exactly the filter specification outline above. I got a performance of 141 million samples per second – that might or might not be enough for you. So here's our question-specific test case (which took me seconds to produce):
Decimation: If you are only going to keep a fraction of the input bandwidth, the output of your filter will be drastically oversampled. Introducing a decimation of $M$ means that your filter doesn't give you every output sample, but every $M$th one only – which normally would lead to lots and lots of aliasing, but since you're eradicating all signal that could alias, you can savely do so. Clever filter implementations (polyphase decimators) can reduce the computational effort by M, this way. In your case, you could easily decimate by $M=50$, and then, your computer would only have to calculate $\frac{1200}{50}= 24$ multiplications/accumulations per input sample – much much easier. The filters in GNU Radio generally do have that capability. And this way, even out of the FFT FIR (which doesn't lend itself very well to a polyphasing decimator implementation), I can squeeze another factor of 2 in performance. Can't do much more. That's pretty close to RAM bandwidth, in my experience, on my system. For
Latency: Don't care about it. Really, don't, unless you need to. You're doing this with typical audio sampling rates? Remember, $96\,\frac{\text{kS}}{\text{s}}\overset{\text{ridiculously}}{\ll}141\,\frac{\text{MS}}{\text{s}}$ mentioned above. So the time spent computing the filter output will only be relevant for MS/s live signal streaming. For DSP with offline data: well, add a delay to whatever signal you have in parallel to your filter to compensate. (If your filter is linear phase, it's delay will be half the filter length.) This might be relevant in a hardware implementation of the FIR filter.
Hardware implementation: So maybe your PC's or embedded device's CPU and OS really don't allow you to fulfill your latency constraints, and so you're looking into FPGA-implemented FIRs. The first thing you'll notice is that for hardware, there's different design paradigma – a "I suppress everything but $\frac1{100}$ of my input rate" filter needs a large bit width for the fixed point numbers you'd handle in Hardware (as oppposed to the floating point numbers on a CPU). So that's the first reason why you'd typically split that filter into multiple, cascaded, smaller, decimating FIR filters. Another reason is that you can, with every cascade "step", let your multipliers (typically, "DSP slices") run at a lower rate, and hence, multiplex them (number of DSP slices is usually very limited), using one multiplier for multiple taps. Yet another reason is that especially half-band filters, i.e. lowpasses that suppress half the input band and deliver half the input rate, are very efficiently implementable in hardware (as they have half the taps being zero, something that is hard to exploit in a CPU/SIMD implementation).


Recover signal from its power spectral density I was wondering whether we can recover a signal $x(t)$ from its power spectral density using signal processing techniques.

Generally (in the absence of prior knowledge) , I don't think this could be done.

The PSD is the absolute value  squared of the DFT of the signal, and the abs function is not 1-to-1. For example abs(x) could be produced by f=x and f=-x.

However, if you have some prior knowledge about the signal then you could reconstruct it from the PSD. Without more knowledge of what you're trying to do I don't think I could help you.


Why does applying RRC to the signal increase bandwidth? Please see this question: Why samples per symbols $\geq 2$ in the GNURadio Constellation Modulator block?.  The first answer here says "filtering the signal, you have bandwidth expansion".  But look at the figure in this question: Relation between Bandwidth and Baud-Rate for 8-PSK - it looks like there is much less bandwidth required for the RRC filtered signal than the square wave / sinc.  So why does the first question have an answer that says that increased bandwidth is necessary and, taking it a step further, why does sps increase, as was the actual question there?

I wouldn't really call it "bandwidth expansion". If you do 2 samples per symbol, from a purely theoretical point of view (i.e. using perfect infinitely long sinc pulses), you could do with exactly half of the Nyquist bandwidth of the output signal.
However, RRC is not infinitely long sincs, so you get excess bandwidth.


FFT error when number of samples growth - Matlab/octave I see when I increase number of samples, the imaginary error part of the FFT growth due to the error of sin(n*pi)!=0 where n is the index of the sample in Matlab or octave.
Is there anyway to force Matlab/Octave to perform FFT without these type of errors?
for example I know if you have sin (sym(pi)*n) =0 in Matlab or octave. Why default FFT in matlab/octave does not use sym(pi) instead of pi?

How could I perform FFT that take care of such PI approximation error?
Thanks

An FFT using symbolic math might be possible, but would be many orders of magnitude slower.  (I'm guessing at least 10,000 times slower, except for a set of exact equation input signals).  You would have to use a symbolic math package instead of Matlab.  (Perhaps Mathematica, Maxima, or Maple?)

Instead you might be able pre- and post-process certain inputs to produce known outputs.  For instance, you can check if the input is strictly symmetrical around input sample 0.  If so, set all the imaginary components of the FFT output to zero.  If strictly odd (anti-symmetric), set all the real result components to zero.  If strictly real, set all the negative frequency outputs to be exact complex conjugates of the positive frequency ones. etc.