Finding the system output by convolution pretty new to this concept, so do bear with me.

A linear dynamic system is described by the following differential equation:

Transfer function H(s) is calculated to be  =

I've already found the transfer function by assuming that the initial conditions are set to 0, getting the calculated stuff as shown above. I will then inverse Laplace it in order to find h(t), which turns out to be (exp^(-4t)t)

I will then use the formula, which is :

This is what I got, and I was wondering if it is correct:

If so, why did the questions give me the initial conditions? What should I do with them?
Also, how do I determine whether the system is stable? Unclear of this.

Thanks :)

As per EDIT:

I've converted them into differential equations where by:

In this case, how do I find the TF? I'm supposed to divide by X(s) in both sides, but the equation on the right will have a problem won't it?

EDIT 2:

You can't compute the output by convolving the input with the system's impulse response if the initial conditions are non-zero. The reason is simple: there is no information about the initial conditions in the impulse response or, equivalently, in the system's transfer function.

What you can do is use the following properties of the unilateral Laplace transform:

$$\mathcal{L}\{y'(t)\}=sY(s)-y(0^+)\tag{1}$$

$$\mathcal{L}\{y''(t)\}=s^2Y(s)-sy(0^+)-y'(0^+)\tag{2}$$

If you use $(1)$ and $(2)$ when you transform the differential equation, you'll get an equation for the Laplace transform of $y(t)$ which takes into account the given initial conditions. The inverse Laplace transform of that expression gives you the desired system output for non-zero initial conditions.

As for stability, causal systems with a rational transfer function are stable if all poles (i.e., the zeros of the transfer function's denominator) are inside the left half-plane.


SVD Based Noise Filter on the Samples in Toeplitz Form Given a sequence of length $2N$
$$
x[1], x[2], \dots, x[N], \dots x[2N]
$$
We can construct the following matrix along with its SVD.
$$
M =
\begin{pmatrix}
x[1] & x[2] & \cdots & x[N] \\
x[2] & x[3] & \cdots & x[N+1] \\
\vdots  & \vdots  & \ddots & \vdots  \\
x[N+1] & x[N+2] & \cdots & x[2N]
\end{pmatrix} = U S V^T
$$
Small singular values can be removed and matrix recomputed. Input signal can be replaced with e.g. 1st and last cols of $M$.
Unfortunately, I do not know the name and motivation behind this procedure.
Where can I find a reference with a detailed description of such a filter?
What is matrix $M$ called?
Example:

Test code in Mathematica:
(* define test signal and noise *)
length = 2^10 ;
amplitude = 1.0 ;
signal = amplitude*Sin[2*Pi*0.2741*Range[2*length]] ;
SeedRandom[1] ;
level = 0.1 ;
noise = RandomVariate[NormalDistribution[0.0,level*amplitude],2*length] ;
input = signal + noise ;
(* filter *)
matrix = Most[Partition[input, length, 1]];
{u,s,v} =SingularValueDecomposition[matrix,Tolerance -> 10^-15 ] ;
keep = 2 ;
s = Take[s,keep,keep] ;
u = Take[u,All,keep] ;
v= Take[v,All, keep] ;
matrix = u.s.Transpose[v] ;
result = {First[matrix],Rest[Last[matrix]]} //Flatten;
(* compare *)
Grid[{
{
ListPlot[signal,ImageSize ->300],
ListPlot[input,ImageSize ->300],
ListPlot[result,ImageSize ->300]
},
{
Periodogram[signal,Automatic,Automatic,KaiserWindow[#,Pi/2]&,PlotRange -> All, ImageSize ->300],
Periodogram[input,Automatic,Automatic,KaiserWindow[#,Pi/2]&,PlotRange -> All, ImageSize ->300],
Periodogram[result,Automatic,Automatic,KaiserWindow[#,Pi/2]&,PlotRange -> All, ImageSize ->300]
}
}]

The matrix M is called a Toeplitz matrix and they have a certain type of decomposition (but X should be the ACF and not the signal itself), so the resulting U and V would be a Vandermonde matrix which are very similar to Discrete Fourier transform.
So the whole process is similar to find spectrum of your signal and putting zero in bins which their power are below a threshold.


What are high frequencies and low frequencies in a signal? I'm new to signal processing. I'm plotting a signal in the time-domain in matlab. I don't understand which parts of the signal are high frequencies and low frequencies, can someone explain what a high and low frequency is, and how to see high frequencies and low frequencies on a graph?

The high frequencies contribute to the fast varying parts of the signal(the sharp transitions),
while the low frequencies contribute to the slow variations of the signal in the time domain.

You might want to take a look over here:
http://cns-alumni.bu.edu/~slehar/fourier/fourier.html

Also, if you have the time you might want to take a look at this online course:
http://academicearth.org/courses/the-fourier-transform-and-its-applications
You won't regret it.


Sampling rate conversion for audio using ffmpegwrapper in python I have many audio files @ 22kHz. I want to convert it to 44.1kHz. Since I have many files to convert, I am looking for an easy solution. I came across ffmpegwrapper for python. I installed the package for python. But I couldnt find out the command to resample. Can anyone please help me ? This is the link for the wrapper. https://pypi.python.org/pypi/ffmpegwrapper

Use SoX to easily convert your files in one go.

Here are some examples:

https://stackoverflow.com/questions/27264156/sox-batch-process-under-debian
https://pwnetics.wordpress.com/2011/07/12/bash-one-liner-using-sox-to-batch-convert-the-sampling-frequency-of-audio-files/

Alternatively, if you want to do this from Python, use librosa. An example:

By default, librosa will resample the signal to 22050Hz.

You can change this behavior by saying:

librosa.load(audio_path, sr=44100)
to resample at 44.1KHz, or

librosa.load(audio_path, sr=None)
to disable resampling.


Subtracting sensors output that have mutual inductnace between them Background: The software that i use to anaylze my signals is matlab. I have two audio signals that were recorded using two magnetic sensors. Lets call one sensor A and the other one B. A and B have mutual inductance between them.

While sensors A and B are working, I would like to subtract the information that was received in sensor A due to mutual inductance from sensor B.

I tried to subtract signal B from signal A, by simply writing in matlab (A-B), but it gives me a wierd answer.
I think that it derives from a shift that I have in phase. The beginning recording time of the two tracks is the same, and therefore I think that it is not a time delay.

I would like to know how to do this subtraction process in theory, and if anyone have any tips how to implement it easily in matlab I would like to know.

I will really appreciate any help.

Thanks in advance.

Below I attached pictures of the signals graphs.
In the first pic you can see the signals that are received from sensor A and sensor B.
In the second picture, I plotted sensor A in red and sensor B in blue, in the range of 2.12:2.16 and zoomed in even more.

This sounds like blind source separation.  In general, you can't unmix things after they've been mixed. If you have two different recordings of two sources with some of each source in each recording, you can sometimes use independent component analysis to separate them.

I have a Python example here.  There's also FastICA for MATLAB.  If they are audio signals, being picked up magnetically, there's probably no appreciable delay between them.  ICA works well in this case.


Aliasing correspondance with respect to DFT Given a sampling frequency Fs lets say we plot the magnitude of the fft of a temporal signal $x$, for different frequencies above the Nyquist frequency, to show the effect of aliasing:

Fs=10;
t=0:1/Fs:1-(1/Fs);
f=6;%7,8,9, ... to see aliasing starting a 6 since 6 is just above the Nyquist freq
x=sin(2*pi*f*t);
freq_range=Fs*((0:length(x)-1)/length(x));
plot(freq_range,abs(fft(x)));

Can someone help to formulate mathematically to which lower frequencies a given freqency f > Fs/2 will be aliased to?

Experimentally, we can see that (for integer frequencies, dont know if it really makes sense to use fractionnal ones?):

6 ->  4

7 ->  3

8 ->  2

9 ->  1

10 -> 0 <- should be zero? but the spectrum goes crazy and there are peaks everywhere.

Why? Then one would think that this would be periodic i.e. f=11 -> f_aliased=1, etc. back to zero in decreasing order, then again some strange things happens at 15 then it goes on cyclically etc. What is the relation?

Also, w.r.t. spectral resolution ($\Delta f = F_s/N$, let $N$ be the nb of samples, i.e. length(x)), one could ask if the next frequency abouve the Nyquist is in fact $f_{Nyq}+\Delta f = (F_s/2)+\Delta f$ ? But in this example the spectral resolution $\Delta f$ = 1 right?

If you do spectral analysis is typically better to use WAY more points then 10 in the time domain
Look at your y axis and make sure your graphs are properly scaled. $f=10$ aliases indeed back to $f=0$. All you see is numerical noise and your y-scale is $10^{-15}$. This would be different if you choose a a different phase, cosine instead of sine.


IDFT of upsampling in frequency domain Suppose I have a discrete time signal $x[n]$ of length 32 and I take it's 32 point DFT $X[k]$. Then I upsample this by a factor of 2 to get $Y[k]=X[k/2]$ and then take 64 point IDFT to get $y[n]$.

What will be the the relation between $y[n]$ and $x[n]$?

By upsampling in frequency domain, I know that $x[n/2]->X[2k]$. So by duality, $X[k/2]->Nx[-n2]$? What will $N$ be here, $32$ or $64$? I am confused on how to apply duality and whether it is required or not?

This answer does not require any elaborate mathematics. Rather it requires a careful interperation of periodicity of DFT sequences considering their lengths.

Lets display the signal chain below:

$$ x[n] \longrightarrow \boxed{ 32DFT} \longrightarrow  X[k] \longrightarrow \boxed{ \uparrow 2 } \longrightarrow Y[k] \longrightarrow \boxed{ 64-IDFT } \longrightarrow  y[n] $$

Now the signals are:

Time sequence $x[n]$ has a length of $32$ samples defined for $n=0,1,...,31$. For any other values of $n$, we interpret $x[n]$ as the periodic extension by the help of modulus operator as $x[n] = x[(n)_{32}]$. For example $x[35] = x[ (35)_{32} ] = x[3]$.

The frequency sequence $X[k]$ has a length of $32$ samples defined for $k=0,1,...,31$. For any other values of $k$, we interpret $X[k]$ as the $32$ point periodic extension of $X[k]$ similarly.

After the expansion block we have:

DFT frequency sequence $Y[k]$ has a length of $64$ samples defined for $k=0,1,...,63$. For any other values of $k$, we interpret $Y[k]$ as the periodic extension by $64$ samples similarly.

The IDFT signal $y[n]$ has a length of $64$ samples defined for $n=0,1,...,63$. For any other values of $n$, we interpret $y[n]$ as the periodic extension by $64$ samples similarly.

Now, simply write the $64$ point inverse DFT equation to define $y[n]$ :

$$
\begin{align}
y_{64}[n] &= \frac{1}{64} \sum_{k=0}^{63} Y[k] e^{-j \frac{2\pi}{64} k n } ~~~&,~~~ \text{ for } k,n = 0,1,...,63 \\
&= \frac{1}{64} \sum_{k=0}^{63} X[k/2] e^{-j \frac{2\pi}{64} k n } ~~~&,~~~ \text{ for } k,n = 0,1,...,63 \\
&= \frac{1}{64} \sum_{m=0, k=2m}^{31} X[m] e^{-j \frac{2\pi}{64} 2m n } ~&,~ \text{ for } m=0,1,..,31 ~~ n = 0,1,...,63 \\
&= \frac{1}{2} \left( \frac{1}{32} \sum_{m=0}^{31} X[m] e^{-j \frac{2\pi}{32} m n } \right) ~~&,~\text{ for } m=0,1,..,31 ~~~ n = 0,1,...,63 \\
y_{64}[n] &= \frac{1}{2} x[(n)_{32}] ~~~&,~~~ \text{ for } n = 0,1,...,63 \\
\end{align}
$$

The expression in the parenthesis above is simply the $32$ point inverse DFT of the sequence $x[n]$. Note that this sequence has a length of $32$ points and we can express it with the modulus notation $x[n] = x[(n)_{32}]$ as in the last line.

And the last line describes the $64$ point signal $y_{64}[n]$ as a two times repetition of $32$ point $x[n]$, divided by $2$. In other words, $y[0] = x[0], y[15]=x[15], y[40] = x[8], y[60]= x[28]$ etc...


Similarity theorem in Fourier analysis I have a homework problem that I'm not quite sure how to complete.  The problem is as follows:

PROBLEM

Write the definition of the Fourier coefficients, and show that

$$f(t + \frac{1}{2}T) = f(t) \rightarrow c_{2n+1} = 0$$
$$f(t + \frac{1}{2}T) = -f(t) \rightarrow c_{2n} = 0$$

Also show that this is a consequence of the similarity theorem in the Fourier analysis.

ANSWER

I have been able to solve the first part.  The definition of Fourier coefficients is:

$$c_n = \frac{1}{T} \int_{0}^{T} f(t) \exp(-2 \pi i n t/T) dt$$

Now let $f(t) = f(t + \frac{1}{2}T)$.  If we plug this in we get:

$$c_n = \frac{1}{T} \int_{0}^{T} f(t + \frac{1}{2}T) \exp(-2 \pi i n (t + \frac{1}{2}T)/T)dt$$

$$ = \frac{1}{T} \int_{0}^{T} f(t + \frac{1}{2}T) \exp(-2 \pi i n t/T) \exp(- \pi i n) dt$$

From this last expression we se that when $n$ is even, $\exp(- \pi i n)$ is $1$, and when $n$ is odd, $\exp(- \pi i n)$ is $-1$.  Thus we see that if we want $f(t + \frac{1}{2}T) = f(t)$ we must have all odd Fourier coefficients equal to $0$, and if we want $f(t + \frac{1}{2}T) = -f(t)$ then all even Fourier coefficients must be equal to $0$.

I don't know quite how to approach the last question though - how to show that this is a consequence of the similarity theorem.  Of course, the definition for  the similarity theorem is that:

$$f(t/ \gamma) \leftrightarrow |\gamma| F(\gamma \omega)$$

and a property is that if we "stretch" the function in the time domain, then we "shrink" it in the time domain, and vice versa.  But how this apply to the problem above?  My intuition tells me that since every second Fourier coefficient is set equal to $0$ we are "stretching" the function in the frequency domain, but I don't know how to show that we then get a "shrinking" in the time domain.

If anyone can give me some help in how to attack this last part of the problem, then I would be extremely grateful!

If you have a $T$-periodic function $g(t)$ with Fourier coefficients $c_n$

$$g(t)=\sum_{n=-\infty}^{\infty}c_ne^{in\omega_0t},\quad\omega_0=\frac{2\pi}{T}$$

and you create a new function $f(t)$ with Fourier coefficients $d_n$ defined by

$$d_n=c_{n/2},\quad n\textrm{ even}\\
d_n=0,\quad n\textrm{ odd}$$

i.e. the spectrum is expanded by a factor two by inserting zeros between the original Fourier coefficients. Then $f(t)$ is given by

$$f(t)=\sum_{n=-\infty}^{\infty}d_ne^{in\omega_0t}=
\sum_{n\textrm{ even}}c_{n/2}e^{in\omega_0t}=
\sum_{n=-\infty}^{\infty}c_ne^{i2n\omega_0t}=g(2t)$$

So you get a compressed time domain signal with period $T/2$, as directly shown in the first equation of the problem formulation in your question.


audio signal reconstruction from MFCC Transforming an audio signal to Mel Frequency Cepstral Coefficients is broadly used in tasks involving learning on audio. I was wondering, is this transform invertible with some good approximability properties? Are there methods to reconstruct an audio signal from its Mel Frequency Cepstral Coefficients?

Well, if the goal is to learn that some recording was similar to some reference, e.g. that the recording sound like the word "waffle", then yes, there's a way to reconstruct things, and that would be speech synthesis of "waffle", for example.

That's not a great example, but it's meant to illustrate one thing:

If your reduced version of the original signal, i.e. your feature vector (MFCCs in your case) is good enough to tell different things apart, and you can restrict the content of the signal to things that are actually well-represented by that feature set (e.g. you say "only voice, not music"), then yes, you can basically use closest-matches from a reference signal database.


Calculating information loss in a signal after smoothing? I have a signal. I applied Gaussian smoothing to it and then for baseline reduction I appied Tophat filter to the smoothed version.

This is the original signal:

This is the final signal:

I read that KL Divergence helps in finding the information loss between two signals but then again a condition was that the elements of the two distributions must be 1 i.e. There are two distributions P & Q then Pi + Qi = 1.

So is there any other way by which I can do that ?

Information loss and KL Divergence are functions that describe the behavior of random variables, in other words a set of random signals from some distribution. When you only have one signal, like you do in this case, calculating the information loss doesn't make sense. You can calculate mean squared error like the other comment suggests.


Why do we use keypoint descriptors? I just studied about SURF and I'm going for its implementation, but I still didn't understand why we use descriptors.

I understand what keypoints are and their purpose, but when we extract the keypoints than why do we need to use descriptors ? What is their importance and role in recognition?

One important thing to understand is that after extracting the keypoints, you only obtain information about their position, and sometimes their coverage area (usually approximated by a circle or ellipse) in the image. While the information about keypoint position might sometimes be useful, it does not say much about the keypoints themselves.

Depending on the algorithm used to extract keypoint (SIFT, Harris corners, MSER), you will know some general characteristics of the extracted keypoints (e.g. they are centered around blobs, edges, prominent corners...) but you will not know how different or similar one keypoint is to the other.

Here's two simple examples where only the position and keypoint area will not help us:

If you have an image A (of a bear on a white background), and another image B, exact copy of A but translated for a few pixels: the extracted keypoints will be the same (on the same part of that bear). Those two images should be recognized as same, or similar.

But, if the only information we have is their position, and that changed because of the translation, you can not compare the images.
If you have an image A (let's say, of a duck this time), and another image B, exactly the same duck as in A except twice the size: the extracted keypoints will be the same (same parts of the duck). Those are also same (similar) images.

But all their sizes (areas) will be different: all the keypoints from the image B will be twice the size of those from image A.

So, here come descriptors: they are the way to compare the keypoints. They summarize, in vector format (of constant length) some characteristics about the keypoints. For example, it could be their intensity in the direction of their most pronounced orientation. It's assigning a numerical description to the area of the image the keypoint refers to.

Some important things for descriptors are:

they should be independent of keypoint position

If the same keypoint is extracted at different positions (e.g. because of translation) the descriptor should be the same.
they should be robust against image transformations

Some examples are changes of contrast (e.g. image of the same place during a sunny and cloudy day) and changes of perspective (image of a building from center-right and center-left, we would still like to recognize it as a same building).

Of course, no descriptor is completely robust against all transformations (nor against any single one if it is strong, e.g. big change in perspective).

Different descriptors are designed to be robust against different transformations which is sometimes opposed to the speed it takes to calculate them.
they should be scale independent

The descriptors should take scale in to account. If the "prominent" part of the one keypoint is a vertical line of 10px (inside a circular area with radius of 8px), and the prominent part of another a vertical line of 5px (inside a circular area with radius of 4px) -- these keypoints should be assigned similar descriptors.

Now, that you calculated descriptors for all the keypoinst, you have a way to compare those keypoints. For a simple example of image matching (when you know the images are of the same object, and would like to identify the parts in different images that depict the same part of the scene, or would like to identify the perspective change between two images), you would compare every keypoint descriptor of one image to every keypoint descriptor of the other image. As the descriptors are vectors of numbers, you can compare them with something as simple as Euclidian distance. There are some more complex distances that can be used as a similarity measure, of course. But, in the end, you would say that the keypoints whose descriptors have the smallest distance between them are matches, e.g. same "places" or "parts of objects" in different images.

For more complex use of keypoints/descriptors, you should take a look at this question -- especially the "low-level local approach" in my answer and "Bag-of-words" approach in @Maurits answer. Also, links provided in those answers are useful.


How to represent the following system function in a block diagram? I am aware that this a simple question, but I haven't been able to find a single reference that explains this situation. The system function is:

$$\frac{1}{1-(0.3)(1+j)D}$$

How can one draw a block diagram with only real components to represent something that is complex?

The system function is

$$\frac{Y}{x}  = \frac{1}{1-(0.3)(1+j)D} \implies Y(1-(0.3)(1+j)D)=X$$
where $Y$ is the output and $X$ is the input.

Therefore $$Y-DY(0.3)(1+j)=X \implies y[n]=0.3y[n-1](1+j)+x[n]$$

It is a simple first-order recursive system like this:

In your case $b[0]=1$, $b[1]=0$, and $a[1]=-0.3(1+j)$. Of course the coefficient is complex if your system is complex. What you can do is split the real and imaginary parts and get two real-valued systems. The output of one is the real part of the output signal $y[n]$, and the output of the other is the imaginary part:

$$y_R[n]=0.3y_R[n-1]-0.3y_I[n-1]+x_R[n]\\
y_I[n]=0.3y_R[n-1]+0.3y_I[n-1]+x_I[n]$$

where $y_R[n]$ and $y_I[n]$ are the real and imaginary parts of $y[n]$, respectively (and the same for $x[n]$). If the input signal is real-valued then of course $x_I[n]=0$.

As for general complex-valued systems also have a look at this figure from the textbook Digital Communication by Lee and Messerschmitt:


Number of DFT (FFT) Points Required for a Specific Frequency Resolution for an Oversampled Signal I have a bandpass signal centered at 2 MHz and bandwidth of 50 kHz (the signal frequency varies from 2 MHz - 25 kHz to 2 MHz + 25 kHz). This signal is being sampled at 10 MHz.
I want a frequency resolution of 100 Hz in FFT.

I want to know:

How many samples of the signal should I take to achieve this frequency resolution.
What should be the number of points in FFT for this frequency resolution.

My thoughts:

Minimum frequency in the bandpass signal = 1.975 MHz.

To complete one period of the minimum frequency = 5.0633e-07 seconds.

No of samples in 5.0633e-07 seconds = 5.0633 ~ 6 samples.

So at least 6 samples should be taken to complete one cycle of min frequency.

Now the frequency resolution is 100 Hz.

Since the sampling frequency is 10 MHz, Maximum frequency can be detected is 5 MHz.

So 5MHz/100Hz = 50000 points will be there in first half of FFT.

Second half of the FFT is redundant (complex conjugation of the first half).

So I should take 2*50000 = 1,00,000 point FFT (2^17 = 131072 point FFT) for the above specifications.

Will this work?

The resolution in the DFT is given by: $ \frac{{F}_{s}}{N} $.
Hence you need 10e6 / 100 = 100,000 samples to get the resolution you want.
You may bring the signal to the baseband (demodulation) and then you'll need a lower frequency of sampling to achieve what you want.
Since the effective bandwidth of the signal is (Two Sided) 50 [KHz], it can be sampled at $ {F}_{S} = 50 [KHz] $ and only 500 samples would be required for the asked resolution.


Weird note structure on guitar audio waveform Why do some times notes take these weird structures (circled below). That is a guitar monophonic audio waveform by the way (3 note melody).

There are 3 notes played in this audio file. These are the beginning of the waveform (0), around sample 52000 and around sample 95000.

My two questions are:

Why did this weird looking negative amplitude peak get created (circled)?
Right after the final note (where the arrow is), why is there such a big change in amplitude (down, up, down again)?

The audio file is here.

Here are some possibilities:

Transverse waves on a guitar string have two polarized components, one parallel and one perpendicular to the guitar's soundboard. User hotpaw2 mentioned energy exchange. It is imaginable that energy transfer could happen between the components. The component parallel to the soundboard is largely inaudible, so it could stealthily store energy until released later.

Or it could be that you have two strings with modes that are slightly out of tune, and they form a beat frequency.

EDIT: Having analyzed the audio sample, the things you have marked are just some low-frequency artifacts, I think from things like someone touching or shaking the mic or moving around and bumping the guitar body. Adding a 100 Hz high-pass filter mostly gets rid of those:

Figure 1. Sample after 100 Hz highpass filtering.


Exact average error probability of QPSK in Rayleigh fading channel In Rayleigh fading channel(SISO), I know the instantaneous SER of QPSK for a given channel realization.
However, I do not know the exact average error probability.
It is very hard to calculate the average error probability because of the square of q function.
Are there any results on the exact average error probability of QPSK in Rayleigh fading channel?

From Tse and Viswanath, "Fundamentals of Wireless Communication" (which I highly recommend): $$p_e = \frac{1}{2}\left(1-\sqrt{\frac{\text{SNR}}{2+\text{SNR}}}\right)$$ Note that this is the probability of symbol error.


Minimum size for positive samples for HOG (Histogram of Gradient) detector training I'm going to train a HOG detector for detect objects (pedestrian head) with SVMLight library. But my question is not about implementation.

I'm wondering how to collect samples. I have to use some videos, and they have 1624x1234px dimension. The pedestrian heads have more or less 20x20px bounding-boxes.

Is the 20x20 pixels too small for a correct HOG feature extraction? Or is suitable for training? The, it has to be power of 2 so.. 36 pixels?

And what about the detection? Can I scale down the videos? As far as I know the detection is done on the training sample and then maybe scaled down. So if I train the detector with 20x20px I can detect objects of 20x20 as minimum, then if I scale down the image as half of the size I can detect object 40x40, but not 10x10 (I should scale up but the noise will be a problem).

And then, what about the computation time?

Any hint will be appreciated..

You should use elements which are a power of 2, for easy and faster computation.
I'll say, start with 8x8, so you can only scale up.
And if you have to look only at some portions of this image, then you can crop the original one.


Difference between DFT and Z-Transform I have searched this question but couldn't find the answer in this network. I know this is very confusing question for DSP beginners. Both DFT and Z-transform work for Discrete signal. I have read that "Z-transform is the general case of DFT, when we consider unit circle then, Z-transform becomes Discrete Fourier Transform (DFT)". What does this mean? Ok, I can understand the mathematical verification but what is the physical meaning of this and how this affect the analysis in DSP?

Actually, the Z transform is not really a proper transform, just a re-interpretation of the sequence of samples as coefficients of a formal Laurent series.

In some cases the formal Laurent series converges, if it does, it does so on an annular region in the complex plane. For useful signals (stable, summable, exponentially decaying) this annulus contains the unit circle, and the evaluation of the Laurent series on the unit circle corresponds to the Fourier series.

The interesting point of connecting a signal sequence to a periodic function on the unit circle is the inverse transformation, that many useful sequences are sequences of Fourier coefficients. And of course that convolution of signals corresponds to point-wise multiplication of the functions.


How to change a continuous signal to form a discrete signal? I have a continuous signal $$ x(t)= 2\cos(200\pi t)+8\cos(250\pi t)$$ and it is sampled at 400Hz.

How is this changed into the discrete time domain in order to determine the two discrete frequency components of $x[n]$ ?

When your continuous-time signal $$x(t)= 2\cos(200\pi t)+8\cos(250\pi t)$$ is sampled at a rate of $F_s = 400$ Hz, (or sampled with a period of $T_s= \frac{1}{F_s} = \frac{1}{400}$ seconds per sample), you can use the following to find the mathematical representation of the discrete-time signal $x[n]$ :
$$ x[n] = x_c(t)|_{t=nT_s} $$
Applying this to $x(t)$ leads to:
$$x[n] = x(nT_s)= 2\cos(200\pi n T_s)+8\cos(250\pi n T_s)$$
$$x[n] = x(n/400)= 2\cos(200 \pi n /400)+8\cos(250\pi n /400)$$
$$x[n] = 2 \cos(\frac {2\pi}{4} n) + 8 \cos(\frac {5\pi}{8} n)$$
Now observing the pattern of sinusoidal signals in discrete time as $\cos(\omega_0 n)$ where $\omega_0$ is the discrete time radian frequency between $-\pi$ and $\pi$ (or $0$ to $2\pi$ if you prefer so) you can deduce your discrete time frequencies as $\pi/2$ and $5\pi/8$
Note that $\omega_0 = \pi$ will be the highest frequency discrete time sinusoid.


Relation between impulse response and Region of convergence? Is there any relation between impulse response and Region of convergence?Means position/location of ROC depends upon impulse response or not??

Please kindly guide about cases for both continuous and discrete time?

For a discrete-time system if the impulse response is IIR and right sided, then the ROC will be outside of the largest pole. Else if the impulse response is IIR and left sided, then the ROC will be inside of the smallest pole.

For all FIR (finite inpulse response) systems, the ROC will be all $z$ except possibly zero and/or infinity.

For the continuous-time systems, for right sided impulse responses, the ROC will be to the right (half plane) of the max pole, and for left sided impulse responses it will be to the left of the minimum pole.


IIR biquad real-time filter just output noises I've implemented a IIR filter but it outputs only noise (like FM radio out of tuning) and I can't see what's wrong.

This is the processing function:

void Filter::process(int size, double *in, double *out)
{
    double xb[2] = { 0, 0 };
    double yb[2] = { 0, 0 };

    for (int n = 0; n < size; n++)
    {
        out[n] = (coeffs->a0 * in[n]) + (coeffs->a1 * xb[0])
               + (coeffs->a2 * xb[1]) - (coeffs->b1 * yb[0])
               - (coeffs->b2 * yb[1]);

        xb[1] = xb[0];
        xb[0] = in[n];
        yb[1] = yb[0];
        yb[0] = out[n];
    }
}

in and out are the buffers of sampled data at 48kHz; and size is 1024.

The coefficients are for $ {Fc \over Fs} = {500 \over 48000 } $ and $ Q = 0.7071 $

a0:  0.00107054
a1:  0.00214108
a2:  0.00107054
b1:  -1.99572
b2:  0.953753

I was unable to get a stable filter using your coefficients.  You did not state, but am assuming this is a lowpass filter.  Using the calculator here, I got:

coeffs = dict(
    a0=0.0010232172047183973,
    a1=0.0020464344094367946,
    a2=0.0010232172047183973,
    b1=-1.9075008174364765,
    b2=0.91159368625535,
)

Filter Code:

I recast to a Transposed Direct Form 2 (Not necessary, but often recommended for floating point) as:

def filter_process(input_data):
    zb = [0., 0.]

    output_data = []
    for data in input_data:

        output_data.append(coeffs['a0'] * data + zb[0])
        zb[0] = zb[1] + coeffs['a1'] * data - coeffs['b1'] * output_data[-1]
        zb[1] = coeffs['a2'] * data - coeffs['b2'] * output_data[-1]

    return output_data

Test Code:

import matplotlib.pyplot as plt
import numpy as np

def plot_f(f):
    Fs = 44000
    sample = 1000
    x = np.arange(sample)
    y = np.sin(2 * np.pi * f * x / Fs)
    out_y = np.array(filter_process(y))
    plt.plot(x, out_y)

for freq in (200, 500, 1000):
    plot_f(freq)

plt.show()

Results:

Update from Comments:

If the input and output buffers are in fact the same buffer, that definitely will be problematic. The transfer function as currently organized, writes a new output to the buffer prior to finishing it use as an input.  But that is easily remedied:

void Filter::process(int size, double *in, double *out)
{
    double xb[2] = { 0, 0 };
    double yb[2] = { 0, 0 };

    for (int n = 0; n < size; n++)
    {
        double in_n = in[n]
        out[n] = (coeffs->a0 * in_n) + (coeffs->a1 * xb[0])
               + (coeffs->a2 * xb[1]) - (coeffs->b1 * yb[0])
               - (coeffs->b2 * yb[1]);

        xb[1] = xb[0];
        xb[0] = in_n;
        yb[1] = yb[0];
        yb[0] = out[n];
    }
}


Point cloud comparison Let me preface this question with the fact that two hours ago I had no idea what a cloud point was. I do however have a solid mathematical background although geometry has never been my focus. I do know a lot about search trees, distances and space partitioning which I feel like would be useful for the potential solution to a problem I'm facing.

I have a 3D scanner that with high accuracy can generate a point cloud. Let's say I have a database of a number of objects in the range of 10,000. The way these are stored depends on the solution, but they could be point clouds themselves scanned by the 3D scanner, actual 3D models or perhaps even another representation of the object.

What I'm looking for is to match the newly scanned object to the 'closest' or most resembling object from the database and preferably tell me if there is no match as well (although this could just be some threshold I have to empirically set). I don't know if this is feasible in theory, feasible computationally or even reliable but I'm very interested if this is, I don't really know what to look for or where to start.

This is termed 3D object retrieval. There are many possible approaches to follow. First, if your meshes are nice, clean and descriptive, you could then consult the academic literature of computer graphics, while computer vision deals with more noisy scenes, partial visibility, non-meshed data. As I am a computer vision oriented person, I would be biased towards such approaches. But keep in mind that Eurographics organizes a workshop on 3d object retrieval every year, called 3DOR. Each year different problems are addressed, but it is very common to find something meeting your requirements.
The following might be good references in the direction you would like to go:

From Low-Cost Depth Sensors to CAD: Cross-Domain 3D Shape Retrieval via Regression Tree Fields

Database-Assisted Object Retrieval for Real-Time 3D Reconstruction

Feature-based Similarity Search in 3D Object Databases

There are many many more works though. For example, here are some, incorporating deep learning:

3DNet: Large-Scale Object Class Recognition from CAD Models

3D Shape Retrieval Using a Single Depth Image from Low-cost Sensors

3D ShapeNets: A Deep Representation for Volumetric Shapes


How to account for sample rate in a FIR Digital Filter? As a learning exercise I am trying to simulate a simple notch filter with Matlab. For test purposes  I use Matlab to simulate sinewaves of various frequencies and plot frequency vs. gain so I can see the effect of the “notch”.

I’m using the difference equation $$y(n) = x(n) - 2  \cos\theta\cdot  x(n-1) + x(n-2)$$ for this simple notch filter where theta is the notch (say $\pi/3$)

I thought x(n) would be the sample input from the sinewave I generate as test input. I'm just not sure how to account for the sample rate when determining the value of y(n).

From what I expect, and with comparison to using freqz([1 -2*cos(theta) 1],1),  I get a good match with a sample rate of 10 and a difference equation as shown below.

$$y(n)=  f_s [x(n) -2\cos(\theta/f_s)\cdot x(n-1) + x(n-2)]$$

but, as I change the sample rate there is a huge effect on the gain.

What am I missing?  My guess is I just don't understand how to account for a change to the sample rate in the difference equation but I haven't seen a good write up on this anywhere

Can anyone help?

The simple code is shown below.  As I change sample rate (from 10 to 100 and 1000) the frequency response appears to be correct but the gain appears to be decreasing by 20dB for each power of 10 I increase the sample frequency. Maybe that’s the clue?

..... CODE ..  (You can tell I'm quite new to Matlab also... apologies for style!)

fs=10;  %Specify Sampling Frequency

Ts=1/fs; %Sampling period.

fb=100;    % Number of frequencies to try

mb=zeros(fb,2);   %Results   (Max / Freq)

w=0.001;           %Frequency of input wave to be samples  (Don't use DC!)
fb_indx=1;         %Index into the results array

THETA=pi/3;          %For a simple test
while fb_indx < (fb+1)

    Ns=round(2*pi*fs/w);  %Nr of time samples to be plotted.

    t=[0:Ts:Ts*(Ns-1)];  %Make time array that contains Ns elements
                     %t = [0, Ts, 2Ts, 3Ts,..., (Ns-1)Ts]

     mb(fb_indx,1)=w;  %Store current test frequency

     y=zeros(1,Ns);     %Result (y(n))

     x=sin(w*t); %create sampled sinusoids at different frequencies

% Here for the notch

  for n = 3:Ns  % loop for number of samples

     y(n) = fs*(x(n)-(2*cos(THETA/fs)*x(n-1))+x(n-2));  % calculate y

  end

% Obtain max. response
  M=max(y);                %Look for highest peak in result

  mb(fb_indx,2)=20*log10(M);  %and store in 20log10 value

    fb_indx=fb_indx+1;

      w=w+(pi/fb);  %Get next frequency

  end;

% Plot result

figure % create new figure

subplot(2,2,1)           % first subplot

plot(mb(:,1),mb(:,2));
xlabel('Freq - rad/s');
ylabel('20log10 Magnitude');
axis([0 pi -100 20]);

grid on;

subplot(2,2,2);
plot(x);
title('Input');

subplot(2,2,3);
plot(y);
title('Output');

The sampling frequency is hidden in your variable $\theta$. The angle $\theta$ determines the frequency of the notch relative to the sampling frequency, i.e.

$$\theta=2\pi\frac{f}{f_s}\tag{1}$$

where $f$ is the notch frequency in Hertz, and $f_s$ is the sampling frequency (in Hertz). So if you want to change the sampling frequency while leaving the notch frequency unchanged, you need to change the angle $\theta$ according to Eq. (1).

The gain of your filter depends on the normalized frequency variable $\theta$. At DC the gain is $2(1-\cos\theta)$, and at Nyquist you get a gain of $2(1+\cos\theta)$. By normalization you can keep one of those two gain factors normalized to $1$, but never both. If you want a gain of $1$ at DC and at Nyquist you need an IIR filter or a higher order FIR filter.


Interpreting magnitude of DFT results I'm working on creating a simple program to render spectrograms like this one. In this plot, the X-axis is time, the Y-axis is frequency, and the color represents the magnitude of the DFT at that frequency and time.

I have most of this working, but I'm not exactly sure how to interpret the magnitude value. I'm computing it as follows:

real * real + imaginary * imaginary

(Note that I'm not square rooting it, since Wikipedia indicates that the spectrogram should be related to the square of the magnitude.)

I've seen certain spectrograms which work with magnitude on the dBFS scale. What's the process to translate these magnitudes to this scale? Or, alternatively, is there some easier way to go about this?

One common practice is to take the log of the squared magnitude of the DFT frames, as in log(re*re+im*im), scale that, and use the scaled log as an integer index into an RGB color lookup table ("hotter" colors for larger magnitudes).  The scale is determined so that the range in your data fits within the color index table.

A log scale is used, since it's a closer match to human perception of loudness than a linear scale.

A square root is un-needed, as that's the same as adjusting the scale after the log by 0.5.


Can I study continuous time Fourier Transform and treat the rest as special cases Say I learned the theoretical result of continuous time Fourier transform. And I want to extends some results(say "convolution rule") to

Lapace transform,
Z transform,
DTFT,
DFT,
Fourier sequence etc.

Then really, I just need to remember how one works and derive the rest. How to do this as a general principle? Is there any textbook that teaches this way? Thanks.

Yes you can, somehow. But that won't be the easiest road.

When teaching basic DSP, I often see students afraid to "learn all the formulae" and properties pertaining to the Fourier transform, Fourier series, discrete-time Fourier transform, $z$-transform and the Discrete Fourier transform with the fast
 Fourier transformation or FFT, up to time-frequency transformations.

They share a lot of common traits: energy conservation, sorts of orthogonality, scalar product formulae, convolution turning into multiplication, Parseval-Plancherel formulae, etc. but also important (and often subtle) differences.
.
One tool to unite such concepts is sometimes called the Pontryagin duality. For a basic overview, the wiki page is ok:

  Pontryagin duality places in a unified context a number of
  observations about functions on the real line or on finite abelian
  groups

For instance, you can view the four primal and dual domains for Fourier "transformations" as in the table:

\begin{array}{lcc}
\textrm{Transform} & \textrm{Original domain}& \textrm{Transform domain}\\
\textrm{Fourier transform} & \mathbb {R}  &\mathbb {R} \\
\textrm{Fourier series}& \mathbb {T} & \mathbb {Z} \\
\textrm{Discrete-time Fourier transform (DTFT)} & \mathbb {Z}  & \mathbb {T} \\
\textrm{Discrete Fourier transform (DFT)}& {\displaystyle \mathbb {Z} /(n)} & {\displaystyle \mathbb {Z} /(n)}
\end{array}

Its provides several ways to derive properties from one domain to another.
But this tool dives deep into complicated concepts such as category theory, see for instance Notes 2: The Fourier transform (Terry Tao) or Category theory applied to Pontryagin duality.

Knowing properties about one transformation, you can derive properties of others, but for that you will have to master some advanced mathematics. There is no free lunch, unfortunately.

You would better learn how to remember a few formulae for a handful of transformations


proof of Autocorrelation property of DFT I am facing problem in proving the Auto-correlation property of Discrete Fourier Transform (DFT), that is

$$\mathcal{DFT}\left\{\sum_{r=0}^{N-1}x[r]x^*[r+n]\right\} = X[k]X^*[k]= |X[k]|^2$$

where $X[k]$ is the DFT of $x[n]$, and $x[n]$ is periodic with period $N$:
$$ x[n+N] = x[n] \qquad \forall n \in \mathbb{Z} $$
Likewise  $X[k]$ is periodic with period $N$:
$$ X[k+N] = X[k] \qquad \forall k \in \mathbb{Z} $$

Below is my Approach:

$$\begin{align}
\mathcal{DFT}\left\{\sum_{r=0}^{N-1}x[r]x^*[r+n]\right\} &= \sum_{n=0}^{N-1} \sum_{r=0}^{N-1}x[r] x^*[r+n] e^{- j 2 \pi n k/N} \\ \\
 &= \sum_{r=0}^{N-1}x[r] \sum_{n=0}^{N-1}  x^*[r+n] e^{- j 2 \pi n k/N} \\ \\
 &= \sum_{r=0}^{N-1}x[r] \sum_{n=r}^{r+N-1}  x^*[n] e^{- j 2 \pi (n-r) k/N} \\ \\
 &= \sum_{r=0}^{N-1}x[r] \sum_{n=0}^{N-1}  x^*[n] e^{j 2 \pi r k/N} e^{- j 2 \pi n k/N} \\ \\
 &= \sum_{r=0}^{N-1}x[r]  e^{j 2 \pi r k/N}  \sum_{n=0}^{N-1}  x^*[n] e^{- j 2 \pi n k/N} \\ \\
 &= \sum_{r=0}^{N-1}x[r]  e^{j 2 \pi r k/N}  \sum_{n=0}^{N-1}  x^*[n] (e^{- j 2 \pi n (-k)/N})^* \\ \\
 &= \sum_{r=0}^{N-1}x[r]  e^{j 2 \pi r k/N}  \sum_{n=0}^{N-1} ( x[n] e^{- j 2 \pi n (-k)/N})^* \\ \\
 &= \sum_{r=0}^{N-1} x[r] e^{ j 2 \pi r k/N} X^*[-k] \\ \\
 &= \sum_{r=0}^{N-1} x[r] e^{ -j 2 \pi r (-k)/N} X^*[-k] \\ \\
 &= X[-k] X^*[-k] = \Big| X[-k] \Big|^2 \\ \\
 & \ne \ X[k] X^*[k] \quad = \Big| X[k] \Big|^2 \\
\end{align}$$

as $$X[k] = \mathcal{DFT} \{ x[r] \} = \sum_{r=0}^{N-1} x[r] e^{- j 2 \pi r k/N}$$

so where is my mistake? and how should I approach??

You did not consider the conjugate property of the DFT.

The complex conjugate property of DFT states that
$\rm{DFT}\{x^{*}[n]\} = X^{*}[-k]$.
This implies that exponent in the second line of your proof should be $-k$. Which then allows for the DFT of $x[n]$ as $X[k]$.
A number multiplied by its conjugate give the square of its absolute value.

And DFT is discrete in the transform domain, hence its parameter should be in brackets rather than parentheses.


Simulating multipath channel I'm trying to simulate a multipath channel to see how it would affect my signal. So considering my channel is LTI, I can model it as an FIR filter. I've been reading Telecommunication Breakdown and looking at some examples and there's one thing I don't quite understand in page 183 (Chapter 9: Stuff happens, Section 9.4: Other Impairments: "More What Ifs")

if cdi < 0.5, % channel ISI
    mc=[1 0 0 ] ; % distortion−free channel
elseif cdi<1.5 , % mild multipath channel
    mc=[1 zeros (1 ,M) 0.28 zeros (1, 2.3*M) 0.11 ] ;
el se % harsh multipath channel
    mc=[1 zeros (1, M) 0.28 zeros (1, 1.8*M) 0.44 ] ;
end

So cdi is a parameter to choose between distortion-free, mild or harsh multipath channel and mc is the resulting FIR filter. What I don't understand is:

What's the point on inserting zeros in between?
How does that make it mild or harsh?

Would you mind explaining that to me? Thanks

You put zeros in between to show delay between the channel taps. I don't know what $M$ is in your code but suppose $M=10$ and say your sample rate is 10 MHz. Then you can interpret mc = [1 zeros(1, M) 0.28 zeros(1, 2.3*M) 0.11];  as you get the first multipath component with zero delay (gain = 1), then 10 zeros later (which is equal to 1 microsecond) you get another multipath component (gain = 0.28), and finally after another 2.8 microseconds you receive the last multipath component (gain = 0.11).
It looks like the "mild" vs "harsh" is just language that the author uses to let you know that mild happens to rarely cause errors where the harsh causes many errors.


Obtain the complex envelope of a two-tone signal I know it is possible to represent a single-tone signal by its complex envelope:
\begin{align}
x(t) &= A(t)\cos(2\pi f_c t + \Phi(t))\\
x(t) &= \operatorname{Re}\left\{A(t)\mathrm{e}^{j \Phi(t)}\,\mathrm{e}^{j 2\pi f_c t}\right\}\\
x(t) &= \operatorname{Re}\left\{\tilde{x}(t)\,\mathrm{e}^{2\pi f_c t}\right\}
\end{align}

I think it should also be possible to represent a two-tone signal by its complex envelope:

$$x_2(t) = A_1(t)\cos\left(2\pi f_1 t + \Phi_1(t)\right) + A_2(t)\cos\left(2\pi f_2 t + \Phi_2(t)\right)$$

Can somebody show me the necessary steps to obtain the complex envelope of $x_2$ ?

You can represent any bandpass signal by its complex envelope.

You must, however, pick a frequency for which you do that. In your first example, that frequency is $f_c$; in your second example, you need to pick $f_1$, $f_2$ or any other frequency as $f_c$ – and then represent your cosines into products of $e^{j2\pi (f_1-f_c) t}$ and $e^{j2\pi f_c t}$.

Any textbook that introduces the concept of complex envelope (and, usually, complex baseband) should introduce that concept; if yours doesn't, it might be worth having a second textbook for cross-literature.


What's the similarities and differences between Wigner transform and wavelet transform? Wigner transform and continuous wavelet transform are both some kind of time-frequency representation of a signal.
What are the similarities and differences between them? Could you give some comparison between them? Let's restrict ourselves in 1D signal at the moment.

There is a nice tutorial: Time-Frequency Toolbox For Use with MATLAB tthat compares many transformations, notably continuous wavelet transform and Wigner and its avatars.
Basically, the CWT belongs to atomic linear representations (a decomposition onto wavelet atoms), while Wigner-Ville is an
energy distribution (to split the energy of signals over  time and frequency).
Since the tutorial comes with a toolbox and many examples, you can easily reproduce it on your signals.


What are some python packages I can use to cut audio files I have a project in which I have a batch of audio files and I need to remove the audio in it from say time 2secs to 5secs i.e. I have a set duration to be removed. How can I do this in python?

The librosa library is usually used. You can find the files here:

https://github.com/librosa/librosa


What does "frequency" mean for various kinds of signals? I've been involved in Signal processing for quiet a time but I'm still so confused of what frequency could possibly is, as it has different meanings in different scenarios, for example,

According to Wikipedia,

  Frequency is the number of occurrences of a repeating event per unit
  time.

  Calculating the frequency of a repeating event is accomplished by
  counting the number of times that event occurs within a specific time
  period, then dividing the count by the length of the time period. For
  example, if 71 events occur within 15 seconds the frequency is:

This is the most basic definition of the frequency that everybody knows. But what is the definition of Frequency in case of Digital Images and Sounds ?

  For example, if a newborn baby's heart beats at a frequency of 120
  times a minute, its period (the interval between beats) is half a
  second.

make sense so far.

Now here is the Sine wave of various frequencies,

the bottom waves have higher frequencies than those above. The horizontal axis represents time.

The above signal makes sense to me that it has a frequency but how about a non periodic signals like a human voice, ?

Have a look,

this signal is not repeated on any instance so how one could say what is its frequency and how one could count the number of cycles repeated ?

aforesaid,

  Frequency is the number of occurrences of a repeating event per unit
  time.

How does this statement can be true in case of the frequency of human voice ? because when we speak we do not repeat anything than how the voice can have some frequency ?

and in case of images,

  the rapid change in the color is the part of the High frequency of
  the image

How this can be counted as frequency ?
if the image has all the different pixel values how there could be any frequency ?

I am so curious to know deep in detail about the exact definition of Frequency that holds for everything mentioned above.

The key insight that Fourier had when he developed Fourier analysis is that any absolutely integrable (thanks Jason R) function can be represented as the weighted sum of sines and cosines.  Explaining why this is true is way beyond the scope of this answer.  I suggest you study Fourier theory to understand this better.


Finding Peaks in an Autocorrelation Function I'm trying to find the period of a signal. I've used FFT to compute the autocorrelation of the signal. As can be seen from the autocorrelation function (plotted below) I obtained, there are 70 sample between peaks which actually indicates the period of my signal.

What is the best way to extract the indices of these peaks from such a data?

Removing the DC offset from your signal will get rid of the triangular "trend" seen here. Another way to detrend data (which is not specific to autocorrelation functions) is to subtract from your function a median-filtered version of itself (the median-filtered version corresponding to the trend).

You can then detect peak by detecting local maxima - if $X(n) = \max_{k \in [k-W, k+W]} X(k)$, then $n$ is a peak. $W$ is a scale factor which indicates how narrow and close to each other you allow your peaks to be.


With regards to multirate processing, what are the benefits of having a 'slower' sampled signal? With regards to multirate processing, what are the benefits of having a 'slower' sampled signal?

Let's say I have a continuous signal

$$x(t) = \cos(2\pi\cdot3t) + \cos(2\pi\cdot 7t)$$

and it is sampled at $10\textrm{ kHz}$. Obviously the sample rate is overkill for the underlying signal. I've started reading into multirate processing, and it often states that the usefulness of resampling (decimation or interpolation) can come in the form of lower processing requirements.

What would be the benefits if we subsampled this $10\textrm{ kHz}$ sampled signal?
So my question:

How exactly do we save processing effort?
Is it simply due to the fact that we can subsample a discrete signal, obtain a new signal with lower time resolution (but still adequate frequency resolution and range), and then only have to perform operations on a smaller set of data?

Also, do the benefits of multirate processing differ depending on whether or not we're doing realtime or offline processing? It seems like in the case of real time processing, subsampling would allow you to use DSPs that operate at slower frequencies.

Please feel free to offer any insight to the topic!

In addition to what hotpaw2 said, it's important to understand that sample rate reduction is the corner stone of actual applicability in many systems.

Software Defined Radio frontends nowadays produce hundreds of complex-valued Megasamples per second – if you don't reduce that amount of data, no commonly available PC can deal with that.

Also, many devices that you meet nowadays actually have software-programmable DSP embedded very deeply. The average mobile phone baseband IC, i.e. what handles the signal, is digital and very broadband, to "get" the whole bandwidth that e.g. LTE demands they can capture; however, at a single time, only a fraction of that bandwidth might be relevant. So isolating/extracting that bandwidth is crucial to functionality.

Not reducing the sample rate in that situation would not only be a pure waste of power – it would render the following steps, such as equalization and synchronization, totally impractical in the telephone form factor.

A movie camera might take up to let's say 16 Million pixels – which are, in the end, nothing but samples of a twodimensional signal! If the producing camera didn't remove that sample rate (i.e. the picture's resolution), there's no way that the resulting movie data would fit through the terrestrial TV channel; no way that a cheap TV might decode it, and no way you could by such a movie on a cheap optical disc at current capacities.

So, this, by far, is not a "historical" problem. It is a very current problem – not being able to handle the full sampling rate that ADC/DAC hardware offers is the core technological problem of pushing communication further.

Now, 10 kHz really isn't something you'd normally worry about – that's something that most algorithms could deal with on the cheapest microcontroller. But think more in the terms of a SDR frontend (e.g. this one) that has 200 MHz of sampling rate per single complex ADC; that means that you can't even get the  two channels of that device into whatever computer you attach at the native sampling rate, because $2\cdot200 \cdot 10^6 \frac {\text S}{\text s} \cdot \frac{16 \text{ b [I]} +16 \text{ b [I]} }{\text S} = 12.8 \frac{\text{Gb}}{\text s} > 10 \frac{\text{Gb}}{\text s}  $, which is the datacenter-grade 10GbitEthernet interface's data rate that usually totally flattens out even the most powerful quad-core PC CPUs, without them even doing any DSP, just shuffling the data into RAM. Now, you realize that for almost all applications, it's very desirable that the device internally reduces the sampling rate of the observed signal as far as possible to correctly represent that signal, for computers to even be able to cope with the raw amount of data. The device does that, by applying cascaded multirate filters internally. The ability to, nevertheless, observe the whole bandwidth at once and "select" any part of these 200 MHz at random is very crucial to designing things like wideband transceivers or frequency-hopping standards; so the ability to use the minimum viable sampling rate is very central to progress in research, application development and implementation, as soon as you're not considering 10 kHz systems any more.


What's a "Fourier filter"? E.g. the constant Q-transform is built by adding so called "Fourier filters".

What's a "Fourier filter"?

People (usually from fields outside signal processing) sometimes use the term Fourier filter for a filtering operation in the FFT domain, which simply works by multiplying the FFT bins of a signal with a given filter function (often just ones and zeros, corresponding to pass bands and stop bands, respectively). Why this is generally not such a good idea is explained here.

Also in Computer Vision, the term Fourier filter is used as explained above.

In the document you linked to in a comment, the term is used to describe the computation of the Discrete-Time Fourier Transform (DTFT) at a given frequency from a finite length portion of a signal. This computation can be interpreted as a filtering operation, because it is a sum of products. The corresponding filter is a band pass filter with center frequency equal to the given DTFT frequency. More more information on the filter interpretation of the D(T)FT have a look at this page.


Magnitude response and DFT normalization Suppose I have an FIR denoted h that represents the impulse response of a system. Using MATLAB syntax (for convenience and brevity),

What does abs(fft(h)) represent?

If the answer is the system's magnitude response, then please proceed to the next question. If not, please correct me.

Is it correct to say that the value of the magnitude response of a system at a given frequency represents the gain applied to an input sine tone at that frequency?

If yes, then please proceed to the next question. If not, please correct me.

If I upsample h by, say, a factor of 2, then the values of abs(fft(h)) are doubled compared to the same computation done at the original sampling rate, at least up to the Nyquist frequency of the original rate. This seems to contradict the answer to question 2, because even though the impulse response has twice the number of samples, it represents the same system and therefore, if abs(fft(h)) represents the magnitude response (with magnitude response defined as mentioned above), then there shouldn't be a discrepancy.

These questions should reveal gaps in my knowledge regarding the meaning of a system's magnitude response, and normalization when computing a DFT (see the answers to this related question). I'd appreciate some clarity on these points. Thanks!

Update 1: Below is a figure showing an example.

The chosen example FIR is a low-pass filter. Blue is the original sampling rate (1024 samples/second). Red is upsampled by a factor of 2. The system magnitude response is a result of doing abs(fft(.)) to each of the IRs shown on the top left. The two plots below (bottom left should actually read "spectrum" and not "response") shows what the system does to an input, unit amplitude sinusoid of frequency 50 Hz.

Update 2: The confusion has been resolved in the comments section of the accepted answer. The source of confusion was the assumption that upsampling a discrete-time impulse response (as shown above) is exactly the same as sampling an analog IR at the higher sampling rate. This is not the case.

There is no contradiction here. The first two points in your question are correct. In the third one the mistake lies in the assumption that "it represents the same system". It doesn't. You've shown yourself that the response to a sinusoid is different for the two systems; the output of the interpolated system is larger than the other one, which corresponds to the larger value of the magnitude response.

As a simple test, note that the value of the frequency response at $\omega=0$ is given by the sum of all samples (assuming a signal of length $N$):

$$X(0)=\sum_{n=0}^{N-1}x[n]\tag{1}$$

Clearly, this sum is different for the original signal and for the interpolated signal. Note, however, that for upsampling and interpolation by a factor of two, the difference is not necessarily a factor of two. The exact factor depends on the signal and on the type of interpolation.


What is the Fourier Transform of a constant signal? I am trying to figure out what the fourier transform of a constant signal is and for some reason i am coming to the conclusion that the answer is 1. Or better yet a step function.

I'll complete a bit the answer given in a comment above.

Intuitively first, to which frequency corresponds a signal constant in time, for exemple $x(t) = 1$ $\forall t$ ? Such a signal shows no variation in time and hence contains only a component with frequency 0 (this is a DC signal). This means that its  Fourier transform must be 0 everywhere, except in $f=0$. Mathematically,
$$X(f) = \delta(f).$$
Now, can we prove this? Yes, simply take the inverse Fourier transform of $\delta(f)$ and use the properties of the Dirac delta $\delta(f)$
$$x(t) = \int_{-\infty}^\infty \delta(f)e^{j2\pi ft} \mathrm{d}f = \int_{-\infty}^\infty \delta(f) \mathrm{d}f = 1.$$


Causal and Non-causal filtering. Strange behavior I'm trying to understand the behavior of a causal filtering exemplified by figure below.

This behavior is making LDA (Linear Discriminant Analysis) misclassify our data, by clearly reverse the signal pattern at the highlighted box.
We need to alter the filtering method from non-causal (MATLAB filfilt function) to a causal (MATLAB filter function), but we face off this strange behavior.

Can someone give us an idea for what is really happening?

EDIT:

Motivation: Real-time processing.
-Note: Delays are expected, but the output causal filtered signal inside the red box seems to be very uncorrelated with the zero-phase filtered signal.
Band-pass cut-off frequency: $[0.16 - 1]\textrm{ Hz}$.
Signal: Electroencephalography @ $128\textrm{ Hz}$

the output causal filtered signal inside the red box seems to be very uncorrelated with the zero-phase filtered signal.

So let's look at this statement. Suppose the original, unfiltered signal is $x[n]$. Filtering this using filtfilt with a filter with impulse response $h[n]$ gives:
$$
y_{ff}[n] = h[-n] \star h[n] \star x[n]
$$
and filtering with just filter is:
$$
y_f[n] = h[n] \star x[n]
$$
so
$$
y_{ff}[n] = h[-n] \star y_f[n]
$$
Regardless of what $x[n]$ is, the 'correlation' between $y_f$ and $y_{ff}$ will look like the anticausal (perhaps anticausal + constant term) $h[-n]$.


Method to rescale signals to mean length I have a set of signals of varying lengths. I have provided an example of the same below -

Their lengths vary between 186 to 202, with a mean length of 197. I am looking to rescale them to the mean length. I am looking for ideas to do the same. Perhaps a good start will be reinterpolation, but I am open to other approaches too.

If you really want to scale data, you can look at the scale transform, detailed in The scale representation, L. Cohen, 1993. There also exists a time-scaling method given in A. Papoulis, "Signal Analysis", McGraw-Hill, p. 290, 1977; one Matlab/Scilab/Octave implementation is provided in the DiscreteTFDs toolbox by  Jeffrey C. O'Neill.For the latter, you just have to provide a scaling factor proportional to the inverse of the signal's length.

However, in the same line as @Florian, and looking at the signals, I would be more prone to 1) find the maxima at fractional steps through interpolation (eg Quadratic Interpolation of Spectral Peaks) 2) realign with fractional filters (Splitting the unit delay: tools for fractional filter design) 3) truncate signals.


Spectral centroid manipulations So, I've created a simple sound analysis application and one of the features I've implemented is the spectral centroid (as explained here https://en.wikipedia.org/wiki/Spectral_centroid). In order to get reliable results, we need to set a threshold on our spectral centroid data because the FFT always returns some energy in bins that really don't have any energy.

Now, I would like to be able to manipulate the spectral centroid (e.g push the values up or down, which to some extent resembles pitch shifting) and then use an inverse FFT operation to convert the values to waveform data. This way, we'd be able to hear the changes but I'm not sure how to do this.

The spectral centroid is calculated by using magnitudes, not raw FFT data. The magnitudes are calculated as given below

mag[i] = 10*log10(sqrt((pFFTLReal[i] * pFFTLReal[i] + pFFTLImg[i] * pFFTLImg[i])/fftNorm));//fftNorm is basically a correction factor that depends on a type of window (blackman, hamming etc)

So, given all this, how would I change the spectral centroid by, say, 10 percent and see that change in the generated (manipulated) waveform?

Is this even possible?

Please, if this is not clear enough, I'll try to provide some more information.

Thank you

It's been some time since I asked this question and after some work done on this subject, I think it's time to revisit it. It should be noted that
while the spectral centroid, pretty much like any other spectral feature, is calculated using the FFT magnitudes, and not raw FFT data, we don't
have to manipulate the magnitudes directly. Instead, we modify FFT data, which in turn, affects the magnitudes. This makes the process
a lot easier than having to manipulate the magnitudes.

The spectral centroid is a single value (per FFT frame) so in the case of STFT we have as many centroid as we have STFT frames. The process of manipulating
the spectral centroid is pretty straightforward. If, for instance, the spectral centroid is 600Hz and we wish to increase its value by 10% so that the new value is
660Hz, we would simply rearrange our frequency bins by a certain amount (a bin step). The actual calculation is qute simple

const binStep = frameLenHalf - frameLenHalf/(1+(float)freqDelta/100);
   //frequency delta % (e.g. 10%, 20%, -10%, -20% etc)

    //for increasing frequency
    float *pCfftLReal = new float[frameLen];
    float *pCfftLImg =  new float[frameLen];
    //for increasing frequency
    for(j=frameLenHalf-1;j>binStep;j--){
            pCfftLReal[j] = pFFTReal0[j-binStep];
            pCfftLImg[j] = pFFTImg0[j-binStep];
                }
        }
    //for decreasing frequency
    for(int j=0; j<frameLenHalf-1-binStep; j++){
            pCfftLReal[j] = pFFTReal0[j+binStep];
            pCfftLImg[j] = pFFTImg0[j+binStep];
                }

Needless to say, this process suffers from rounding errors (primarily calculating the binStep), but in general, the larger the bin resolution, the larger the errors.
This can be alleviated by having a greater frame length (but for STFT this is usually undesirable). Other more complex approaches are also possible but I haven't investigated any of these.

Also, sometimes the results can be skewed if the amount of frequency shift exceeds the Nyqiust frequecy and we keep the original sample rate (upsampling would help here).

Before we can actually listen to our manipulated sound, an inverse FFT shoud be carried out on our modified data.

What I've described here is just a rough, naive approach to spectral centroid manipulation. It should also be said that this kind of manipulation
doesn't retain the original harmonics ratio, given that the frequency delta shift is applied equally on every bin. Keeping the harmonics ratio intact would require
calculating a new bin step for every frequency, which is only slightly more complex that the algo given above.

I've tried both approaches (with a constant bin step and a changing bin step) and perceptually, the first doesn't really change the original sound ( I guess
this has to do with higher frequencies affecting timbre more than lower ones and the first approach changes the higher frequencies less than it does the lower ones).

I'm not done with this yet, so I will update this post as soon as I have something interesting to share.


What are the fundamental reasons that suggest us to minimise transmission power for energy-efficient radio link design? My question in short:

Under which scenarios minimising transmission power is beneficial in terms of transmitter energy-efficiency? is there any fundamental reason to keep a low-power oriented design strategy?

In particular, small, low-cost transceivers (such as those used for M2M/IoT communications) are designed according to this strategy. I'd like to understand why.

Background

It is quite natural to think that energy efficiency in radio link design is directly linked to a low transmission power. This, however, cannot be concluded from a simple analysis (see baseline model below). For instance, we may feel tempted to decrease the data rate in order to reduce the required transmission power. This in turn increases the transmission time, and the total energy consumed does not change.

Naturally, this question requires a look to the whole transmitter circuitry (which is a little outside my field). After some research, I noted that there must be a trade off regarding clock rate (high-precision and high-rate is payed with extra power) and amplifier efficiency (lower input powers are usually less efficient). But this suggests that minimising power is not always the best option.

The baseline model (optional)

Consider a wireless radio system where a transmitter $A$ wants to communicate with a receiver $B$, at distance $d$. The total attenuation (due to to signal propagation) can be simply modelled as $L=d^\alpha$, where typically $2 \leq \alpha \leq 4$.

The carrier-to-noise density ratio $\frac{C}{N_0}$ at $B$'s receiver input is

\begin{equation}
 \frac{C}{N_0} = \frac{P_t G_A G_B}{d^{\alpha}N_0},
\end{equation}

where $P_t$ is the power at the output of $A$. When designing these systems, there is usually a performance constraint expressed as a maximum tolerated bit error rate, which turns into a minimum $E_b/N_0 = \gamma_0$ requirement.

So in order to achieve this $\gamma_0$, it seems that there are only two system design parameters to play with (assuming all other parameters are fixed): transmission power and data rate, i.e.

\begin{equation}
 \gamma_0 \propto \Big( \frac{P_t}{R_b} \Big).
\end{equation}

Thus, we observe that increasing the transmission power has the same effect that decreasing the data rate,  in terms of $E_b/N_0$.

Now consider that the total amount of data to be delivered is $L_b$. The total energy consumed by $A$ to deliver $L_b$ bits is therefore

\begin{equation}
 E_t \propto  P_t\frac{L}{R_b}=\gamma_0 N_0 L,
\end{equation}

which does not depend on the operating power neither on the data rate.

Note that there are a host of practical secondary effects of increasing transmission power that do not scale with potentially less required transmission time, such as: battery impedance, voltage converter startup time, oscillator stabilization and carrier synchronization time, leakage or resistive losses due to larger required transistor sizes and/or higher voltage capacitors, more EMI shielding due to any higher digital edge rates required for processing, antenna shielding/spacing to limit human RF exposure, larger required heat spreaders/sinks, and etc.  The goal is often to minimize, not transmission power, but all these secondary power losses or inefficiencies, and/or their associated device costs, assuming some lower bandwidth or latency of delivering messages meets design goal targets.


What is the theoretical bit error rate for QPSK Matlab simulation when constellation points are transmitted (no symbol duration)? I'm simulating a QPSK system in Matlab and I'm trying to calculate the theoretical bit error rate to compare it with the simulated ber. The formula that I'm aware of involves the term Eb/N0, but how do I calculate it, knowing just the SNR and the energy per bit? I seem to need to know the bandwidth of the signal but since I'm just sending constellation points with no carrier (so just complex numbers), there is no symbol duration, so I can't calculate N0 from N0 = Pn/B.

Any help would be greatly appreciated.

TL;DR The theoretical BER of an ideal QPSK system is $\displaystyle Q\left(\sqrt{\frac{2E_b}{N_0}}\right)$ where $E_b$ is the energy per bit and $\frac{N_0}{2}$ is the two-sided power spectral density of the additive white Gaussian noise channel.

The receiver in a standard QPSK system is essentially two BPSK receivers using phase-orthogonal carriers. With a standard matched filter, the signal output in each BPSK receiver has value $E_b$ while the noise variance is $\sigma^2 = E_b N_0/2$.  Putting in the complex-value stuff, the signal output in the QPSK system is a complex number taking on values $(\pm E_b \pm jE_b)$ perturbed by complex Gaussian noise $N_i+jN_Q$ whose real and imaginary parts are independent zero-mean Gaussian random variables with variance $E_b N_0/2$. The decisions that each BPSK receiver makes are independent (because the noise in the receivers are independent), and since each receiver needs to decide between $E_b+N$ and $-E_b+N$, each receiver makes errors with probability
$\displaystyle Q\left(\frac{E_b}{\sqrt{E_bN_0/2}}\right) = Q\left(\sqrt{\frac{2E_b}{N_0}}\right)$. Notice the complete lack of information about the symbol duration $T$ or the channel bandwidth $B$.

"But, but, but..." you splutter, "All this is nonsense.  My receiver
has signal output $(\pm A \pm jA)$ and the noise variance is measured to be $(12+j12)$, not this $E_b$ stuff that you keep on yammering about, and what I want to know is how do I get $E_b$ and $N_0$ when I simulate the system!"

Not to worry. The BPSK receiver is a linear receiver (except for the sampler and decision device), and so the output that you are measuring is really some constant $K$ times $E_b$, and the noise variance $12$ is just
$K^2E_bN_0/2$. So, when you run your simulations, generating gazillions of complex numbers $(\pm A \pm jA)$, perturbing them with complex Gaussian noise of variance $(12+j12)$, your simulated BER should be quite close to $$Q\left(\frac{A}{\sqrt{12}}\right) = Q\left(\frac{KE_b}{\sqrt{K^2E_bN_0/2}}\right)=Q\left(\sqrt{\frac{2E_b}{N_0}}\right).$$
In short, you don't need to know the value of $K$ or to calculate the value of $E_b$ or $N_0$ explicitly; just use your measured values and go from there; it all comes out in the wash.


Converting Analog Filter into Digital Filter, why Bilinear Transform? What's the advantage of using the Bilinear Transform?
$$H_d(z) = H_c(s)\bigg|_{s=\frac{2}{T_s}\frac{z-1}{z+1}}$$
When you can just use this equation?
$$H_d(\omega) = H_c(\Omega)\bigg|_{\Omega=\omega/T_s}$$
In other words, why does bilinear transform exist? These two equations look almost the same to me...What are the tradeoffs between using one equation verses the other?  Is there a case where you would use one verses the other?

$$\left . H_d(z) = H_c(s) \right |_{s = \frac{2}{T\,s}\frac{z-1}{z+1}}$$ describes a transfer function in the $z$ domain that you can easily translate into a difference equation and realize in software.
$$\left . H_d(\omega) = H_c(\Omega) \right |_{\Omega = \frac{\omega}{T\,s}}$$ describes an idealized frequency response that you would like $H_d$ to have when you are done realizing it physically in software.
They're different.  Note particularly the "would like to have" -- any translation from the continuous-time domain to the discrete-time domain is an approximation; part of your job is to make sure it's both practically realizable and good enough.
Note that there are other ways of approximating $H_c$ with some $H_d$ -- the bilinear transform is just one way.  It has a lot of currency because it's conceptually simple, and works pretty well.  It also has a lot of currency because it's easy to do with pencil, paper, and a slide rule -- today, there's numerical optimization techniques that can get you closer to some desired filtering goal, for less work (but -- I can never remember the search terms :( )


How to find the best Q format representation? I'm currently working in a project in which I use the coefficients of a IIR Elliptical digital filter (9th order) and the voice signal from a recording of my cellphone but I'm having issues finding the best Q format to work its implemenation, any suggestions? Till now I've been analysing which are the maximum and minimum values of those two sources in orden to determine which would be the biggest or smallest result from their product, is this a good approach to solve this problem? How can anyone determine the best Q format to work in a system?

You need to calculate the transfer function from your input to each individual state variable.

This depends A LOT on how you implement your filter: I strongly recommend splitting it in second order sections and using either Direct Form I or transposed Form II for each section. Section order and pole/zero pairing can also make a big differences. Unless you have a REALLY good reason not to, follow this recipe: https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.zpk2sos.html

Once you have your filter topology sorted out you need to calculate all the transfer functions from the input to the state variables and look at the  frequency responses. At this point you can also play around with the individual section gains to make this as smooth as possible and manage outliers and peak frequencies.

Then you can pick the Q based on the maximum gain over all frequencies in each transfer function. Add another 3 dB of headroom for transients (sine waves are NOT the worst case signal).


Effects of DC offset on a signal What are some of the reasons why DC offset is considered undesirable for a timeseries signal other than wastage of power during transmission ? Specifically, How does it cause spectral leakage and out of band emissions ?

In a mathematical sense it doesn't but (almost) any nonlinear interaction, say modulation or demodulation, will cause additional signals, leakage/out of band, signals.
The point is that when you do the evaluation you must include the modulator or demodulator; as well as things like transformers which will be biased off center and thus generate extra distortion.  It's the non-linearities in the overall signall path that induce extra signals.  By definition of the Fourier/Laplace transform superposition and linearity properties a linear system will not generate extra signals.


Spatial Audio (Ambisonic) and calculating sound source position I have 4ch Spatial Audio (Ambisonic) WAV file and I want to calculate sound source position and put it on 2D plane.

What I know so far is the following:

WAV file consists of signals that support of 360 video playback

Channel 1 = omnidirectional sphere (named W: $w(t) = p(t)$)
Channel 2 = left-right dipole (named Y: $y(t) = \sqrt{2}p(t)\sin(a)\cos(b)$)
Channel 3 - has no data
Channel 4 = front-back dipole (named X: $x(t) = \sqrt{2}p(t)\cos(a)\cos(b)$)

How to calculate sound source position at any point of time using these data?

Edit:

calculation like this gives mess:

theta = Math.Atan2(y, x);
px = r * Math.Cos(theta);
py = r * Math.Sin(theta);

Edit 2:

So imagine we have Spatial Audio/Ambisonic/B-Format WAV file
and we need to draw/visualize a microphone in the center of screen and a sound source around the microphone at any time reading WAV file

Edit 3:

We use H2n SPATIAL AUDIO mode that gives 1 WAV file with 4 Channel
Using the suggested wiki-formula here is our result for first 1000 ticks using MatLab (does not look proper):

I know nothing about Ambisonic decoding, but the Wikipedia page suggests:
$$
p(\theta_n) = w(t) + x(t) \cos(\theta_n) + y(t) \sin(\theta_n)
$$
Perhaps you could form that and vary $\theta_n$ to find the direction of highest power?

Note: I realize that's not what the page is really saying. It's saying that this is the signal to feed to a speaker at $\theta_n$... but that seems similar to finding the directionality of the source.  Also note that it does not find the distance of the source, just the $x-y$ plane angle.


Reconstructing Signal From Its Cyclic Autocorrelation Can a signal be reconstructed from its cyclic autocorrelation? Specifically, if we know

$$
R^{\alpha}(\tau) = \int{x(t)x^{\ast}(t-\tau)e^{-j2\pi\alpha t}\mathrm{d}t},
$$

can we reconstruct $x(t)\in\mathbb{C}^{N}$? It seems like the answer should be yes, but it's not obvious to me how to do it. I know there has been work on reconstructing a signal from its autocorrelation, and these algorithms of course require some other knowledge about the signal (e.g., Discrete signal reconstruction from its autocorrelation function and one sample). I guess I am hoping knowing the cyclic autocorrelation will at least relax the required assumptions. Moreover, how does the problem change if we also had the conjugate version:

$$
\bar{R}^{\alpha}(\tau) = \int{x(t)x(t-\tau)e^{-j2\pi\alpha t}\mathrm{d}t}
$$

Edit: Let me try rewording the question. I understand that autocorrelation is not an isomorphism. Is this also true of cyclic autocorrelations? Specifically, can we recover $x\in\mathbb{C}^N$ from

$$
R(\alpha,m)=\sum_{n=0}^{N-1}x(n)x^{\ast}(n-m)e^{-j2\pi n \alpha}.
$$

If we only have knowledge for $\alpha=0$, then we have just the regular autocorrelation and the answer is in general no. However, does knowing $R$ for $\alpha=0,1/N,2/N,\dots$ enable one to recover $x$ or at least reduce the required assumptions necessary for doing so? If not, what if the conjugate version is also available? Can we recover $x$ then?

No, you cannot reconstruct the original signal from the cyclic autocorrelation. The fundamental reason is that it results from an averaging operation. Like the autocorrelation and PSD, the cyclic autocorrelation is a many-to-one functional (you already noted that I think). For example, all BPSK signals with identical timing parameters (symbol clock phase and carrier phase), identical pulse functions, and IID symbols have the exact same cyclic autocorrelation functions. Yet there are an infinite number of such signals obtained by all the different sequences of transmitted bits. In general, there is no single unique signal that gives rise to a particular cyclic autocorrelation.

I note that none of the equations you wrote are the cyclic autocorrelation. In general, it is the infinite time average of the second-order lag product of the signal multiplied by a complex sine wave. Aside from the missing scale factor $1/N$, your third equation is an estimate of the cyclic autocorrelation for discrete-time signals.

There are special cases, though, such as periodic signals. Such signals are non-random. The time-varying autocorrelation function $R_s(t, \tau)$ for such signals is identically equal to the second-order lag product itself $s(t+\tau/2)s^*(t-\tau/2)$. Since the time-varying autocorrelation is composed of the sum of sine-wave-weighted cyclic autocorrelations, one might recover $s(t)$ from the set of cyclic autocorrelations, to within a scale factor of $\pm 1$.

The cyclic autocorrelation is sensitive to phase in the sense that the cyclic autocorrelation for $s(t)$ differs from that for $s(t-D)$ in general. This is a basis for synchronization algorithms.

I'm curious why it is obvious to you that there "should be a tremendous amount of redundancy in the cyclic autocorrelation." Redundant what?


On the symmetry of a $2$-dimensional discrete-time signal Can we distribute the minus sign as follows?

$$-h[n_1,n_2] = h[-n_1,-n_2]$$

What you have displayed is true for which it's odd symmetric for only one of the variables $n1,n2$ and even-symmetric for th other.

For 1D signals, an even-symmetric signal has the property that
$$ h[-n] = h[n] $$

and an odd-symmetric signal has the property as:
$$ h[-n] = - h[n] $$

For 2D signals, like $h[n_1,n_2]$, more symmetries are possible:

For even signals, the signal with the property:
$$ h[-n_1,n_2] = h[n_1,n_2]$$
is even-symmetric only in $n_1$ variable, whereas the signal with the property
$$ h[n_1,-n_2] = h[n_1,-n_2]$$
is even-symmetric only in $n_2$ variable. Consequently the signal withe property that
$$ h[-n_1,-n_2] = h[n_1,-n_2] = h[-n_1,n_2] = h[n_1,n_2]$$
is even-symmetric in both $n_1$ and $n_2$.

For odd signals, the signal with the property:
$$ h[-n_1,n_2] = -h[n_1,n_2]$$
is odd-symmetric only in $n_1$ variable, whereas the signal with the property
$$ h[n_1,-n_2] = -h[n_1,n_2]$$
is even-symmetric in $n_2$ variable. Consequently
$$ h[-n_1,-n_2] = -h[n_1,-n_2] = -h[-n_1,n_2] = h[n_1,n_2]$$
is odd-symmetric in both $n_1$ and $n_2$.

So the property you have displayed
$$ h[-n_1,-n_2] = -h[n_1,n_2]$$

implies an even-symmetry on $n_1$ and odd-symmetry on $n_2$ OR an odd-symmetry on $n_1$ and even-symmetry on $n_2$


Differences between filtering and polynomial regression smoothing? What are the differences between classical low-pass filtering (with an IIR or FIR), and "smoothing" by localized Nth degree polynomial regression and/or interpolation (in the case of upsampling), specifically in the case where N is greater than 1 but less than the local number of points used in the regression fit.

Both low pass filtering and polynomial regression smoothing could be seen as approximations of a function. However, the means of doing this are different. The key question to ask here is "Can you do one in terms of the other?" and the short answer is "not always", for reasons that are explained below.

When smoothing by filtering the key operation is convolution where $y(n)=x(n)*h(n)$, which in the frequency domain translates to $y=F^{-1}(F(x)F(h))$ where $F$ denotes the Discrete Fourier Transform (and $F^{-1}$ the inverse). The Discrete Fourier Transform (e.g. $F(x)$) offers an approximation of $x$ as a sum of trigonometric functions. When $h$ is a low pass filter, a smaller number of low frequency components are retained and the abrupt changes in $x$ are smoothed out. This sets low-pass filtering in the context of function approximation by using trigonometric functions as the basis functions, but it is worth revisiting the convolution formula to note that when filtering, y(n) (the output of the filter) depends on $x(n)$ as well as a weighted sum of past samples of $x$ (the weighting here determined by the "shape" of $h$). (similar considerations hold for IIR filters of course with the addition of past values of $y(n)$ as well)

When smoothing by some n-degree polynomial though, the output of the interpolant depends only on $x(n)$ and a mixture of (different) basis functions (also called monomials). What are these different basis functions? It's a constant ($a_0x^0$), a line ($a_1x$), a parabola ($a_2x^2$) and so on (please refer to this for a nice illustration). Usually though, when dealing with equi-distant samples in time and for reasons to do with accuracy, what is used is Newton's form of the polynomial. The reason i am citing this is because through that it is easy to see that when performing linear interpolation you could construct a filter kernel that returns a linearly weighted sum of available samples, just as a low order interpolation polynomial would use "lines" to interpolate between two samples. But at higher degrees, the two approximation methods would return different results (due to the differences in the basis functions).

As i wrote above, not taking into account past values of $x(n)$ is not strict. This is a subtle point. Because usually, when building a polynomial the values outside the given interval ("past" and "future" of a signal) are not considered. It is however possible to include these by fixing the derivatives at the edges of the interval. And if this is done repeatedly (like a non-overlapping sliding window) then effectively, the "past samples" of x(n) would be taken into account. (This is the trick that splines use and in-fact there is a convolution expression for bicubic interpolation. However, please note here that the interpretation of $x$ is different when talking about splines -note the point about normalisation-)

The reason for using filtering as interpolation some times, say for instance in the case of "Sinc Interpolation", is because it also makes sense from a physical point of view. The idealised representation of a band-limited system (e.g. a (linear) amplifier or lens in an optical system) in time domain is the sinc pulse. The frequency domain representation of a sinc pulse is a rectangle "pulse". Therefore, with very few assumptions we expect a missing value to be more or less near its neighbours (of course, within limits). If this was performed with some n-order polynomial (for higher n) then in a way we "fix" the way that a missing value is related to its neighbours which might not always be realistic (why should the sound-pressure values of a wave-front hitting a microphone be fixed to have the shape of a $x^3$ for example? It puts an assumption on how the sound-source behaves which might not always be true. Please note that i do not imply any suitability of an interpolation scheme from a psychophysics point of view here, which involves the processing of the brain (see Lanczos resampling for example). I am strictly speaking about constraints imposed by interpolation when one tries to "guess" objectively missing values.

There is no universal "best method", it pretty much depends on the interpolation problem you are faced with.

I hope this helps.

P.S. (The artifacts generated by each of the two approximation methods are different as well, see for example the Gibbs Phenomenon and overfitting, although overfitting is "at the other side" of your question.)


Record and playback has echo/feedback I didn't mess around with the code here, it was straight from their website, I just want to see if it would work. However, whenever I run it, it has feedback in it, sometimes I can't hear anything at all. It is not the speaker because I can record and play things separately. Here is the code if you want.

I tried changing the frames per buffer to 5000 and it lost the feedback and started an echo. I also tried to have it record and then play back later, but it said that the input had overloaded.

NEW DEVELOPMENT: the feedback is the echo, for some reason when the frames per buffer or the CHUNK is decently high, there is an echo, when it is like the code below, the echo compounds and speeds up to result in feedback, multiple echos playing at once to become a high shrill noise.

PLEASE HELP!!!

import pyaudio

CHUNK = 1024
WIDTH = 2
CHANNELS = 2
RATE = 44100
RECORD_SECONDS = 5

p = pyaudio.PyAudio()

stream = p.open(format=p.get_format_from_width(WIDTH),
                channels=CHANNELS,
                rate=RATE,
                input=True,
                output=True,
                frames_per_buffer=CHUNK)

print("* recording")

for i in range(0, int(RATE / CHUNK * RECORD_SECONDS)):
    data = stream.read(CHUNK)
    stream.write(data, CHUNK)

print("* done")

stream.stop_stream()
stream.close()

p.terminate()

Feedback squeal or echo is pretty much the expected result of the given code.
This program just takes audio in from the microphone and sends it back out through the speaker.  If the speaker is loud enough, the audio the microphone picks up from the speaker will be louder than what was sent to the speaker - with each go around, it gets louder.  This produces the infamous feedback squeal.  If the delay is long enough, you an echo that gets louder and louder instead of a squeal.

You may also have a problem with your sound card configuration.  Some sound cards allow you to send the microphone input straight to the speaker output so that you can use headphones to monitor the recoring live.  If that is active on your sound card, then that may contribute to the feedback.  This will more likely cause a high pitched squeal that is independent of the length of the buffer.

To break the feedback, do the following:

Check that your sound card isn't routing the microphone straight to the speaker. (Stop your program and talk into the microphone.  If you hear your voice coming out of the speaker then you need to fix this in your sound card mixer.)
User headphones instead of a speaker.
Turn the speaker volume down or lower the microphone sensitivity (both can be done through the mixer, speaker can also be done on the speaker with the volume control.)


Examples where non-linear phase filters are used I stumbled apon this old question:
Why is a linear phase important?

There the explanation of why linear-phase processing (filtering) is important is very clear. Also the effects on waveforms due to linear and non-linear phase processing are very well illustrated. My question is the following: why then would someone use non-linear phase filters at all? Can someone give examples where one does not care at all about phase? Are there any advantages in using non-linear phase filters and why is not everybody using only linear-phase signal processing?

Two examples:

Audio: The human ear isn't very sensitive to nonlinear phase (probably because the world filters sound with nonlinear phase filters).  It is, however, sensitive to percussive sounds that "pre-ring"*; i.e., that start making sound before the main "bang".

So for high-fi audio, one often wants to use minimum phase filters (which tend to have nonlinear phase response).  In particular, the use of IIR filters is common (or at least commonly advocated).

Control Systems: In a control system, absolute phase shift (or delay) rather than nonlinear phase shift is the enemy.  FIR filters that have make nice pretty-looking output pulses in response to input pulses make systems that cannot be tuned to the highest possible bandwidth without oscillating -- and systems with pure delay (as you get from a FIR filter) are particularly hard to tune.

So with exceptions (i.e., comb filters can be used to advantage, if you're careful) control engineers want any filtering done with IIR filters, and they tend to want that to be kept to a minimum.

* Or at least there are audio equipment designers that believe this -- I'm parroting things I've read.


Measurement delay versus control delay For a typical conventional control loop (no delay, measurement transfer function is 1) the closed-loop transfer function is
$$H(s) = \frac{G_{comp}G_{process}}{1+G_{comp}G_{process}}$$
With an actuator delay of $T$, the transfer function is as follows.
$$H(s) = \frac{G_{comp}G_{process}e^{-Ts} }{1+G_{comp}G_{process}e^{-Ts} }$$
However, suppose that the delay is in the measurement, i.e. in the feedback loop, the transfer function would then be as follows.
$$H(s) = \frac{G_{comp}G_{process}}{1+G_{comp}G_{process}e^{-Ts}}$$
I have the following two questions.

How do I analyze the stability of the last case? As far as I understand, I cannot use the margin and allmargin Matlab functions as the delay is not in the open loop transfer function ($G_{comp}G_{process}$).
Now assume that I'm stuck with a total delay of $100 \; \mu s$. Say, I can split the delay, $50 \;\mu s$ in the measurement and $50 \; \mu s$ in the actuator. Can I improve the performance, somewhat, of my controller compared to having an actuator delay of $100 \; \mu s$?

Edit Corrected the 2nd transfer function, with the actuator delay only. Thanks to @TimWescott for finding my mistake.

1 - How do I analyze the stability of the last case? AFAIK, I cannot use the margin and allmargin Matlab functions as the delay is not in the open loop transfer function ($G_{comp}G_{process}$).

In both cases the open-loop transfer function is $1+G_{comp}G_{process}e^{-Ts}$

The open-loop transfer function is the transfer function from the output of the summing junction, back to the input of the summing junction.  It includes the effect of all of the blocks -- controller, actuator, plant, and measurement.

You're being mislead by the fact that the closed loop transfer function doesn't include the delay in the numerator -- but it still does.

As to how to make that work with the Matlab functions -- I dunno.  Can they take fixed delays, or must they be presented with rational transfer functions?  If they can only be given rational transfer functions, then you need to calculate the frequency response of the open-loop transfer function and make a Bode or a Nyquist plot and assess stability margins yourself.

  2 - Now assume I'm stuck with a total delay of 100 µs. Say I can split the delay, 50 µs in the measurement delay and 50 µs and the actuator delay. Could I improve somewhat the performance of my controller compared to having an actuator delay of 100 µs?

In this regard, having less actuator delay is fundamentally better, because you have less delay from command to output.  Actually designing a controller that can take advantage of that may be a challenge, however.  Certainly the response to a disturbance cannot be made better -- only the response to a command.  In general pure delays like that just make life hard, no matter what.