Is it possible to make 3D point cloud with RADAR? Is it possible to make 3D point cloud with RADAR? Both LIDAR and RADAR are EM waves, but with different wavelength. So, if it is technically possible to make point clouds with LIDAR, why is it not possible to do the same with a RADAR?

Is it because, the RADAR cannot have narrow beams like LIDAR? Or the problem is more related to the antenna? Also, why isnt doppler effect experienced on LIDARs?

In short, the question is how are RADAR and LIDAR technically different with respect to sensors?

A point cloud is related to how data is acquired and stored.  One can usually transform between point clouds and rasters and say, since they can be be used in similar ways, the answer is essentially yes.  The most similar type of RADAR to a LIDAR point cloud is a Synthetic Aperture RADAR which is just called SAR.

and contrast it with a LIDAR Image , centered on the same place.  This might not be the best comparison but adequate for a 3 minute Google search.

clearly LIDAR has higher resolution.  The SAR image is formed from a coherent process of data acquired over a scan. LIDAR points are one time, travel time shots.  Each measurement adds a point.  I'm not aware of coherent LIDAR processing but I'm not State of the Art LIDAR researcher.  Doppler has an effect on the phase of the signal, and an incoherent energy detector is insensitive to phase.  You could take a SAR image and turn it into something like a point cloud but the way the data is collected and processed is different, so maybe your answer is no.

LIDAR, used like this are usually looking down over a narrow range of angles under the aircraft. SAR can more offset over a wider range of angles.  The SAR image is probably from a single flight path.  A LIDAR Image of the same extent would take a number flight paths.  One can use SAR in orbit.  I'm not sure about LIDAR through a lot of atmosphere.

SAR will probably have the advantage with respect to transmission power for some time.

The highest resolution global digital elevation maps that are publicly available are derived from SAR data.

https://en.m.wikipedia.org/wiki/Shuttle_Radar_Topography_Mission

LIDAR derived elevation data is growing.

requested edit:

A raster is a regular sampling, like uniform sampling.  Each point in a SAR image is produced by a calculation involving data that was collected over a time window, that for a pulsed waveform, would span multiple pulses. One needs to select those solution points.  Adjacent solution points overlap in the raw collected data input to the calculation.  SAR raw data can include polarization.

A point cloud is a set of individual measurements and the points aren’t regularly spaced.

SAR is a interferometric or coherent processor. SAR has some history to it and was first processed using analog processing.  Digital processing is the norm today.

Im not familiar with coherent Lidar but the term heterodyning, in my experience is analog processing. Cloud sampling isn’t a necessary theoretical requirement for Lidar.  The article that is cited in the comment speaks of continuous wave coherent Lidar which would seem to me a raster sort of sampling.

Addressing what one can do that the other can’t,  has categories like can’t now, and can’t never.  The can’t now category shifts with technology. The can’t never ignores the possibility of an ultra Wideband, ultra agile, hybrid system.

The missing comparison to Lidar is stereo vision and even in that case, they aren’t mutually exclusive.


Calculation of sampling rate and its effect on an Image? I know that for a given signal, the sampling frequency $F_s$ must be twice or more than maximum frequency of the signal $F_m$. It is easy to understand the concept for a 1D signal. But I don't know how to calculate sampling frequency or Nyquist rate for a 2D image.

Also what is effect of over-sampling and under-sampling in case of image shown below (dimensions $\text{width} \cdot \text{height}=204\cdot226$)?

Depends on how the image is created.  If by photography, the focus cone and Airy disk for the given lens have a somewhat low-pass filter effect.  Figure out the effective transition frequency at the imaging plane, and have a sampler with a bit over twice that resolution to diminish aliasing effects during image capture.

If you undersample (either during capture or with later processing) without low pass filtering, then you may get alias artifacts such as Moire patterns on high contrast fine detail.  Classic example is a pin striped shirt at SD TV resolution.


Does speaker recognition system need the same sample rate audio? I'm a new man in speaker recognition project. I've train a model with some .wav files, its sample rate is 16000.

I want to test this model with my record, but my record's sample rate is 8000

Doed it need the same sample rate for training and testing speaker recognition system?

Hope to have your answers,

Thanks

Does it need the same sample rate for training and testing speaker recognition system?

I suppose that you are asking whether or not the difference in sampling frequency ($Fs$) will affect your classification result.

The short answer is "It depends on the data". The slightly longer answer is:

The $Fs$ has a direct impact on the captured bandwidth. At 8kHz, you are capturing (nominally) 0-4kHz. At 16kHz, you are capturing (nominally) 0-8kHz.

Speaker recognition is essentially trying to discriminate between the spectral differences in the sound of different speakers.

Therefore, a trivial example is one where if your $Fs$ is 8kHz and the key variation of spectral differences occurs between 6-7kHz, then your classifier will completely fail to discriminate between its classes.

Another detail is being sure that the $Fs$ is suitable, i.e. it includes the frequency range of interest, but, the recordings are awful and some other sound source is masking the range of frequencies the classifier depends on.

Therefore, "It depends on the data".

But, if we assume that the recordings are not problematic and a lower $Fs$ does indeed contain the frequency band of interest, then what might be worth considering more is the size of the transforms that are involved in the derivation of the MFCC that your classifier depends on. In other words, lowering the $Fs$ might require that you increase the length of the transforms to maintain the same (or similar) frequency resolution. Otherwise, adjacent frequency bins might get lumped together and reduce the clarity of the observed spectrum.

Finally, a more practical note would be to make sure that any existing code you might be using does not contain hard coded $Fs$. If the original author has assumed a fixed $Fs$ throughout the implementation of the model, it is likely that the classifier will have unpredictable behaviour if you feed it with sound files of different characteristics.

Hope this helps.


Learning about inverse-z-transform and how to apply it to a rational transfer function I have been studying IIR filters and know that a rational transfer function:

$$
H(z) = \frac {b_0 + b_1 z^{-1} + ... + b_N z^{-N}}{1 + a_1 z^{-1} + ... + a_N z^{-N}}
$$

has a finite difference equation:

$$
y[n] = b_0 x[n] + b_1 x[n - 1] + ... + b_n x[n - N] - a_1 y[n - 1] - ... - a_N y[n - N]
$$

What I am confused about is how the transformation from transfer function to difference equation is performed.

I have read about the z-transform and the inverse z-transform and understand that the inverse z-transform is the key, but none of the internet sources I have read provide a joined up walk though of how to take a transfer function and derive a finite difference equation from first principles.

Could anyone recommend a good source, or explain the end-to-end approach?

You do not really need the inverse $\mathcal{Z}$-transform to derive the difference equation corresponding to a given rational transfer function. You just need to know that $z^{-k}$ corresponds to a delay of $k$ samples. The inverse $\mathcal{Z}$-transform of $H(z)$ would give you the corresponding impulse response, but you don't need the impulse response to write down the difference equation.

Knowing that

$$H(z)=\frac{Y(z)}{X(z)}$$

you can rewrite the equation given in your question as

$$Y(z)\left(1+a_1z^{-1}+\ldots +a_Nz^{-N}\right)=X(z)\left(b_0+b_1z^{-1}+\ldots +b_Nz^{-N}\right)\tag{1}$$

With the correspondences $Y(z)z^{-k}\Longleftrightarrow y[n-k]$ and $X(z)z^{-k}\Longleftrightarrow x[n-k]$ you can directly write down the difference equation from $(1)$:

$$y[n]+a_1y[n-1]+\ldots +a_Ny[n-N]=b_0x[n]+b_1x[n-1]+\ldots +b_Nx[n-N]$$

Here's a way to show that $\mathcal{Z}\{x[n-k]\}=X(z)z^{-k}$:

$$X(z)z^{-k}=\sum_{n=-\infty}^{\infty}x[n]z^{-n}z^{-k}=\sum_{n=-\infty}^{\infty}x[n]z^{-(n+k)}=\sum_{n=-\infty}^{\infty}x[n-k]z^{-n}$$


Clear examples of Fixed Point Algorithm Specification Block Diagram for FPGA Implementer I have been working on converting a simple, floating point algorithm to fixed point and need to create a block diagram to help the guy who will implement this on an FPGA.  As often is the case, the output can be represented as a linear combination of (the real and imaginary parts) of some input signals.
I would like to create a clear block diagram that shows bit-widths, rounding vs truncation, saturation at each point, bit shifts (due to multliplying integers with fractions) etc. (Do you use Qmn notation in the diagram?). Most of what I've seen online is either so detailed with HW registers, etc that I can't actually tell what's going on or so high level, that it's essentially a conceptual representation of the algorithm with some annotations containing bit widths.
While I recognize that there is no one right answer, I'm curious to know how others create good diagrams to hand the algorithms off to implementers.
Thanks.

I designed a lot of DSP algorithms and I would usually be the one implementing them in an FPGA. I rarely had to explain the internals of an algorithm to other people. That beind said, it is a good idea to represent the number of bits (I usually used the s:m:f notation from Sony Playstation 2) in a block diagram/design document for every components. You should also specify whether or not the FPGA designer must truncate or round.
Otherwise if you leave the number of bits used for the internal calculations up to the FPGA designer, you're looking for trouble. Using an insufficient number of bits can render your algorithm useless or unstable.
It is also a good idea to create a fixed-point model of your algorithm in Matlab/Python/Octave/C++ for 2 reasons. First, you will be able to validate the number of bits needed so that your fixed-point algorithm behaves similarly to the floating-point arithmetic algorithm. Secondly, you will be able to compare the FPGA implementation to your fixed-point model. The FPGA implementation should behave exactly the same as your fixed-point model.


Using phase information to extract sample values I have one cycle of the following signal stored in a memory:

 x(n)=sin(2 * pi * n/N + theta)

 where  theta= 2*pi*q/N and q, N are integers.

How do I use these values to obtain sinusoides (samples) with same frequency but different phase?

If the values are stored in a linear array indexed by $n$, read the values out
sequentially from a different starting point in the array, that is beginning with a nonzero value of $n$, and have the read out "wrap" around the end of the array, that
is, the incrementing of the index is done modulo $N$. This gives sinusoids
with possible initial
phases $0, 2\pi/N, (2)2\pi/N, (3)2\pi/N, \ldots, (N-1)2\pi/N$ depending on the
choice of starting point of the read out.  If you want other values
for the initial phase, you will need to do some computations, and not just
a plain reading out from memory.


Combining HRTF with echo I am experimenting with HRTFs.

Specifically I'm trying to create 12 impulse responses that the brain can easily differentiate.

The more distinct they are from one another the better.

So I'm trying to figure out ways of making them as distinct as possible.

My best attempt so far is to use HRIRs on the horizontal plane with varying azimuths.

I could use HRIRs for 12 equally-spaced points (i.e. azimuth 0 30 60 90 ... degrees),  but human beings can perceive tiny differences in angle directly in front (as opposed to the sides), so a more sensible spacing would be something like 0 10 20 45 90 135 160 170 180.

That's the sort of thing I'm looking into...  I give it as an example as I would be grateful for guidance in this more general enquiry, and also it shows my motivation for this question: maybe I should drop this line of enquiry as it won't help?

Anyway, I'm wondering whether HRIRs from an echoic environment may present more cues, and hence provide better localisation.

And if this is the case, is there any shortcut to actually recording the HRIRs in a cathedral (which has actually been done: http://www.ind.rwth-aachen.de/en/research/tools-downloads/aachen-impulse-response-database/)?

i.e. is it possible to take an impulse response for a cathedral,  and convolve it with each HRTF?

π

Choosing points in the horizontal plane with proper spacing seems like the right approach. Aneochic HRTFs are probably the easiest but typically lead to "inside the head" localization. Depending on your application this may or may not be a good thing. If you want externalization it's best to use binaural room impulse responses measured in a space that's similar to where the actual playback happens (there are psychological reasons for that). However these binuaral room impulse responses are a lot more complicated and require significantly more computational effort.

You cannot binauralize a measured monaural impulse response after the fact. A real binaural room impulse response consists of the direct sound plus many reflections that all come from different directions. So each part of the room impulse response contains many different anechoic HRTFs. There are however computer programs that allow you to put in a geometry, acoustic properties and then they will calculate a binaural room impulse response.


What is the need for prediction filter in PCM and DPCM? The DPCM works by the difference between actual samples and predicted samples.

If we already have the actual samples, we can quantize it and encode it later.

But why we use prediction filter?
And by the definition "predicts the present samples with previous samples" how does the filter do that?

The DPCM works by the difference between actual samples and predicted samples.

  If we already have the actual samples, we can quantize it and encode it later.

  But why we use prediction filter?

the idea is to transmit the D (or $e[n]$) and let the receiver add that to the previous PCM sample it has (or can predict from the previous PCM samples using the same predictor).  and the reason was that, if the prediction was good, then $e[n]$ would be small and a smaller word width was needed.

  And by the definition "predicts the present samples with previous samples" how does the filter do that?

the idea is that a good prediction of $x[n]$ is a linear combination of $x[n-1]$, $x[n-2]$, $x[n-3]$, $x[n-4]$... $x[n-N]$.  so some guys in the sixties just suggested that the prediction is:

$$ \hat{x}[n] \triangleq a_1\, x[n-1] \ + \ a_2\, x[n-2] \ + \ a_3\, x[n-3] \ + ... + \ a_N\, x[n-N]   $$

the prediction error is

$$ e[n] \triangleq \hat{x}[n] - x[n] $$

and this $e[n]$ is what is transmitted and the next PCM sample can be constructed from the previous PCM samples and the error as

$$ \begin{align}
 x[n] &= \hat{x}[n] \ - \ e[n] \\
      &=  a_1\, x[n-1] \ + \ a_2\, x[n-2] \ + \ a_3\, x[n-3] \ + ... + \ a_N\, x[n-N] \ - \ e[n]
\end{align} $$

you can see that the prediction is a linear FIR filter and that the reconstruction has feedback and is a linear IIR filter (that happens to have all of the zeros at the origin, sometimes called an "all-pole filter" which is slightly imprecise).

then the idea is to, with knowledge of the statistical spectral characteristics of $x[n]$, namely its autocorrelation, find the values of $a_1$, $a_2$, $a_3$, ... such that the prediction is the best, that is such that $|e[n]|^2$ is minimized in a statistical manner.  this is what linear predictive coding (LPC) is at the most fundamental.  smaller $e[n]$ means that the word width needed for that signal is smaller and that the number of bits you need to send $e[n]$ over a transmission channel is less.  and, it means you can add Shannon and Huffman to the mix and send even fewer bits over the channel.  that is one way to do lossless coding.


Perform non-power-of-two FFT using ARM CMSIS library I have 1280 points which I would like to perform an FFT using a M7 cortex with hard floating point engine.

As the title suggest I would like to use the ARM CMSIS library to perform a non-power-of-two FFT.

As I am not familiar with the library functions or FFTs, I would like to know if this is possible based on the available CMSIS FFT functions.

Can i combine or use the available radix-2, 4 and 8 functions to 'build' what I want? What particular FFT algorithms should i be looking at?

Other options I gather would be to zero pad the signal from 1280 to 2048.

The application is for 3-phase Power System Voltage and Current monitoring (50/60 Hz) at 6.4 kHz sampling rate.

I would like to average the Harmonics over time (only the integer multiple of fundamental comp) and compare their magnitudes against industry regulated threshold.
Perform a THD on the harmonics.
I would like to extract the fundamental component for further processing (symmetrical components
calc and unbalance ratios).

Here is some code that calculates a 1280 tap FFT based on 5 256 tap FFTs

%% 1280 FFT based on five FFTs of 256 each
n = 1280;
% Create a piece of noise
x = randn(n,1);
% calculate FFT using MATLAB native fft() function.
% We'll use this as a reference to prove it works
fx = fft(x);

% Break down into five signals of 256 points each, interleaved
p = x(1:5:end);
q = x(2:5:end);
r = x(3:5:end);
s = x(4:5:end);
t = x(5:5:end);
% FFT each of those. This is a 256 power-of-two standard FFT
fp = fft(p);
fq = fft(q);
fr = fft(r);
fs = fft(s);
ft = fft(t);
% Do five times periodic extension (just repeat it 5 times)
fp5 = repmat(fp,5,1);
fq5 = repmat(fq,5,1);
fr5 = repmat(fr,5,1);
fs5 = repmat(fs,5,1);
ft5 = repmat(ft,5,1);
% calculate the 1280 twiddle factors
k5 = (0:n-1)';
W5 = exp(-1i*2*pi*k5/n);
% assemble the result
fy5 = fp5 + W5.*fq5 + W5.^2.*fr5 + W5.^3.*fs5 + W5.^4.*ft5;
% calculate the error
ferror = fy5-fx;
fprintf('Error = %6.2f dB\n',10*log10(sum(ferror.*conj(ferror))./sum(fx.*conj(fx))));


Graphical fourier series of a square wave This is probably off-topic since it isn't really a question, but I thought that this GIF of the fourier series of a square wave was too cool not to share.

Here is an interactive version of the GIF, where you can test how partials add up to shape a certain wave: http://toxicdump.org/stuff/FourierToy.swf


Downsampling and Gaussian Filtering in the Context of Scale Space Pyramids In the context of scale space image pyramids using Gaussian filters I noticed that it's common to downsample the image after blurring with $\sigma = 2*\sigma_{init}$ . My question is: What is the relation between variance and downsampling? What is the proof that I'm not losing information when the image is downsampled after blurring with this factor?

I think this can easily be shown using Nyquist and the cut-off frequency of the Gaussian filter, but I'm not sure how.

You're correct, it has to do with the Cut Off frequency of the Gaussian Blur Filter in its Frequency Domain.

In order to see it, just apply a DFT (Using MATLAB it can be achieved by fft / fft2) and look on the absolute value.
Look for the -3dB point and you'll see.

There is also an intuitive explanation on the original article which say that blurring using Gaussian Kernel with $ \sigma $ means removing any detail which its size is smaller than $ \sqrt{ \sigma } $.


denoising floor function Suppose that

$y(t) = \lfloor x(t) + \epsilon(t) \rfloor$,

where $\epsilon(t)$ is zero mean, independent noise.

Is there any techniques on recovering $x(t)$ from $y(t)$?

The problem is a bit more complicated that mere denoising, since a highly non-linear operator affects the original signal. Restoration or recovery would be more appropriate.

The floor function acts as an instance of data quantization (and possibly saturation or clipping). And adding noise is also known as dithering, a method to  reduce sensitive quantization effects (dequantization).
Those are keywords you can use along with: restoration, recovery, inference... Without more information on the signal's dynamics, models, and the noise, a couple of pointers that I am aware of (to be updated):

Wavelet denoising of coarsely quantized signals
quantized compressing sensing
piecewise regular approximations


Apply "non-causal" filter buffer-wise, a.k.a "soft real-time filtering" I am dealing with digital filtering of signals, both offline and in real-time. Typical filtering purposes are highpass filter or bandpass filter.
So far I worked on prerecorded signals (e.g. wav files) so I could use a non-causal filter, such as Python's SciPy filtfilt or MATLAB's filtfilt, which produces zero-phase filtering (see What is zero-phase filtering and forward-backward filtering?). This has the desired effect of canceling any phase distortion. So far it worked very well off-line.
However, since according to theory this kind of filter is non-causal (e.g. the implementation of filtfilt has a stage in which the signal is filtered backwards from its end to its start), it can't be used for "hard" real-time filtering (i.e. reading single sample at a time and immediately filtering it).
I ask if this non-causal filter can be used for "soft" real-time filtering, that is: I read a buffer of n samples, then filter the buffer using filtfilt in the time-domain or by using its frequency response (which is purely real, obtained from the filter by the SciPy function freqz) in the frequency domain (via FFT), and then I read the next buffer and so on.
Another question is what is the recommended way to apply this filter: in the time-domain via filtfilt or in the frequency domain (multiplying the FFT coefficient $F(k)$ with the real positive frequency  response $|H(k)|$)? The purpose is to have a minimally-distorted filtered signal and fast and computationally efficient filtering (to run in "soft" real-time).
Added in edit:
In my application a filtering in the frequency space is desired, and multiplying $F(k)$ by $|H(k)|$ seems to work, even in real-time. But multiplying by the absolute value is like performing zero-phase filtering, which by theory is non-causal. However, since I apply it buffer-wise, this seems to be the loophole. Am I correct?
Thanks.

Overlap-add or overlap-save/scrap with zero-padded data are the common methods of using block based convolution on streaming data.  Pad the convolution (FFT/IFFT fast, or linear) by at least the length of the impulse response above your desired noise floor, minus 1.  The basic idea is that these methods save the remainder of the impulse response that doesn't fit in the current block's filtered results, and applies that remainder to blocks that are processed later.  For linear phase FIR filters, this still incurs as least half the length of the impulse response in latency.


Digital Signal Quantized Processor Complexity In “traditional” digital signal processing, the complexity is computed as a number of multiplication the operation requires, e.g. the computation of the N-point DFT via the decimation-in-frequency FFT requires (N/2)log_2(N) complex multiplications (and the same number of additions, but in the DSP the additions are assumed to require less time and their contribution is traditionally omitted, unless specifically required). However, we typically do not define which numbers we multiply. My question is how we define the complexity in the DSP sense, if we apply some quantization rule for the numbers to be multiplied: e.g. we proceed from multiplying two real numbers in all operations of the processing devise to multiplying the pairs of int8, float 16, etc. The latter operations are evidently “simpler” that the multiplication of two reals, but how can we quantify and express this complexity reduction for the DSP element?

First, if there's a more formal way to do this than what I'm going to show, I don't know it.

Call out the MAC operations along with data path widths.  I.e., "the FFT algorithm used takes $4 n \log n$ operations with 24-bit wide data paths"

A corollary to this would be that sometimes you may have different precision in your coefficients and your actual data.  So you may call out "24 bit signal path, 16 bit coefficients".

Call out the LUTs used (in an FPGA).  This isn't direct, and has a lot to do with the cleverness of the implementation, but within a specific FPGA technology it has meaning.

Note that these days, unless you're doing full-custom ASIC design, data paths come in chunks - processors give you multiples of their native word width, and FPGAs come with "DSP blocks" that come in fixed widths.  I'm only guessing, but I'd bet money that if you go to design an ASIC then at the first level of design complexity there will be signal processing libraries with their own fixed-width data paths, and that libraries with variable-width data paths would take more engineering resources to use.


Why is signum function used to calculate Fourier transform of unit step function I read in a standard textbook that the Fourier transform of unit impulse function is calculated with the help of approximations and signum function as the integration of unit impulse does not converge. What's so special about signum function that it is used to calculate Fourier transform? I tried to find out an approximation as:

$$
\lim_{ a \rightarrow 0 }
\int_{-\infty}^{+\infty} e^{-at} u(t) e^{-j\omega t} dt
$$

But I am getting wrong result. Why is this so?

If somebody you trust told you that the Fourier transform of the sign function is given by

$$\mathcal{F}\{\text{sgn}(t)\}=\frac{2}{j\omega}\tag{1}$$

you could of course use this information to compute the Fourier transform of the unit step $u(t)$. Using

$$u(t)=\frac12(1+\text{sgn}(t))\tag{2}$$

(as pointed out by Peter K. in a comment), you get

$$\mathcal{F}\{u(t)\}=\frac12\left(\mathcal{F}\{1\}+\mathcal{F}\{\text{sgn}(t)\}\right)=\pi\delta(\omega)+\frac{1}{j\omega}\tag{3}$$

However, you don't need the sign function to compute the Fourier transform of the step function. As suggested in your question, using the function $e^{-at}u(t)$ and taking the limit $a\rightarrow 0^+$ will also result in the expression given in $(3)$.

You can see this as follows. The Fourier transform of $e^{-at}u(t)$, $a>0$, is given by

$$\int_0^{\infty}e^{-at}e^{-j\omega t}dt=\frac{1}{a+j\omega}\tag{4}$$

Taking the limit $a\rightarrow 0^+$ appears to give $1/j\omega$, but this is only valid for $\omega\neq 0$. Splitting the result $(4)$ in its real and imaginary part gives

$$\frac{1}{a+j\omega}=\frac{a}{a^2+\omega^2}+\frac{\omega}{j(a^2+\omega^2)}\tag{5}$$

The real part of $(5)$ is known as ($\pi$ times) a nascent delta function. It has the same form as the Poisson kernel, which in the limit becomes a Dirac delta impulse. So for $a\rightarrow 0^+$ the limit of $(5)$, and hence of $(4)$, is actually given by

$$\lim_{a\rightarrow 0^+}\frac{1}{a+j\omega}=\pi\delta(\omega)+\frac{1}{j\omega}\tag{6}$$

which equals the expression in $(3)$.


Differential BPSK vs BPSK what are the advantages and disadvantages? Recently I have come across different variation of the sane modulation techniques. I am quite confused about the advantages of different constellations. In particular, what are the advantages of using

Differential BPSK over BPSK?

Any thoughts?

Differential BPSK

As the information is encoded in the phase difference of two consecutive symbols, phase estimation can be omitted in systems where the carrier phase can be assumed to be constant during two symbols. Specifically, it can be used in a self-coherent receiver which mixes one symbol with its predecessor in order to reconstruct the phase difference between the two symbols. In such a scheme, no local oscillator is required which simplifies the system (e.g. frequency synchronization is not needed).

Non-differential (coherent) BPSK

In differential BSPK, two symbols are "added" and as a consequence the noise of the two symbols adds up. This results in a 3 dB sensitivity gain of coherent BPSK. In other words, coherent BPSK is more robust against noise.


Implications of $X( -j\omega ) = X^*(j\omega)$ What are the implications of:

  If $x(t)$ is real and $x(-t) = x^*(t)$, then $X(-j\omega) = X^*(j\omega)$ and $X(j\omega)$ is real.

I am trying to understand it and I would like to research it further (what's the theorem that guarantees this?).

I'd like to understand why it's important aside from the fact that it guarantees a symmetric Fourier transform or frequency response if we're talking about that of the impulse response.

Thanks,

The fact is that if $x(t)$ is real and then $X(-j\omega) = X^*(j\omega)$ (it's easy to prove) and, because of duality, if $X(j\omega)$ is real, so also then $x(-t) = x^*(t)$.  in your case, you have both $x(t)$ and $X(j\omega)$ being real.  keep in mind that anything that is purely real is also equal to its complex conjugate.

so you can say:

"if $x(t)$ is real and even (that is $x(-t) = x(t)$) so also is $X(j\omega)$ real and even ($X(-j\omega) = X(j\omega)$)."


What is the physical meaning of the cross-spectrum of two signals measured at two locations and different moment? I am working with steady state pressure field measurement where sound pressure at multiple locations in the field needs to be measured. In my case, more than one microphone existing in the field at the same time will intrude the field itself. As a result, the pressure of each location needs to be measured at different moments.

My question is is it still physically meaningful to calculate the cross-spectrum of signals of two microphones? If so, what is the physical meaning of the phase of the cross-spectrum at one given frequency?

I understand that for two signals that are measured at the same time, the phase of the cross-spectrum is the phase shift due to the time lag. But when both time and space is involved, I am confused on the meaning of the phase from cross-spectrum calculation.

I figured this out myself. The answer to this is that for the steady state, phase from the cross-spectrum is a contribution of both time and space. And if two signals are measured at the same time or location, the contribution of the other factor to phase will be eliminated. The conclusion is made based on formula deduction, which I will not discuss further here.


Relation between Laplace and Fourier transforms I know that  $$X_L(s=j\omega)=X_F(\omega)$$ if $x(t)$ is one sided and absolutely integrable and hence the imaginary axis of the Laplace transform is the Fourier transform.

But Fourier transform also has imaginary and real parts. So how could this be right?

The Laplace transform evaluated at $s=j\omega$ is equal to the Fourier transform if its region of convergence (ROC) contains the imaginary axis. This is also true for the bilateral (two-sided) Laplace transform, so the function need not be one-sided.

As for real and imaginary parts, since $s$ is a complex variable, both the Laplace and the Fourier transform generally have real and imaginary parts. Take as a simple example the function $x(t)=e^{-at}u(t)$, with $a>0$, where $u(t)$ is the unit step function. The Laplace transform is

$$X_L(s)=\frac{1}{s+a}\tag{1}$$

Since $a>0$, the ROC of $X_L(s)$ contains the imaginary axis, and the Fourier transform of $x(t)$ is simply obtained by evaluating $X_L(s)$ on the imaginary axis $s=j\omega$:

$$X_F(\omega)=X_L(j\omega)=\frac{1}{j\omega+a}\tag{2}$$

Since $s=\sigma+j\omega$ is generally complex, not only the Fourier transform but also the Laplace transform $(1)$ has a real and an imaginary part:

$$X_L(\sigma+j\omega)=\frac{1}{\sigma+j\omega+a}=\frac{\sigma+a}{(\sigma+a)^2+\omega^2}-j\frac{\omega}{(\sigma+a)^2+\omega^2}\tag{3}$$

Only when evaluated on the real axis $s=\sigma$ ($\omega=0$) does the imaginary part vanish.


Doubt in Frequency Modulation code While searching for something FM, I came across this code. I am not able to understand a particular line of code. I have written a comment against that line of code.

fs = 1000;
fc = 200;
t = (0:1/fs:0.2)';
x = sin(2*pi*30*t)+2*sin(2*pi*60*t);
fDev = 50;
int_x = cumsum(x)/fs;  %%% This line of code
xfm = cos(2*pi*fc*t).*cos(2*pi*fDev*int_x)-sin(2*pi*fc*t).*sin(2*pi*fDev*int_x)   ;

I understand the cumsum function but I am not able to understand why is he dividing the integration by Fs i.e. sampling frequency.

Link to the code - https://in.mathworks.com/matlabcentral/answers/318969-fm-modulation-without-using-fmmod

When generating the sine wave, the code uses a sampling frequency to establish discrete points in the time domain. It is not continuous. This means the wave is broken up into bins that are exactly 1/fs in size/width.

So if you want a true integral of the wave, you have to take the average density of all the points in a bin by dividing by the sampling frequency.

Hope that helped.


Creating 2D array to play audio in left and right channel in matlab I'm new to matlab sound processing and I am trying to create a 2D array for the sound code below. I have to make the sound have 400 samples(i.e 10ms) at 20hz, but having trouble getting the code to work. Any help would be greatly appreciated.
The steps are:
(1)    Construct a two-dimensional array (n x 2  array, i.e., 2 columns, n rows, n is the number of samples per sound channel).  The array is the sound signal, the first column is the sound for the left channel (left speaker), the second column is the sound for the right channel (right speaker)
(2)    Then play the sound.

Fs = 44100;
t = [0:4/Fs:4-4/Fs];
freq = 20;
f1 = sin(2*pi*freq*t);
f2 = sin(2*pi*2*freq*t);
%Sound 1
sound(f1,Fs);
%Sound 2
sound(f2,Fs);
%Play Consecutive
f12 = [f1 f2];
sound(f12,Fs);
%Play together
f_12 = [f1+f2];
sound(f_12,Fs);

Since 20 Hz was a little hard to hear with my speakers, I've changed the main frequency to 600 Hz. I guess it should work:

clc;
clear all;

Fs = 44100;
N = 40000; % sample size
t = (1:N)*(1/Fs);
freq = 600;

sound_array = zeros(N, 2); % initializes a N x 2 matrix
sound_array(:,1) = sin(2*pi*freq*t);
sound_array(:,2) = sin(2*pi*2*freq*t);

% Sound 1
sound(sound_array(:,1),Fs);

% Sound 2
sound(sound_array(:,2),Fs);

% Play the first column at left channel and the second column at the right
% channel
sound(sound_array,Fs);

In fact, if you use the entire array as input in sound, the first column will go to the left channel and the second column to the right channel.


DWT on short signals I have signals with a very short rate of ~100e-6 seconds, sampled at 10MHz,so about 1k samples.
I am wondering, is there a practical limit to the usefulness of time-frequency methods concerning a minimal duration of the signal or a certain sampling frequency for that matter?
i am talking about chirp radar signals which fires 1000 linear chirps each 100us long which I sample at 10 MHz. so the whole frame takes only about 100ms

With such short (known) stimuli, does it make sense to use generic analysis methods, or is it «better» (less work, results that are more usable) to tailor your processing to the transmitted processing (matched filter?) or apply a parametric model and some numerical optimization?


Why is there little emphasis on phase distortion? I have a question regarding the frequency response of digital filters. Specifically the phase aspect of the frequency response.

Why do we not seem to care about the nonlinear phase response in certain applications?
Is it possible to correct the phase distortion if we are filtering offline (i.e. not in realtime)?
Given that we have a model of our filter, can we take the outputted phase distorted signal, and then individually shift the frequencies within the signal?
It also seems that for applications where we can't allow phase distortion, we often implement FIR filters. What are some applications where FIR filters would be absolutely necessary? Essentially anything being filtered in realtime?

Why do we not seem to care about the nonlinear phase response in certain applications?

Very likely because these certain applications don't care about phase!

  Is it possible to correct the phase distortion if we are filtering offline (i.e. not in realtime)?

With another filter that, combined with the distorting filter, has the desired phase properties, yes. The question whether such a filter exists is actually a complicated one and depends on the distorting filter!

  Given that we have a model of our filter, can we take the outputted phase distorted signal, and then individually shift the frequencies within the signal?

This is where the complications come from: If you pick "single" frequencies only, you could only represent perfectly periodic, infinite signals (limit of bandwidth-time product in the Fourier transform). So you'll need some kind of system that is "continuous enough".

  It also seems like for application where we can't allow phase distortion, we often implement FIR filters. What are some applications where FIR filters would be absolutely necessary? Essentially anything being filtered in realtime?

No; real time isn't the point here.

So, there's many applications where you want a linear system.

Which practically reduces your choices to IIR and FIR filters.

However, the question whether to choose either of these is often one of design possibilities (I personally find it relatively easy to design a FIR with the window methods, and relatively hard to design good IIRs, but that's me, who's grown up in a FIR-friendly environment); it's often also a question of how stable an IIR design can stay under limited bit width assumptions, or how long a FIR becomes.

Overall, the choice of digital filters is a along topic, that actually fills shelves in libraries. It's impossible to just say "we use that", without knowing a lot of the requirements set for a filter.

Having said that, in many cases you just really need something that is easily adaptable. FIRs lend themselves to that excellently, because you'd usually "just" have to exchange some tap values in software, and depending on the actual implementation of the FIR, you might even have something like a continuous response during that exchange.

Note that all questions regarding which type of filter you use typically involve tradeoffs – your filter response isn't perfect, but you can't implement infinite length FIRs anyway. Your latency is low, but you paid for that with a suboptimal suppression in the stop band. Your decimating filter has adjustable decimation, but in exchange its CIC response is rather round.

A lot of these tradeoffs involve the actual thing you need to implement your filter on. You really shouldn't assume that an experienced DSP hardware designer would just take your Matlab FIR and put it in his FPGA/ASIC design. First of all, there's typically no floating point in those designs, so he'd first need to check dynamic range, convert to fixed point, estimate the error incurred by that, adjust for that, and secondly, things that might be hard in matlab might be easier in hardware; changing the path that samples travel based on the sign bit at Megasamples per second? Will totally thrash CPU performance if you constantly exchange taps in CPU caches, but for FPGA designs, that doesn't matter, since the hardware needs to be there, anyway. Using an IIR in software? What's the point, linear access to memory is preferable and caches are fast, so if the filter's significantly long, it'll probably be fastest to use a FFT-based FIR implementation. In hardware, that'd be a horrible waste of space if the same can be done with an IIR of half the length.

Coming back to the question of phase distortion and filter choices: We know it's very simple to test whether a FIR filter has linear phase: iff its taps are symmetric (hermitian in the complex case), then it has linear phase. That is a strong tool, also for designing linear phase systems (don't confuse with linear systems). Approximating linear phase with an IIR can get pretty complicated / impossible, and getting linear phase with a nonlinear filter – I'm a communications engineer; non-linear filters are pretty much frequency mixers, so we, most of the time, avoid those, if we don't intend to drop phase information anyway.


Non-uniqueness of continuous time sinusoids? I have read the DSP literature and I get why:
$$e^{j(\omega+2 \pi m)n}=e^{j\omega n}$$

But can't we prove the same in continuous domain also in the following way:
$$e^{j\omega t}$$
$$e^{j(\omega+2 \pi m)t}$$
$$\left( e^{j\omega} \cdot e^{2j\pi m}\right) e^t$$
$$e^{j\omega t}$$

If my above reasoning is correct, then does not it mean that analog frequencies separated by $2\pi$ are same? But I know this is not the case.

Please help me in this case.

In discrete-time we have

$$e^{j\omega_0n}=e^{j(\omega_0+2\pi k)n},\quad \forall n$$

because $2\pi kn$ is an integer multiple of $2\pi$ if $k$ and $n$ are integers, and, consequently, $e^{j2\pi kn}=1$. In continuous time you generally have

$$e^{j\omega_0t}\neq e^{j(\omega_0+2\pi k)t}$$

because $2\pi k t$ is (for arbitrary $t$) not an integer.


How up sampling of a signal is done in real world? I know that up sampling by L means padding L-1 zeros between each samples. My question is how it is done practically? Do we use MATLAB or any other software, or there are hardware circuits to do so.

If padding L-1 zeros were the only thing that's done for upsampling, nobody would actually store the end result with zeros that way, because the zeros would be redundant.
In practice, padding with L-1 zeros is always followed by low pass filtering. If that filtering is done with an FIR filter, a dumb implementation would have many coefficients being multiplied by 0. A smart implementation will instead use a polyphase approach where only those coefficients that won't see a guaranteed 0 will be multiplied and summed.
So that partially answers your question: most hardware or software implementations won't even do any padding. They will instead cycle through a set of polyphase filter coefficients, and do something with the result of these successive FIR filter outputs.
Things might be different when doing Matlab or a NumPy implementation: in that case, quick prototyping is often more important than the computational or storage efficiency of the implementation, and rather than using a polyphase filter, people will indeed insert zeros in between samples and convolve the result with the non-polyphase version of the FIR filter.


Calculating IIR Filter gain at given frequency Let's consider an IIR filter with transfer function: $H(z)$.
Given the sampling frequency $F_s$ how can I calculate gain at say $F$ ?
When I was dealing with analog systems when I wanted to calculate gain at say $F=0 Hz$ (DC gain) I would simply put $S=0$ and I would find out what DC gain is.
Now I'm not sure if the same technique applies here.
My inital thought was to put $z= e^{j2\pi}=1$ into the equation but it won't work.
The problem boils down to finding $|{H(z)}|$ which I don't know how to do.
How can I approach this problem?

For continuous-time systems, you obtain the frequency response by evaluating the transfer function $H(s)$ on the imaginary axis $s=j\omega$ (assuming stability). In discrete-time, you get the $\mathcal{Z}$-transform instead of the Laplace transform, and the imaginary axis is replaced by the unit circle. So the frequency response is obtained by evaluating the transfer function $H(z)$ for $z=e^{j\omega}$. Note that for discrete-time systems, $\omega$ is a normalized frequency (in radians):
$$\omega=\frac{2\pi f}{f_s}\tag{1}$$
where $f_s$ is the sampling frequency.
Coming to you example, a frequency of $8000$ Hz with a sampling frequency of $4000$ Hz is the same as DC (Eq. $(1)$ should make that clear). Computing the frequency response at any frequency means to determine the desired $\omega$ from Eq. $(1)$, and then evaluate $H(z)$ with $z=e^{j\omega}$. In the case of DC ($\omega=0$) and Nyquist ($\omega=\pi$) this becomes especially easy: just evaluate $H(1)$ for DC, and $H(-1)$ for Nyquist.
For most other values of $\omega$, the frequency response is most likely complex-valued, so after computing $H(e^{j\omega})$ you just compute the magnitude of that complex number to obtain the gain at the given frequency. You don't need to find a general expression for $|H(z)|$.


Can I see from the PSD if a system is a LTI system? I am studying the meaning of the power spectral density (PSD) and I am wondering whether I can see from the PSD of the input signal and output signal if a system is a LTI system?

Suppose that the input signal to the system has power spectral density (PSD)
$S_i(f)$ while the output has PSD $S_o(f)$. Then, if the system is an LTI
system, then it must be that
$$S_o(f) = |H(f)|^2 S_i(f) \tag{1}$$
where $H(f)$ is the transfer function of the system. So, for a given pair
of PSDs, is there a $H(f)$ such that $(1)$ holds? One obvious case when
$(1)$ cannot hold for any choice of $H(f)$
is when the supports of $S_i(f)$ and $S_o(f)$ are different (in the sense
that there are frequencies $f$ for which $S_o(f) > 0$ while $S_i(f) = 0$).
So in this
case we can say unequivocally that the system is not an LTI system.

What about when the supports of $S_i(f)$ and $S_o(f)$ are identical or
the support of $S_o(f)$ is a subset of the support of $S_i(f)$ (think
bandstop or notch filter)?
Well, in this case, we do have a $|H(f)|^2$ from which we can
deduce $|H(f)|$. Throw in a conjugate-symmetric $\angle H(f)$ and we
have constructed the transfer function $H(f) = |H(f)|e^{\angle H(f)}$
of an LTI system whose output
PSD is $S_o(f)$ when the input PSD is $S_i(f)$. Note that there can
be infinitely many LTI systems that will provide this specific
input-output relationship.

But, the above does not allow us to say that the system that actually
produced the given output PSD from the given input PSD is an LTI
system, let alone one having the specific transfer function that we have
constructed. The actual system that we are looking at might be nonlinear
or time-varying, etc. We cannot say. After all, the
$|H(f)|e^{\angle H(f)}$ that we have come up with
might be a physically unrealizable transfer
function.


Auto-covariance of the product of deterministic and wide-sense stationary signal Anybody give me an advice how to find the auto-covariance of the product of deterministic and wide-sense stationary signal. I couldn't find how to solve this, I have looked and searched the internet, I found just 1 similar topic here so I think someone in here may solve it.

The mean and autocorrelation of $X(t)$ are $m$ and $R$, respectively, and the $X(t)$ process is wide-sense stationary. $g(t)$ is a deterministic function, $Y(t)=X(t)\cdot g(t)$ is defined.

A. Find the mean, autocovariance and autocorrelation functions of the $Y(t)$ process.
B. Is the $Y(t)$ process wide-sense stationary ? Explain your answer.

Let me give a clue
$$ E\{ y \} = E\{ x g \} = E\{ x \} g
$$


A query on Power spectral density (PSD) Say a narrow band signal $n(t)$ has Power Spectral Density (PSD) $S(f)$.

If the signal $n(t)$ got multiplied by $\cos(2\pi Ft)$, then what will the PSD of resulting signal in terms of $S(f)$ be? Will this be $\frac{S(f+F) + S(f-F)}{2}$  ?

Almost, but let's start at the beginning. If the random process $N(t)$ has a power spectrum then it is at least wide-sense stationary (WSS), i.e., its mean and its autocorrelation function do not depend on time. However, the process

$$Z(t)=N(t)\cos(2\pi f_ct)\tag{1}$$

is not stationary, and, consequently, it has no power spectrum (the mean of the random process $Z(t)$ is given by $E[N(t)\cos(2\pi f_ct)] = E[N(t)]\cos(2\pi f_c t)$ which clearly depends on time even if $E[N(t)]$ does not, a similar reasoning can be made for the autocorrelation function).

Luckily, we're usually not interested in the power spectrum of the process $(1)$, but rather in the power spectrum of the process

$$Y(t)=N(t)\cos(2\pi f_ct+\theta)\tag{2}$$

where $\theta$ is a random phase which is independent of $N(t)$ and which is uniformly distributed on $[0,2\pi)$. This random phase reflects the uncertainty of the carrier phase with respect to the process $N(t)$. So adding a random phase is not only a "trick" that works well, but it actually reflects the nature of a modulated random process much better than the process $(1)$, which assumes a known relationship between the carrier phase and the signal.

The process given by $(2)$ is WSS, and its autocorrelation function can be computed as follows:

$$\begin{align}R_{YY}(\tau)&=E[Y(t)Y(t+\tau)]\\&=E[N(t)\cos(2\pi f_ct+\theta)N(t+\tau)\cos(2\pi f_c(t+\tau)+\theta)]\\&=\frac12 E[N(t)N(t+\tau)(\cos(2\pi f_c\tau)+\cos(2\pi f_c(2t+\tau)+2\theta))]\\&=\frac12R_{NN}(\tau)\cos(2\pi f_c\tau)+\frac12 R_{NN}(\tau)\underbrace{E[\cos(2\pi f_c(2t+\tau)+2\theta)]}_{=0}\\&=\frac12R_{NN}(\tau)\cos(2\pi f_c\tau)\tag{3}\end{align}$$

where I've used the independence of $N(t)$ and $\theta$, and where $R_{NN}(\tau)=E[N(t)N(t+\tau)]$.

Finally, the power spectrum of $Y(t)$ is given by the Fourier transform of $(3)$:

$$S_{YY}(f)=\frac14\left[S_{NN}(f-f_c)+S_{NN}(f+f_c)\right]\tag{4}$$

where $S_{NN}(f)$ is the power spectrum of $N(t)$.

In sum, if we add a random phase to the carrier then the modulated process is also WSS (if the baseband process is WSS), and its power spectrum is given by $(4)$.


Noise power transient at output of RC filter Consider a simple first-order RC filter with AWGN input noise with single-sided density of No. We know the output noise power is No/(4RC). What is rarely noted is that this is true after steady-state has been reached. I simulated the noise from such a filter starting from rest at t=0 and found that the ensemble average of the instantaneous noise power seems to increase from 0 to No/(4RC) exponentially with a time constant of RC/2, not the circuit time constant of RC.
Can anyone derive the equation for this?

Just to clarify this is input noise that is switched on by a unit-step at t=0 at the input to the filter.
So I am guessing the usual analysis based on the assumption of stationarity probably does not apply.

Thanks,
John

Consider a simple RC filter whose step response is $1-e^{-t/RC}$ for $t > 0$,
and impulse response is
$$h(t) = \begin{cases} \frac{1}{RC}e^{-t/RC}, & t \ge 0,\\ 0, &t < 0.\end{cases}$$
Then, when the input is white noise, the steady-state output variance is
$$\sigma^2 = \frac{N_0}{2} \int_{-\infty}^\infty h^2(t)\,\mathrm dt
= \frac{N_0}{2}\int_0^\infty \frac{1}{(RC)^2}e^{-2t/RC}\,\mathrm dt = \frac{N_0}{4RC}.$$
However, if the input noise is also $0$ for $t < 0$, then at time $T$, the
noise variance at the output is
$$\begin{align}
\sigma_T^2 &= E\left[\int_0^T N(t)\frac{1}{RC}e^{-(T-t)/RC}\,\mathrm dt
\int_0^T N(s)\frac{1}{RC}e^{-(T-s)/RC}\,\mathrm ds\right]\\
&= \int_0^T\int_0^T E[N(t)N(s)]\frac{1}{(RC)^2}e^{-(2T-t-s)/RC}
\,\mathrm dt\,\mathrm ds\\
&= \int_0^T\int_0^T \frac{N_0}{2}\delta(t-s)\frac{1}{(RC)^2}e^{-(2T-t-s)/RC}
\,\mathrm dt\,\mathrm ds\\
&= \int_0^T \frac{N_0}{2}\frac{1}{(RC)^2}e^{-(2T-2s)/RC}
\,\mathrm ds\\
&= \frac{N_0}{4RC}e^{-2T/RC}\int_0^T \frac{2}{RC}e^{2s/RC}\,\mathrm ds\\
&= \frac{N_0}{4RC}e^{-2T/RC}[e^{2T/RC}-1]\\
&= \left.\left.\frac{N_0}{4RC}\right[1 - e^{-2T/RC}\right]
\end{align}$$
showing that $\sigma_T^2$ increases to its asymptotic value
$\displaystyle \sigma^2 = \frac{N_0}{4RC}$ more rapidly than the unit step
response increases to its asymptotic value $1$.

The time constant of the RC circuit is $RC$, not $1/RC$ as the OP
claims, and the time constant of the increase in the noise variance is
$RC/2$, half that of the RC circuit. The smaller the time constant,
the more rapid the decay to $0$  (or rise to the asymptotic value).


When do we have a minimum phase couple I am a DSP newbie and I am really confused after reading Chapter two of

"An Introduction to Digital Signal Processing by J. H. Karl".

On page 32, first paragraph, the book defines a minimum phase couplet as one with zeros OUTSIDE the unit circle of the z-plane. But I have read everywhere else that a minimum phase system has all zeros and poles INSIDE the unit circle.

Is this a contradiction? When do we have a  minimum phase couplet?

The book uses the convention that $Z$ represents a unit delay (this is explained on this page). The more common convention is to denote the unit delay by $z^{-1}$. Consequently, minimum-phase systems have all their zeros inside the unit circle of the $z$-plane, but outside the unit circle of the $Z$-plane.


window function between program and math I'm newbie on signal processing.
I know when we program python with window function. (just signal multiply with window)
According to Harris paper

I know we change w with different window(blackman).
However in python fft(signal*window(blackman)).
Is there still window(rectangular) function when we do signal capture(fft) in programming?
I mean that dft(signal * blackman * rectangular).

The fact that you are using a time limited sample of the waveform is the rectangular window that is referenced. Since the fft operates on a time limited sequence, it is mathematically identical to a infinitely long duration of that sequence multiplied the rectangular window.


How to create a test to measure end to end audio delay Lets say I have a program like Skype and I call to my friend. I want to create a test to figure out what is the average delay from the audio, from my mic to his speakers.

I want to measure:

delay in one instance of time
delay over longer period (60 s)

Assumptions that i can make:

at least 8000 samples per second.
devices are in the same room and network propagation delay is negligible.
external mic and photo-detector/camera should be used.

I recently ecountered this problem too (about a month ago). This is what I found..

Network delay measuring protocol (NDMP) - little useful free-ware program. I used the same program to measure other VoIP call delay.

You can download it from here: http://www.unforgettable.dk/42.zip
ReadMe file is included in .zip file.


OFDM - why do we need IFFT? This question is really related to - is there any need for IFFT in the first place?

If we want to send a bunch of complex-numbers (vectors), then most OFDM discussions involve presenting those numbers to an IFFT block, and the output of the IFFT block will then have a bunch of complex numbers for which to transmit in a queued fashion through a QAM block.

Can I ask ----- is the IFFT block really necessary? Instead of using the IFFT block, is it possible to just place the original set of complex numbers (vectors) on the queue....and then transmit those numbers (real and imaginary components) - one at a time - through a QAM block?

The signal (without using IFFT) is still going to be OFDM, right?

Thanks for any comments/advice!

The “O” in OFDM stands for orthogonal.  The IFFT takes a complex number (1 bin of the input) and turns it into samples of a sinusoid of a frequency that, over a certain length (of time), is orthogonal to (and thus, under ideal linear conditions, won’t interfere with) any other frequency subcarrier output by the IFFT.

If you just transmit the original bits, the impulse response of the channel can smear the bits together, interfering with the adjacent bits, and leave the decoder/demodulator the interesting job of unscrambling this mess.  Depending on the channel response and how it changes over time, this may or may not be more difficult to do.


How to build a phase shifter with arbitrary phase shift Fred, a DSP engineer, goes to his favorite DSP store to do some shopping.

Fred: Hi, I'd like to buy a phase shifter.

Shop assistant: Hmm, what exactly do you mean?

Fred: Well, you know, if you put in a sinusoid like $x(t)=\sin(\omega_0t)$ you get $y(t)=\sin(\omega_0t-\theta)$ at the output, for any $\omega_0$. And of course, $\theta$ must be adjustable.

Shop assistant: Oh, I see. Sorry, no, we don't have those. But I remember other guys needing the same thing, and they always buy a Hilbert transformer, a couple of multipliers, and an adder, and they somehow connect all these things together to make an adjustable phase shifter.

Fred: Oh yes, right!

Fred pretends to understand what the guy is talking about. Of course he has no idea how to do that. He buys everything the guy said he needed, and thinks by himself that he might figure it out at home, or, everything else failing, he could ask it at DSP.SE.

How can Fred build a phase shifter with adjustable phase shift $\theta$ using the components he got at the store?

Nice question! It uses one of my favorite trig identities (which can also be used to show that quadrature modulation is actually simultaneous amplitude and phase modulation).

   The Hilbert transform of $\sin(2\pi f_0t)$ is $-\cos(2\pi f_0t)$. Also, $$\sin(2\pi f_0t+\theta)=a\sin(2\pi f_0t)+b\cos(2\pi f_0t)$$ (constrained to $a^2+b^2=1$), with $\theta=\text{atan2}(b,a)$. This suggest one possible approach. Say Fred needs $\theta=2.1$ radians. He calculates $\tan(2.1)\approx-1.71$. Then, he needs to find $a$ and $b$ such that $a^2+b^2=1$ and $b/a=-1.71$, with $a<0$ and $b>0$, which is a simple algebra problem. Set $a_0=-1$, $b_0=1.71$, $n=\sqrt{a_0^2+b_0^2}$, $a=a_0/n$, and $b=b_0/n$. Then, Fred can easily generate a sine with the desired phase by using a Hilbert transformer, two multipliers, two DC sources (one set to $a$ volts and the other to $-b$ volts, to take care of the sign of the cosine), and one adder.

The impulse response of the system described above is given by:

   $a\delta(t)+\frac{b}{\pi t}$

Block diagram:


Getting the $\mathcal Z$-transform of the moving average summation I'm getting started with signal processing, and I'm stuck on a problem that asks me to take the $\mathcal Z$-transform of the following causal, DT-LTI system:
$$\sum_{k=0}^M b_k  x[n-k] $$

I'm not really sure how to take the $\mathcal Z$-transform with the summation in there. Can anyone point me in the right direction?

You can use the linearity property of the $\mathcal Z$-transform:

$$a_1x_1[n] + a_2x_2[n] \overset{z}\longleftrightarrow a_1X_1(z) + a_2X_2(z)$$

together with the time shifting property:

$$x[n-k] \overset{z}\longleftrightarrow z^{-k}X(z)$$

Expand the summation like below and apply the two properties.

$$\sum_{k=0}^M b_k  x[n-k]=b_0x[n]+b_1x[n-1]+\ldots+b_Mx[n-M]$$


aliasing folding phenomena in frequency domain matlab Hello i am trying to create the folding phenomena of undersampling in matlab,
When i undersample the sampling frequency is 135 less than the Nyquist frequency for 70 Hz signal
ou will see that it is shifted back by the amount of this new sampling frequency (105-70=35) Hz.
as shown in the matlab plot bellow.
I want to create the folding effect shown in diagram bellow (shifted copy).
Where did i go wrong.
This is my code which producedthe plot in the first post.
my code is shown as one line for some reason.

f1=10;
f2=30;
f3=70;
% twice the sampling rate
Fs=1.5*70; % sampling frequency is a bit above 2 times to get all the peaks.
Ts=1/Fs;
Tn=0:Ts:1;
fft_L=length(Tn);
y4_samples=10*sin(2*pi*f1*Tn)+10*sin(2*pi*f2*Tn)+10*sin(2*pi*f3*Tn);
%stem(Tn_new,y4_samples);
ff=fft(y4_samples);
ff1 = abs(ff/fft_L);
fft2 = ff1(1:floor(fft_L/2)+1);
fft2(2:end) = 2*fft2(2:end);
f = Fs*(0:fft_L/2)/fft_L;
plot(f, fft2)

Your plot is correct. You are sampling three waves $f_1 = 10\textrm{Hz}, f_2 = 30\textrm{Hz}$, and $f_3 = 70\textrm{Hz}$, with a sample frequency of $F_s = 1.5\times70\textrm{Hz} = 105\textrm{Hz}$. This means that your Nyquist frequency is $F_N = F_s/2 = 52.5\textrm{Hz}$, and corresponds to the maximum value of your frequency axis.

As such, the signals $f_1$ and $f_2$ will be correctly sampled without problems, but the signal $f_3$ is higher than the Nyquist and so should appear as an alias (folded back into the sampled frequency range) at the frequency $f_{\textrm{alias}} = F_N-(f_3-F_N) = F_s-f_3 = 35\textrm{Hz}$. This is what you obtained.

If you want to create a plot similar to the one in the figure, then you could try to first sample your waves with a much higher sampling rate, in order to fully capture them all faithfully. Then, on the same plot, display what happens if you undersample with $F_N<f_3$.


Why is random noise amplitude inversely proporional to the frequency? http://astarte.csustan.edu/~tom/SFI-CSSS/info-theory/info-lec.pdf compares a bacteria DNA vs a random string:

  My first step was to generate for myself a "random genome" of
  comparable size to compare things with. In this case, I simply used
  the Unix `random' function to generate a file containing a random
  sequence of about 4 million A, C, G, T. In the actual genome, these
  letters stand for the nucleotides adenine, cytosine, guanine, and
  thymine.

  We want to find "interesting" sections (and features) of a genome. As a
  starting place, we can slide a "window" over the genome, and estimate
  the entropy within the window. The plot below shows the entropy
  estimates for the E. coli genome, within a window of size 6561 (= 3^8
  ). The window is slid in steps of size 81 (= 3^4 ). This results in
  57,194 values, one for each placement of the window. For comparison,
  the values for a "random" genome are also shown.

  At this level, we can make the simple observation that the actual
  genome values are quite different from the comparative random string.
  The values for E. coli range from about 5.8 to about
  5.96, while the random values are clustered quite closely above 5.99 (the maximum possible is log2 (64) = 6).

  There are various "data massage" techniques we could use. For example,
  we could take the fourier transform of the entropy estimates, and
  explore that. Below is an example of such a fourier transform. Notice
  that it has some interesting "periodic" features which might be worth
  exploring. It is also interesting to note that the fourier

  transform of the entropy of a "random" genome has the shape of
  approximately $1/f = 1/f^1$ (not unexpected), whereas the E. coli
  data are closer to $1/f^{1.5}$

I must draw you attention that the fourier transform is applied to the entropy rather than to the random signal itself. Why is it expected that resulting spectrum is pink rather than absolutely random? I do not understand why it is expected that the random frequency amplitude is inverse to the frequency, why not $amplitude=1/freq^{1.5}$? The higher low frequencies mean that the next number has higher probability to be closer to the current one and random is not quite random, right?

I have also noted that, in the end of the book, author says that "random" is naturally band-limited by the power law, $1/f^\alpha,$  and, in SETI project, we better watch after the white noise because optimal use of the channel is a sign of intelligent life.

This looks like a pink noise model, where the power spectral density of the noise is proportional to $\frac{1}{f}$. This model is appropriate for some real-world phenomena; I can't speak to whether it's the right choice for this particular application.

Your interpretation of the meaning of this sort of spectrum is correct, although the noise is still considered random. Instead, we consider consecutive noise samples to be correlated with one another (in contrast to white noise, where all samples are jointly uncorrelated). That doesn't mean that the signal isn't random.


Determining Parameters of S Function for Image Global Contrast Enhancement I know that S-function (membership function from fuzzy set theory) is a commonly used function for image enhancement. I am reading a paper titled - A Novel Fuzzy Image Enhancement Using S-Shaped Membership Function. The S-shaped membership function is defined as follows:

$$\begin{eqnarray*}
S(x;a,b,c)=\begin{cases}
0, & x\leq a\\
\frac{(x-a)^2}{(b-a)(c-a)}, & a<x\leq b\\
1-\frac{(x-a)^2}{(c-b)(c-a)} & b\leq x\leq c\\
1, & x>c.
\end{cases}
\end{eqnarray*}$$

where the parameters are to be found. In the above mentioned paper, it given that the values of a, b, c are entropy(I), median(I) and (max(I)+mean(I))/2 respectively. But there's no explanation for how and why these values were chosen. I tried to look for it in many papers referred by the above paper but they are of no use. Can someone give me some idea please? Thank you in advance.

Okay, I am convinced that there is a typo in Eq. (4) of the cited paper.  It should be:

$$ S(x;a,b,c)=\begin{cases}
0 & x\leq a\\
\frac{(x-a)^2}{(b-a)(c-a)} & a<x\leq b\\
1-\frac{(x-c)^2}{(c-b)(c-a)} & b\leq x\leq c\\
1 & x>c
\end{cases} $$

This sigmoid function or smoothstep function is continuous everywhere and the first derivative is continuous everywhere, but the second derivative is discontinuous when $x=a$ or $x=b$ or $x=c$.

You must have $a < b < c$.

$a$ is defined where $S(x)$ begins to transition from 0 upward.  and $c$ is defined for when the transition is over and $S(x)=1$.  $b$, which is in between $a$ and $c$ is the point of inflection for when the curvature changes from a concave curvature to a convex curvature (or maybe it's the other way around).

There are better sigmoid or smoothstep functions than this that the author of the paper uses.  I dunno why he/she is so enamored with this definition.


Finding the Width of a Point Spread Function (PSF) I have a binary image that has been convolved with the following PSF:

$$\mathrm{PSF}(x,y)=\cosh^{-2}\left(\frac{\sqrt{x^2+y^2}}{w}\right)$$ where $w$ is unknown.

The image binary image generally consists of contiguous black and white features that are larger than $w$. What is the most efficient way to estimate $w$? I typically don't have a problem to look at the image and provide an initial guess for $w$, and my goal is to get it accurately.

One simple way is to identify zeros in the Frequency Domain.
This was one of the legitimate method for Blind Deconvolution few decades ago.

Draw the connection between the zeros in frequency domain of the Convolution Kernel to the parameter $ w $.
Learn the zeros in frequency domain of the given image.
Find the $ w $ which explains most of the zeros of the image.

It's pretty simple and for a parameterized family of PSF's can be quite efficient.


Radar Pulse Compression Gain (PCR) Assume a monostatic radar is transmitting a pulsed waveform. The pulse has been modulated by a Barker code with$\ N$ chips or bits. The waveform chip rate is $\tau_c $, the radar receiver's detection bandwidth is $\beta_c = 1/\tau_c$, and the waveform's pulsewidth is $\ N*\tau_c$.

I've heard conflicting opinions about a signal processing gain,$\ P_g$, due to compressing the pulse that is proportional to the waveform's time-bandwidth (or pulse compression ratio (PCR)) product (or equivalently the product of the uncompressed pulsewidth $\tau_p$ and the compressed bandwidth $\beta_c$) of the pulse: $\ P_g = \tau_p*\beta_c=N\tau_c*\beta_c$

Assuming that you start with the power form of the radar range equation:

Is it appropriate to add a$\ P_g$ term to the power form of the SNR radar range equation?

(Pulse compression allows a radar to transmit a longer pulse to get more energy on the target while simultaneously achieving the range resolution of a much shorter pulse.)

FWIW, I have the same question about the energy form of the radar range SNR equation:

My question is similar to this older question but more specific

The answer is yes but one has to specify $B_n$ properly to avoide possible confusions. In case if one uses a pulse compression, the bandwidth through which the receiver collects the noise will normally be $B_n = \beta_c$. Then, the "new" signal-to-noise ratio should be written as:

$SNR = \dfrac{P_TG_TG_R\lambda^2\sigma{P_g}}{(4\pi)^3R^4(kT_{sys}\beta_c)} = \dfrac{P_TG_TG_R\lambda^2\sigma{\tau_p\beta_c}}{(4\pi)^3R^4(kT_{sys}\beta_c)} = \dfrac{P_TG_TG_R\lambda^2\sigma{\tau_p}}{(4\pi)^3R^4(kT_{sys})}$

As you see, the signal-to-noise ratio remained the same as for the pulse without modulation (for which $B_n = 1/\tau_p$).

So, the pulse compression helps to increase the range resolution but the signal-to-noise ratio remains the same as for the unmodulated pulse. This holds true if you do not change $\tau_p$ and the peak power $P_T$ when you modulate your pulse.

The confusion (and as you mentioned "conflicting opinions") comes from the fact that it is quite often forgotten that if you increase the bandwidth of your pulse then you also increase the bandwidth through which the receiver collects the noise since radar's receiver is matched to the pulse and therefore has a bandwidth $\beta_c$.

Your last equation is correct. But since you put $F_n$ - the noise figure- then it should be specified that this is the signal-to-noise raito at the output of the receiver that has a noise figure $F_n$.


Expanding piecewise polynomial using Daubechies wavelet What is the best Daubechies wavelet (i.e. the number of vanishing moment) to expand a signal $\boldsymbol{x} \in \mathbb{R}^n$?
$\boldsymbol{x}$ consists of $m$ pieces of polynomial with $d$ degree. The criterium is to make the DWT signal as sparse as possible.

Update:
The goal of sparsifying the signal in wavelet domain is denoising. Let $\boldsymbol{W}$ denote DWT.
$$
\boldsymbol{y} = \boldsymbol{Wx}
$$
Apply a soft-thresholding to $\boldsymbol{y}$,
$$
\hat{\boldsymbol{y}} = \text{sign}(\boldsymbol{y})(\max(\boldsymbol{y}-\lambda,0))
$$
Choose $\lambda=\sqrt{2\log n}$ according to Ideal spatial adaptation by wavelet shrinkage.
The sparsity is defined by $\| \hat{\boldsymbol{y}} \|_0$.

According to your formula, you also apply soft-thresholding to the approximation coefficients, which is not standard. Aside, your operator $W$ does not seem to specify the number of wavelet levels used. Finally, your class of signals does not seem to address the regularity at piecewise junctions.

I believe in this case very unlikely that in a discrete implementation, without further relation between $m$ and $n$, you can find, theoretically,  a best wavelet in all cases, because $\|\hat{y}\|_0$ is a quite sensitive index (and it is not a norm).  Of course, a Daubechies with $d$ vanishing moments would seem appropriate.

But since the DWT is quite fast, in your context, you could just find a "generally best wavelet",  by simulating random signals, and iterating over some levels, and each Daubechies wavelet with moments in $[d-2,\ldots,d+2]$ for instance.


After applying a Hamming window to your data I get very similar autocorrelation values:

R =
   1.6765e+00
   1.5828e+00
   1.3697e+00
   1.0688e+00
   7.4186e-01
   4.6660e-01
   2.5672e-01
   1.2121e-01
   4.9779e-02
   1.6453e-02
   3.1701e-03

You can check your LPC filter coefficients (not the reflection coefficients) by computing the matrix-vector product $\mathbf{R}\mathbf{a}$ where $\mathbf{a}$ is the vector of LPC coefficients and $\mathbf{R}$ is a Toeplitz matrix with the autocorrelation values:

$$\mathbf{R}=\begin{bmatrix}R(0)&R(1)&\ldots && R(p)\\
R(1)&R(0)&R(1)&\ldots&R(p-1)\\
\vdots& & \ddots & &\vdots\\
R(p) & R(p-1) & \ldots && R(0)\end{bmatrix}$$

This product must result in a vector with the first element being non-zero, and all the other elements must be zero (up to numerical accuracy).

Using the Octave function [lpc,v,k]=levinson(R(1:11)); with the autocorrelation function given above, I get the following LPC filter coefficients:

lpc =

   1.000000
  -1.370890
  -0.386146
   0.813691
   0.992372
  -0.949576
  -0.642884
   0.427407
   0.515102
  -0.310107
  -0.039544

And these are the corresponding reflection coefficients:

k =

  -0.9440538
   0.6826730
   0.2730461
  -0.1825513
  -0.4886926
   0.1882461
   0.3690086
  -0.0056205
  -0.3648883
  -0.0395440

So we get very similar values from which you can conclude that your code for computing the reflection coefficients from the autocorrelation function is correct.


What's the physical meanning and principle of applying anisotropic diffusion in transform domain？ I am reading several papers applying anisotropic diffusion in the transform domain like Shearlet, curvelet. All of them have got excellent results.

The physical meaning and principle of applying anisotropic diffusion in space domain is quiet straightforward -- the gradient.

I am confused as to the physical meaning and principle of applying anisotropic diffusion in the transform domain.

Papers, for example:

Shearlet-Based Total Variation Diffusion for Denoising
Glenn R. Easley, Member, IEEE, Demetrio Labate, and Flavia Colonna
IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 18, NO. 2, FEBRUARY 2009

Physically, there is usually not much implied meaning (but this may depend on a case-by-case basis). Transforms like shearlets are actually used to embed the data into anisotropic or edge-aware transform domains.

Mathematically, it may make more sense. Isotropic diffusion tends to create blur (smear) in whatever domain it is applied. Thus, if you want to obtain a result that is compact in the considered domain, be it transformed or not, applying non-isotropic diffusion is very tempting.

Another way to see it mathematically, is to understand that isotropic diffusion corresponds to solving a least square problem that will create many small non-zero coefficients by trying to spread the error everywhere. However, a helpful and widely available assumption about signal is their sparsity in the transform domain, i.e., they have few non-zero coefficients in a suitably chosen transform domain that concentrate the energy. This corresponds to solving problems where the exponent of the norm of the prior term is less than 2 (typically 1) and thus the use of anisotropic diffusion to find a solution.


What should I do with FFT result array from microphone input? I currently want to create a small pitch detector and I have to say that I am very new in this topic of signal processing.
At the moment I implemented the FFt algorithm and I got an array of window size (= 1024) out of it. But what should I do to define for example C4 or other music tones.

I know that FFT is not the best way to do pitch detection, but I think it is a good starting point to detect some single tones.

Is the best approach to sing the wanted tone than to find the peaks in this array and map that to the tone or do I need to do some other calculations first?

I hope you can help me out with some information.

Thank you and best wishes,

MusicMagician

EDIT:
I know that here are many similar posts, but I want to understand what exactly the FFT does with my frequency and how I should work with that result. I think that for some guitar tuners the FFT is used to detect the pitch. So I wanted to start with that simple approach. I hope you can understand me.

Musical tone or pitch is different from FFT magnitude peaks.  A musical waveform is rarely a pure sinusoid, integer periodic in aperture.

Please, instead research pitch detection/estimation algorithms (several listed in wikipedia) that do a lot more than simply look at FFT results to determine the tone or pitch of the human singing voice.


Phase Modulation and Intermodulation Anomaly ? (Spectrum plots inside - **Updated**) I have a MATLAB script for modelling a two tone signal (a sum of two sinusoids) going through a non linear transfer function (such as an amplifier). The amplitude of a signal is amplified non linearly, and also the phase is modified depending on the input power of the signal.

The two sinusoids used are at 20 Hz and 21 Hz.

With just amplitude modulation, the output spectrum looks as expected (first two plots) Its non linear so 3rd order and fifth order intermods can be seen around the frequency of the two sinusoids, aswell as harmonics. I know the amplitude modulation is correct because I obtain the correct third order intercept for the amplifier model. Upper graphs are in volts, lower graphs in dB Watts

However, when phase is adjusted depending on the input signal power, the spectrum looks like this... and my fifth order intermods are lost (again upper graph in volts, lower graph in dBW)

I have an idea but I am not confident with the explanation of this and was hoping someone with a trained eye can assist or offer advice

How its coded --> My script looks at each point on the input signal, reads its amplitude and adds a phase adjustment to it during reconstruction of the signal, as shown in the phase transfer plot below (which shows the effect the equipment will have on any signal going through it). So if the amplitude (input power) of the signal is -70 dBW then a phase of 0.01 radians is added to the reconstructed signal at that exact point... then it looks at the next sample point and adds a phase to that. So if the amplitude of the wave is at -70dBW then the reconstructed waveform is cos(2*pi*f + 0.01).  Adding 0.01 in phase is an I and Q modulation when cos(2*pi*f + 0.01) is written with the trig identity to get the IQ form.

Confident I have not made any mistakes with unit conversion (again because third order intercept is correct).

Look at the input signal, output signal and the phase deviation added to the signal in time domain... it looks good

A couple of points to think about:

1) Looking at the scales of your FFT plots, it appears that the fundamental is around -25 dBW, whereas the 5th harmonic is around -125 dBW. This is a very large difference, implying the system is only weakly nonlinear. In practice you'd have problems detecting the 5th harmonic even without the phase modification just due to thermal noise. If you change your scales to have a dynamic range of say 100 dB, you would see the two plots look nearly identical.

2) The phenomenon you are modelling looks to be AM-PM conversion. Under this interpretation, it can be thought of as "adding" some phase noise to the output signal. (The tails of the phase noise spectrum would droop down and mask weaker signals.) This could be an explanation for why the result looks like a noise floor.

3) A sanity check could be-- what does the spectrum look like when the peak phase shift is halved? Or cut to 1/4th, etc? Does the noise floor drop by some commensurate amount?


What are the units of Es/N0? I want to make sure that my units are all correct. I've been using the formula

$$
\frac{P_S}{P_N} = \frac{E_S}{N_0} \frac{R_S}{B}
$$

The dimensions that I'm confident of are

$P_S$, $P_N$ : Watts (Joules/sec)
$N_0$ : Joules/(sec*Hz)
$B$ : Hz

The last two are then (I think)

$R_S$ : symbols/sec
$E_S$ : Joules/symbol

If those are all correct, and given that

$$
\frac{E_S}{N_0}  = \frac{P_S}{P_N} \frac{B}{R_S}
$$

it seems like the units of $\frac{E_S}{N_0}$ are

$$
\frac{\text{Hz}}{\text{symbols/sec}} = \frac{1}{\text{symbol}}
$$

I was under the impression that it should be dimensionless. Is this correct?

The "symbols per second" (or, more strictly, the baud) unit is a bit weird. It is really a unit that is used for human convenience and not for mathematical rigor. Just like "samples per second", bauds can be interpreted as hertz and still make sense, since they measure the "number of times something happens per second".

So, in your calculation, the units of $B/R_S$ are $\text{Hz}/\text{Hz}$, making it unitless.

I know this is somewhat unsatisfying (at least, it is to me), but this is the best explanation I can give.

(All your unit definitions at the top of the question are correct, BTW).


Farrow fractional delay filter - Range of delay? I am implementing a fractional delay for resampling purposes, using the Farrow structure to enable continuously variable delays. For now I am simply using windowed sinc functions as my low-pass filters. My method was originally as follows:

Create a series of $K$ $N$-sample windowed sinc functions covering the delay range: $N/2-1$ to $N/2$, and all the fractional delay values ($d$) in between.
For each sample across all the filters, create a $K-1^{th}$ order polynomial in d to replace the coefficients.
Stick the poynomial coefficients into a Farrow structure filter, and Bob's your uncle.

However, on reading some papers this for example, it became apparent that I need to consider the "overall filter"; all the separate polyphase branch filters from the farrow structure lumped together. Now, taking this into account, if I crate my individual filters as above, ranging from $N/2-1$ to $N/2$, then I have a problem, as shown in the image. The "overall filter" is not a nice lowpass function anymore. Essentially it the two filters at $N/2-1$ and $N/2$ are the same, just delayed 1 sample, so we get the reptition of values shown in the plot.

My attempt at solving this is instead to create my initial filters over the range N/2-1+0.5/K : N/2-0.5/K. This way, they are still evenly distributed in terms of delay and now the "overall filter" looks great:

BUT... now my delays do not cover a full sample of fractional delay, but instead $1-1/K$ samples. The only way to approach a full 1-sample range is to increase $K$, the number of filters created at the start, and this is not ideal as it just increases computation and I am VERY restricted in that sense.

I also tried a second approach, creating $K$ original filters, then approximating with an $L^{th}$ order polynomial, where $L < K$. This allows me to create a large ($K$) number of filters to start with, so I can get very close to the full N/2-1 : N/2 range, then approximate with a low order polynomial. This is OK, but it does introduce some unwanted artefacts into the "overall filter" response.

So, my questions is this: Is there a way to have it all; a low-order farrow fractional delay filter, that covers the full range of a sample delay, without introducing any nasties into the frequency response of the "overall filter"?

Any thoughts welcome :)

P.S. I know the sinc functions in the plots are rectangular windowed and thus have a high sidelobes, but this is not my problem.

Do you think this link could help you?

http://users.spa.aalto.fi/vpv/publications/vesan_vaitos/ch3_pt1_fir.pdf

I'm not sure I understand 100% of your problem, but IIRC when implementing fractional delays with say a FIR filter of 13 coefficients, you should stick with a group delay between 6 and 7.


MFCC feature vector from wav file I am trying to implement a spoken language identifier from audio files, using Neural Network. I need to generate one feature vector for each audio file.

From what I have read the best features (for my purpose) to extract from the a .wav audio file are the MFCC.

Does anyone know of a Python code that does such a thing?

Have a look at these two python libraries that provide a number of audio features easily from WAV files, including MFCC.

Librosa: MFCC docs, github
Madmom: MFCC docs, github

Good luck!


confusion about poles and zeros of Lead compensator? I was reading about lead compensator from a website but i am not able to understand how they have extracted zero and pole from transfer function of lead compensator as shown in attached photo?

or they have extracted/calculated wrong values of zeros and poles?

You're right, the given pole and zero are wrong. They should be

$$s_0=-\frac{1}{\tau}$$

and

$$s_{\infty}=-\frac{1}{\beta\tau}$$

because for $s=s_0$ the numerator becomes zero, and for $s=s_{\infty}$ the denominator becomes zero.


position of spike in the spectrum of a sin wave I have a sin wave of 16Hz, sampled at 2048 samples persecond. I performed the FFT of the sequence and did a spectrum plot. My plot gave me two spikes at both ends of the graph. But i was told the spike is meant to be at 16. how exactly am i supposed to do the plotting in order to get this exact spike?

A sinusoid is the sum of two complex exponentials ("negative" and "positive" frequency) so that their imaginary parts cancel out each other:

$$cos(f \cdot t) = \frac{1}{2} \left( exp(f \cdot t \cdot i) + exp(-f \cdot t \cdot i) \right)$$

The Fourier Transform decomposes your input into a weighted sum of complex exponentials. Therefore, you get to see two pulses in the spectrum. One at -f and another at +f, basically. Typically, the "negative half" of the spectrum is located after the positive half because in discrete signals there is no difference between -f Hz and fs-f Hz where fs is the sampling frequency. So, your FFT result's ordering is 0 Hz to almost fs Hz where you can think of the second half as the negative portion of the spectrum.

Since real-valued signals have an imaginary part of zero, you'll notice that the Fourier Transform of real-valued signals always exhibits a symmetry between negative and positive frequencies with a certain kind of phase relationship (same amplitude and negated phase) so that the complex exponentials' imaginary parts cancel between negative and positive frequency.

Complex-valued signals are more general in that their spectra don't have to be symmetric around 0 Hz. The imaginary part helps you distinguish between positive and negative frequencies:

$$exp(f \cdot t \cdot i) = cos(f \cdot t) + i \cdot sin(f \cdot t)$$

As you can see above, the sign of f does not change the real part. But it does affect the sign of the imaginary part.


Receiver figure of merit So I am solving receiver figure of merit. I have gain in dB and noise temperature in K.

I know

[G/T] = [G] - [T] dB/K

Here should I need to take T in dB or just as K?

The relationship [G/T] = [G]-[T] dB/K is incorrect.

G/T when in units are already given in dB/K is simply the gain in dB divided by the temperature in K.


How to find a noisy line in a picture Hi I would like to find a line in a noisy image as showed bellow. I tried gabor filter and cany edge detection, but it didn't produce anything usefull. I am looking for lines of specific orientation and I don't want to segment them, just to find if the line exist in the picture.

Example pictures with lines

Example pictures without lines

What about an algorithm in image domain?

For every pixel, check into the direction you desire to know (up-down in your images above, I guess) if there is a pixel with "the same" value (i.e. a value inside an acceptance band, of course. Let's say +/- 10%) as the pixel under investigation. If so, go further into that direction. Count every pixel that you find. If the number of counts exceeds a certain threshold (i.e. 1/4 FOV or so in your case), then there is a line.

I would first try to use something simple like that instead of using more complicated algorithms that might be suitable for different kinds of images. Yours pretty much look like MRI phase images that can be nasty ;) )


What is the most precise frequency analysis spectrography algorythm? Everyone uses Fast Fourier Transform, which is fast at the detriment of precision. The input audio has sample accuracy and the FFT has 1/64 sample accuracy.

What algorithms can output high resolution frequency analysis of audio that varies over single samples, and that can divide input audio into i.e. 1024 bands of different frequency sines, each of the 1024 tracks is 22k sample rate?

You don't have a loss in time precision when using FFTs because the FFT is fast. The FFT is just a fast algorithm for implementing the discrete Fourier transform (DFT), nothing more. Instead, there is an inherent tradeoff in time and frequency resolution due to the Heisenberg uncertainty principle. While its statement is explicitly focused at quantum mechanics, the same underlying principle remains true: the more precisely you know the frequency of a signal, the less able you are to localize it in time.

With that said, there are another class of techniques known as bilinear time-frequency distributions that are appropriate for some applications. One example is the Wigner-Ville distribution. In short, these techniques can provide simultaneously high resolution in time and frequency. The cost, however, is the presence of spurious features in their resulting outputs. There do exist modified versions of these distributions that can reduce the magnitude of the artifacts, however.


C/C++ library for digital modulation Is there a library written in C/C++ that offers functions like QPSK modulation?

I have created a C++ library that I haven't officially released yet that could be used for modulation.  It is called NimbleDSP.  It doesn't have a function or method called "modulation", but it does have interpolation and supports complex numbers, so it could easily be used to do modulation.

Liquid-dsp appears to be a mature SDR library, and it is used by the government framework RED HAWK that was recently open sourced.  My main problem with the library is that it uses the GPL license (may not be a problem for you) and the interface seems complicated.  Having said that, I looked at it fairly briefly so I could be wrong about that.

There is also SPUC (Signal Processing Using C++) and Aquila.  SPUC is also GPL I believe, but looks to be simpler than liquid-dsp while still being very comprehensive.  Aquila has a very permissive license but has some odd holes in its functionality.


Dithering adds harmonic distortion after filtering and upsampling I have an issue which I wasn't able to solve for a few days. I have implemented a polyphase FIR interpolation filter within an FPGA which does interpolate by a factor of 8 or 16 times. I do accept an audio word of up to 32 bits and coefficients are quantized on 35 bits resulting in MAC units of 32x35 and an accumulator with width of 67 bits.

The filter itself works flawlessy. There aren't any issues in general. However, I did add a dithering algorithm which uses LFSR to generate a random number two times in a row which is then added to the output from the filter (67 bits wide) in order to reduce bit-depth from 67 bits to lets say 18 bits. LFSR does generate a uniform distribution, so adding two pseudorandom numbers generated by it results in a triangular distribution (TPDF dithering).

For the sake of simplicity I did implement following equation within VHDL:

output_lch_sample_dithered_int <= std_logic_vector(signed(output_fifo_lch_sample_int) + signed('1' & dither_rnd1_int(47 downto 0)) + signed('0' & dither_rnd0_int(47 downto 0)));
output_rch_sample_dithered_int <= std_logic_vector(signed(output_fifo_rch_sample_int) + signed('1' & dither_rnd1_int(47 downto 0)) + signed('0' & dither_rnd0_int(47 downto 0)));

It has a forced sign, so hopefully the triangular distribution is around 0 (from -1 to 1 by simplifying the idea). Keep in mind that lch_int has 67 bits. Basically dither_rnd0_int is the first generated output from LFSR and dither_rnd1_int is a previous dither_rnd0_int which is replaced with a new value each tick (the same clock is used for both). I have tried different ways of summing random values (signed, unsigned, addition, subtraction) to match the generic requirement it all does result in the same issue.

As soon as I add those pseudorandom values generated by the LFSR to low leveled signal like -96 dBFS several harmonic distortion (mostly 2nd, 3rd, 4th and 5th) do show up. This is not happening with signal -60 dBFS and above. In fact, there are no changes at all within higher signals (no distortion changes and no more noise).

Following are two measurements of a -96 dBFS signal with and without dithering to visualize the issue:

Without dithering:

With dithering:

The most interesting part is that a signal of level -70 dBFS does in fact gain a bit from using the same dithering:

Without dithering:

With dithering:

Somewhere around -75 dBFS is where dithering changes nothing and the lower we go it gets worse with dithering since it does bring strange harmonic distortions.

It's an old R-2R 18-bit DAC from late 80 / early 90, so don't mind the generic THD, just the difference.

The question is where do those extra harmonic distortions come from? Are they hiding within LSBs which I'm truncating and I do not see them unless I enable dithering which simply reveals them? If so, why would there be there? How can FIR coefficients introduce harmonic distortions? Considering the fact that I do quantize them on 35 bits with rounding I would not expect any issues above -150 dBFS or so.

This issue happens only if dithering is done AFTER filtering. If I dither input signal (before filtering and upsampling) the results are correct and quantization error is transformed into random noise. I did try to implement different coefficients with different windowing methods and types (linear or minimum), but the result is always the same.
This issue happens only when I want to reduce bit depth of filter output (67 bits) to 18 bits and only at low leveled signals. If I add dither to an input word of 32 bits effectively reducing it to 18 bits I do get correct results.

It seems like this is an expected behavior because -96 dBFS signals are too much correlated with the quanitzation error due to 18 bit trancation.

Adding more dither, like on 2 LSB does in fact get rid of those extra harmonics. However, the noise floor is getting much larger than expected. I don't think it's possible to have the best SNR in audible range while shaping the noise away to the end of the spectrum. I've tried that with 2nd noise shaper (TDPF dither with high-pass function) and it did not give any good results.


How to relocate face points in opencv / face distortion? I have to simulate facial expressions on a face image ( say open mouth ). For this I first extracted the facial feature points and found the corners of the lips. Now I need to deform the image by moving the points down.

In the above image I need to move the points ( 2 ) and ( 3 ) to some distance left and right respectively. And point ( 18 ) and ( 19 ) littele bit down. So that I will get an expression like opened mouth.

My Questions:

1) Is this the above way right to proceed to simulate facial expression?

2) If it is right how to move the points and create a new image in opencv?

Sure, that is one way of accomplishing your task so long as you do blending/interpolation to smooth the pixels in between. See these slides (esp. the last slide) or this survey for some literature on image warping and registration.

Thin-plate spline or "surface spline" problem
One reasonable approach would be establish point correspondences $r_1, \ldots, r_N$ and  $w_1, \ldots, w_N$ as you have done already then solve the following problem:

Find the 2D transformation model $\mathbf{p}=\{c_1,...,c_N,A,t\}$ s.t.
$$w_i=f(r_i)=Ar_i+t+\sum_{k=1}^{N}c_k\phi(d_i),$$
where $\phi(d_i)$ are radially symmetric basis functions (RBF). An example RBF choice is $\phi(d_i) = d_i^2\ln[d_i^2]$, where $d_i$ is a distance of your choice, e.g. Euclidean $d_i=\|r_i-r_k\|_2$.

Under certain thin-plate spline conditions not detailed here, this problem can be cast as a linear system $\mathbf{L}\mathbf{p}=\mathbf{w}$:
$$\begin{bmatrix}\mathbf{K} & \mathbf{P}\\\mathbf{P}^T & \mathbf{0}\end{bmatrix}\begin{bmatrix}\mathbf{c}_{N\times2}\\\mathbf{A}_{3\times2}\end{bmatrix}=\begin{bmatrix}\mathbf{r}_{N\times2}\\\mathbf{0}_{3\times2}\end{bmatrix},$$
where $\mathbf{K}_{ij\,\,[N\times N]}=\phi\left(\left\|(r_i^x,r_j^y)-(r_j^x,r_j^y)\right\|_2\right)$ and $\mathbf{P}_{N\times3}=\begin{bmatrix}\mathbf{1}_{N\times1}\left|\mathbf{w}_{N\times2}\right.\end{bmatrix}$, and then solve for $\mathbf{p}$.

You can then apply the transformation model $\mathbf{p}$ to every pixel in the original/control image, and use 2D interpolation to get the final result.

And to answer your second question, assuming you are using OpenCV's CvPoint structure, you can just do the following:

cv::Point newPoint( oldPoint.x + xoffset, oldPoint.y + yoffset )


Adjust the power of a digital signal to a given value I have two signals: $x(n)$ is a white-noise signal with a given variance, $y(n)$ is the sum of white noise plus a sum of sinusoids.
$$
x(n) = v_1(n)
$$
$$
y(n) = \sum_{i=1}^{q} [ a_i \cos(\omega_i n) + b_i \sin(\omega_i n)  ] + v_2(n)
$$
$$
v_{1,2}(n) \sim W.N.(0,\sigma^2)
$$

And that is how I generate the two signals on Matlab. Now I'd like to normalise $y(n)$ so that its power is the same as $x(n)$, how can I do that?


Ideal sampling of audio signals-Periodic spectrum When applying ideal sampling on audio singals, why is the spectrum periodic ?

Discritizing a signal in the ideal sense is equivalent to multiplying by a Dirac delta impulse train.  The Fourier transform of an impulse train is also an impulse train.  Lastly, a multiplication in the time domain is equivalent to a convolution in the frequency domain.  So when we discritize a signal, we are convolving that signal’s frequency transform with an impulse train, resulting in a frequency spectrum that repeats at intervals equal to the sample rate, making it periodic.


Frame acquisition vs symbol timing synchronization For a known length of frames  and using a known sequence in a preamble that is also  of known length can auto correlation at the receiver  be enough to acquire timing or is timing synchronization at the symbol level still needed?

Auto-correlation will be enough to determine the symbol center of 1 symbol at the end of the preamble, so that is sufficient for acquisition.

Symbol center tracking throughout the frame is still required to keep the receiver symbol clock synchronized to the transmitter's symbol clock.  Your receiver's symbol clock frequency will never be exactly the same as the transmitter's symbol clock frequency in real implementations, therefore continual estimation and correction or tracking is needed.


Is is necessary to learn digital signal processing before starting research in digital image restoration? I am going to start my research on digital image restoration, I just want to know that should i learn the course digital signal processing before starting my research.

My advice would be: do not delay the start of your research, but integrate the study of signal processing basics in the time you have for your research. Basic signal processing knowledge is essential for a proper understanding of the methods of digital image restoration. Divide your time between learning any basic knowledge you might be missing, and studying specific literature on image restoration.