Negative instantaneous frequency with hilbert transform using scipy hilbert When I produce a chirp signal with a duration of 2 sec, sampling frequency of 200Hz, with frequencies ranging from 20Hz-100Hz, I get negative instantaneous frequencies using the scipy(Python) Hilbert transform. Here is my signal code:
    self.duration = 2
    self.fs = 200.0
    samples = int(self.fs * self.duration)
    self.time = np.arange(samples) / self.fs
    carrier = chirp(self.time, 20, self.time[-1], 100)
    self.signal = carrier * (
        1 + 0.5 * np.cos(10 * 2.0 * np.pi * self.time)
    )

Here is my Hilbert transform code:
demean = data - np.mean(data)
hil = hilbert(demean)
inst_amplitude = np.abs(hil)
inst_phase = np.unwrap(np.angle(hil))
inst_frequency = np.diff(inst_phase) / (2 * np.pi) * sampling_frequency
inst_frequency = np.insert(inst_frequency, 0, 0)
regen_carrier = np.cos(inst_phase)

Here is the plot of the Hilbert properties:

Why are there negative instantaneous frequencies? I have been Googling the issue with zero luck finding an answer. Any input would be greatly appreciated.
UPDATE: I applied the exact same Hilbert transform code as above to a "real life" dataset with a sampling frequency of 800Hz. I am getting negative instantaneous frequencies in the middle of the data, not just at the edges. I don't think I am allowed to show the data as it is proprietary, so sorry on that front. Any insight as to why I am seeing negative frequencies in the middle of my data and not just at the edges?
UPDATE #2 I went ahead and generated a synthetic chirp signal with multiple frequency sweeps using scipy.signal.sweep_poly and I am getting negative instantaneous frequencies in the middle of my data.
Here is the code for my new signal (duration of 50sec and sampling frequency of 100):
    self.duration = 50
    self.fs = 100.0
    samples = int(self.fs * self.duration)
    self.time = np.arange(samples) / self.fs
    p = np.poly1d([0.025, -0.36, 1.25, 2.0])
    self.signal = sweep_poly(self.time, p)

Here is the plot of the new signal's Hilbert properties:

Here's the same image, but zoomed in where there are negative instantaneous frequencies in the middle of my data:

Why are there negative instantaneous frequencies in the middle of my data? There are definitely no edges in the middle of my data, so it's not an edge effect.
UPDATE #3 I added random noise to the exact same signal using the following code:
    self.duration = 50
    self.fs = 100.0
    samples = int(self.fs * self.duration)
    self.time = np.arange(samples) / self.fs
    p = np.poly1d([0.025, -0.36, 1.25, 2.0])
    carrier = sweep_poly(self.time, p)  # Carrier signal
    noise = np.random.normal(0, 1, carrier.shape)
    self.signal = carrier + noise

Here is the plot of the noisy Hilbert transform properties:

Even more negative frequencies! Why is this happening?
UPDATE #4 I decided to multiply 2 sine functions with two different constant frequencies together to see what I get. There is no noise in this instance. Here is my code:
    self.duration = 2
    self.fs = 100
    samples = int(self.fs * self.duration)
    self.time = np.arange(samples) / self.fs
    self.signal1 = np.sin(1 * 2.0 * np.pi * self.time)
    self.signal2 = np.sin(9 * 2.0 * np.pi * self.time)
    self.signal = self.signal1 * self.signal2

Here is the plot of the resulting Hilbert properties:

If I change the sampling frequency, the magnitude and sign of the instantaneous frequency spikes change! Here are a few examples (the sampling frequency used is in the title of the image):

Something about the phase calculation is unstable and the sampling frequency seems to alter things. No I am even more curious as to what is happening.
Thanks!

It looks like artifacts due to the derivative. I used this code in Octave:
fs=300;
t=[0:1/fs:2];
c=chirp(t, 20, 2, 100);
s=c.*(1+0.5*cos(2*pi*10*t));
h=hilbert(s);
m=abs(h);
a=diff(unwrap(arg(h)))/2/pi*fs;
plot(m,"",a)

and this is what comes out:

I also tested the equivalent of this in LTspice and, without any form of unwrapping, this is the result:

where V(mag) is the result of hypot() (m), V(signal) is s, V(freq) is the frequency for the FM, and V(demod) is a. All are delayed (for matching I/O) by 23/2/200 which represents the order of the Hilbert transformer (Hamming window), divided by two (group delay) and by the sampling frequency (200).
My guess: don't take for granted the result of the differentiator and try to take a finer comb to it (is this a pun? not intended).


BW of a modulated signal I come across different types of bandwidth. In digital modulation should it be null to null always?
when do we use 3 db bandwidth?

It's whatever makes the most sense for the problem at hand.  Which can get confusing, and the same signal may have multiple definitions of "bandwidth" -- sometimes even in the same document, when, for example, one is reading up on or designing a communications system.
Useful Bandwidth
A signal's user is interested in the practically useful bandwidth of the signal, i.e., if it's analog can I fit comms-quality audio (300-3000Hz) on it, or hi-fidelity audio (20-20000Hz), or video (0-4.5MHz for old-style NTSC, and yes, that's debatable), and if it's digital what symbol rate could I achieve.
That 3dB point you're wondering about is often taken to be a signal's useful bandwidth.
Occupied (or assigned) Bandwidth
How much bandwidth the signal takes up in reality.  This is the amount of bandwidth that you need to assign for the signal in the real world, so that signal A doesn't interfere with signal B when you're using real-world equipment.  Typically occupied bandwidth is considered to be the bandwidth within which some large percentage (99%, 99.9%, etc.) of the total energy is confined, and in regulatory documents its usually defined in a way that makes it easy to test with reasonably-expensive instruments (i.e., a high-quality but not necessarily lab-grade spectrum analyzer).
Equivalent Noise Bandwidth
This is more a property of filters or systems, but the noise bandwidth of a filter is -- assuming the filter gain is normalized to 1 -- the ratio of the total noise power out of the filter to the total noise power in.  So, for example, the equivalent noise bandwidth of a 1$^{st}$-order lowpass filter with a 3dB point at $f_0$ is $\frac{\pi}{2}f_0$.
Whatever
Sometimes, in the context of your work, you'll come across some other bandwidth that you wish to define (such as null to null).  In that case, you state so clearly, and you use that as "bandwidth".  So if you run across some screwy bandwidth definition, accept that the author has their reasons, understand it, and roll with it.


Verifying the result from lowpass filter I have a event record from acceleration censer.
I'm trying to remove gravity from the event record by low/high pass filter as following site shows.
https://developer.android.com/reference/android/hardware/SensorEvent#values
As the first step I tried to extract Gravity by using following formula as a low pass filter.
G[t] = 0.1*G[t-1]+(1-0.1)*a[t]

G[t] is extracted gravity as the result.
G[t-1] is previous result .
a[t] is censer value in current row of event record.

I calculated it on Excel sheet for each X,Y,Z axis.
And I'd like to know how to verify the result.
I thought I can verify by summarizing  G value for X,Y,Z from low pass filter will be 1023 based on following picture.

Green is X
Blue is Y
Red is Gravity
G only weighs on Y at the first.
After rotate X-Y,G will weigh on both X and Y.
There will be G on X and G on Y
However, Gx+Gy+Gz is not 1023 during censing so far.
Here is  top 8 line of the data.
time-x  x   time y  y   time-z  z   time-zz zz
0   904 0   -20 0   504 0   -320
0.006   892 0.004   -16 0.006   504
0.008   888 0.015   -4  0.015   488
0.02    888 0.018   8   0.02    492
0.022   892 0.022   16  0.022   492
0.031   876 0.029   20  0.032   500
0.034   876 0.034   4   0.038   496
0.049   884 0.051   4   0.056   488

To eliminate gravity you need a high pass filter. From your source

This can be achieved by applying a high-pass filter. Conversely, a low-pass filter can be used to isolate the force of gravity.

You implemented a low pass filter. Instead try
$$G[n] = p \cdot G[n-1] + 0.5 \cdot (1-p) \cdot (a[n]-a[n-1])$$
where p is the location of your pole. I would start with $p = 0.9$

And I'd like to know how to verify the result.

Same you should do with any algorithm

Write down the math of the algorithm
Select some suitable test signals where you can pre-calculated the result
Run the test signals through your implementation and check whether the results match.
Verify all other quantitative requirements

So in your case, you should run a constant signal through and check whether the output is zero and you should run a sine wave at different frequencies to it and verify that the output is indeed a sine wave at the same frequency and that the gain at the corner frequency is $\sqrt{1/2}$


Compressive sensing: numerical generation of RIP matrices The restricted isometry property (RIP) states that:
\begin{equation}
  (1-\delta_K)||x||_2^2 \le ||A x||_2^2 \le (1+\delta_K)||x||_2^2
\end{equation}
for any $K$-sparse vector $x$ of length $N$. The corresponding restricted isometry constant is $\delta_K$, $0 < \delta_K < 1$.

I would like to perform a numerical experiment based this expression to produce a RIP-compliant matrix $A$ (Matlab code is supplied below).
The matrix is constructed with Gaussian-distributed RV's (zero mean, variance $1/\sqrt{M}$). Having specified $N$ and $K$, the number of measurements $M$ is obtained with $M \geq K \log(N/K)$. This is a hit-or-miss procedure, I let the code run until the isometry constant falls between zero and one.
Is this the proper way to generate $A$? Am I missing anything? Further analysis on $A$ (eigenvalues, etc.) doesn't seem to behave properly...

function [delta] = RIP_numerical(n,s)

% delta: restricted isometry constant; n: # columns; s: sparsity

% choose m (# measurements/rows) based on measurement criterion

m = ceil(s*log(n/s))

flag = 0; incr = 0; % initialize computing parameters

while flag ~= 1  % run until condition met

    incr = incr + 1
    % generate random, m x n normalized matrix N(0,1/m)
    A = (1/sqrt(m))*randn(m,n);
    sumA = sum(A);
    for k = 1:n
        A(:,k) = A(:,k)/sum(A(:,k));
    end

    % generate vector of length n, sparsity k
    x = zeros(n,1);
    list = randperm(n);
    for l = 1:s
        x(list(l))= randn;
    end

    y = A*x;    % obtain measurements

    norm_x = norm(x,2);  % l_2 norms
    norm_y = norm(y,2);

    delta = 1-(norm_y^2/norm_x^2);  % look at lower bound

    if (delta >0 & delta < 1) % condition met; save vars and kick out of routine.
        flag = 1
        save('RIP_config.mat', 'A','x','n','m','s');
    end
end

You can't prove RIP through numerical exploration of all possible cases. If you are interested in numerical analysis I suggest to use Coherence instead,  however Coherence is not as strong condition as RIP. Note that RIP is a tool which is used for Mathematical proofs.
I explain how to check for RIP, in case your are interested.
As I said Coherence is a weaker condition in comparison with RIP, however, it is numerically feasible to check. Coherence of matrix A (which is sensing matrix multiplied in sparsity basis) is defined as below:

where $a_i, a_j$ are columns of matrix A, and $u_A$ is coherence factor. Its computational complexity is of $O(n^2)$ which in many cases is feasible.
For checking coherence, you simply need to compute the given type of correlation between each two different columns of A matrix and among those values chose the max. Based on Welsh Inequality we'll have:

In CS sampling $M<<N$, so the above relation simplifies to $u_A>sqrt(1/M)$, therefore, to satisfy RIP using coherence, number of measurements ($M$) must be of order of $O(K^2)$ where $k$ is number of non-zeros. This number of measurements is so greater than $Klog(N/K)$.

Another conditions instead of RIP, is stRIP which is much much simpler to prove.

Ref: [A.Amini,"Deterministic Compressed Sensing" PhD Thesis, Sharif University of Technology, Iran, 2011]


Inverse DTFT low pass filter can someone explain why the ideal frequency response of a low pass filter cant be implemented using the inverse discrete time fourier transform.

i understand the concept of why it cant be implemented because a digital filter is infinite and non-causal, but i can't explain it using the equation. can someone help?

looking at this equation for an ideal low pass filter

how can i explain using this that the ideal frequency response cannot be implemented in practise?

Since the impulse response of an ideal low pass filter extends to plus and minus infinity, the corresponding filter is not causal, i.e., in order to implement it (exactly) you would need to know the complete input signal (including its future) to compute the current output value. Furthermore, the sinc function only decays as $1/n$, which makes the filter unstable. In sum, the ideal low pass filter is non-causal and unstable, and for that reason it can't be implemented. Of course, we can approximately implement it by making it causal and stable in such a way that its desirable characteristics are not changed too much for the given application. That's the topic of the field filter design. An little overview for the discrete-time case can be found here.


Derivative Notation as a factor In the figure below, from the book "Principles of Linear Systems and Signals - Lathi", the author uses the notation D to replace the Leibniz's notation $d/dt$, and after that he rewrites the equation, disconnecting the symbol of the derivative, $D$, of its function, considering each as a factor. I don't understand how he did that. Is that mathematically possible?

As stated in the text, it's called operational notation, i.e. $D$ is an operator being applied to a continuous-time signal:

$$D^mx(t)=\frac{d^mx(t)}{dt^m}$$

In discrete time you have the operator $T$, defined by

$$T^mx[n]=x[n+m]$$

(see e.g. here).

What is often done by people is that they use the complex Laplace transform variable $s$ instead of $D$, or the $\mathcal{Z}$-transform variable $z$ instead of $T$. This is what Robert Bristow-Johnson referred to in his comment as 'sloppy abuse of notation'. I would even go further and would call it simply wrong. The use of different letters for $D$ and $T$ may seem an unimportant triviality, but the real problem here is that you mix different domains (Laplace/$\mathcal{Z}$ transform domain and time domain), which can be a great source of misunderstandings and errors.

I personally do not see a great advantage in the use of this type of operational notation for engineering applications. It rather makes things more confusing for beginners because they have a tendency to mix domains anyway, and operational notation seems to suggest that this is legitimate.


How does MATLAB's $\tt upfirdn$ select extra points when downsampling? I had a signal 512 samples long and I used upfirdn to filter and downsample it by a factor of 2. I thought this would result in a 256 sample long variable, but instead it comes out as 261.

The MATLAB help says this:

  The length of the output, yout, is
  ceil(((length(xin)-1)*p+length(h))/q)

(where xin is the signal, h is the filter, p is the upsample ratio, and q is the downsample ratio) which, if you work it out (filter length is 10), does in fact come out to 261.

What I want to know is how it selects points when downsampling so that I can reproduce the upfirdn behavior in another language (long story, not important). Where does that formula come from?

The function performs upsampling, filtering and downsampling. The reason for the extra samples is due to the fir filter delay. If you want to reproduce this behaviour for downsampling only you should do the following steps:

Take your input signal, xin and add length(h)-1 zeros at the end. This is done because of the filter's delay.
Use filter(h,1,xin_padded); to filter the signal.
Downsample the filtered signal.


Impulse invariance vs. DT representation of a CT system: Where is the inconsistency? Suppose you have a continuous-time (CT) system $h_c(t)$, bandlimited to $B$. Your goal is to represent the system as a discrete-time (DT) system $h[n]$, sampled at $f_s \leq 2 B$. Clearly $h[n]$ won't be a faithful representation of $h_c(t)$ because we do not satisfy the sampling theorem. But I am interested in the exact relationship between $h_c(t)$ and $h[n]$ in this case. There are two ways to think about this:
Representing the CT system as DT one
Based on this question (and answer): Representing a continuous LTI system as a discrete one, $h[n]$ represents a low-pass filtered version of $h_c(t)$:
$$
h[n] = \tilde{h}_c(n T_s) ,\\
\tilde{h}_c(t) = \int h_c(\tau) \operatorname{sinc}\left(t-\tau\right) \operatorname{d}\tau
$$
Impulse Invariance
According to https://en.wikipedia.org/wiki/Impulse_invariance, the sampled $h[n]$ corresponds to an aliased $h_c(t)$:
$$
H(e^{-j\omega}) = \frac{1}{T_s} \sum_{k=-\infty}^{\infty} H_c\left(j\frac{\omega}{T_s} + j\frac{2\pi}{T_s}k \right)
$$
I could take the IDTFT:
$$
h[n] = \sum_{k=-\infty}^{\infty} T_s e^{-j2\pi k n} h_c(n T_s)
$$
which is not the same as the low-pass filtered version above.
There no low-pass filtering.
These two approaches seem to give me two different answers on how to relate $h[n]$ and $h_c(t)$ (when sampling theorem is not obeyed).
What is wrong?

The Impulse Invariant method does not promise to represent the frequency response of the  continuous-time system as a lowpass version.  What the Impulse Invariant method does is frequency-alias the frequency response by sliding it by multiples of the sampling frequency and adding up all of the translated copies.


What is the magnitude of the impulse function at t=0? The unit impulse function, in a few textbooks that I have referred, has a value of 0 at t≠0 , and an area of unity (1). The height of the impulse function also tends to infinity at t=0. But since it's height cannot be it's magnitude (as we generally use it's unit area in calculations), what would be the value/amplitude/magnitude of the unit impulse function at t=0?

Before you say that my question has answered itself, and that it's value is one(1) at t=0, note that this is not defined as such in any standard textbooks. This question has been bothering me for quite a while now, so any attempt at elaborating a solution will be highly appreciated.

The "impulse function" is not actually a function but a distribution, and its value at $t=0$ is not defined.

The impulse only really makes sense when multiplied by a signal inside an integral, where it is used to "select" the value of the signal at specific time instants.

Personally, I try to avoid impulses whenever possible. Look for Lapidoth's book on digital communications (available for free on his website) for a complete treatment of Fourier analysis (as relevant for communications) without a single delta. Not many people share my squimishness with deltas, though, and if you're taking a signals course they are pretty much unavoidable.

Note also that, once you go to the discrete time domain, everything becomes nicer: instead of the Dirac delta, you use the Kronecker delta, which is a regular function and has a value of 1 at $n=0$.


Why do digital images use the concept of ‘Stride’? … why are Width and Height not sufficient to work with digital images?

"Stride" is something completely different then width and height of an image; it's a concept related to the logical storage layout (memory addressing) of the image data array on the physical memory accessed by a program, and not to the geometric properties of the image.

Stride means "the logical memory address distance between two successive pixels of the image on a given axis" (it's the value that should be added to the pointer (which points into the buffer of the image data on memory) that would then access the data of the next sample on the chosen image axis)

For example, consider a row-major image with width $w=3$ and height $h=2$:

$$\mathbb R^{2\times 3}\ni A = \begin{bmatrix}
a_{11} & a_{12} & a_{13}  \\
a_{21} & a_{22} & a_{23}  \end{bmatrix}$$

with the internal storage looking like

$$
\begin{matrix}
\text{storage loc.} & \text{pixel} \\
0 & a_{11}\\
1 & a_{12}\\
2 & a_{13}\\
3 & a_{21}\\
4 & a_{22}\\
5 & a_{23}\\
\end{matrix}
$$

Here, the stride between elements of the same row (i.e. along the horizontal axis) is simply $s_x=1$.

The stride along the vertical axis is exactly the width ($s_y=w=3$) of the image here, since the storage is row-major and compact.

Now, it's pretty common to see formats specify a step size for row lengths, because that typically has computational advantages. For example, if we took the same $A$ and said "the row length is padded to multiples of $4$", then the storage format for $A$ would look like
$$
\begin{matrix}
\text{storage loc.} & \text{pixel} \\
0 & a_{11}\\
1 & a_{12}\\
2 & a_{13}\\
3 & \text{padding (don't care)}\\
4 & a_{21}\\
5 & a_{22}\\
6 & a_{23}\\
3 & \text{padding (don't care)}
\end{matrix}
$$

Now, the stride (i.e. how far are two elements successive on the same column apart?) is $s_y=4$, but the image width is still $w=3$, hence $s_y \mathbf \ne w$. I think this illustrates nicely why stride and dimensions are very different things.

This gets even more important for higher dimensional data.

For example, imagine a video. We can think of it as sequence of images like $A$. Let's call it $B$:

$$B = \left(A_i\right)_{i=1,\ldots N}\quad\text,$$

with $A_i \in \mathbb R^{2\times 3}$ being our frames, and thus, $B$ having the first axis signifying "frame number", the second "column", the third "row index",

$$ B \in \mathbb R^{N\times 2\times 3}\,\text.$$

Now, let's again assume the padded storage format.

The row- and column-strides stay the same.

But, what is the storage distance between two pixels with the same position, but in successive frames?

$$\begin{align}
s_t &= h s_y\\
&= 2 \cdot 4\\
&=8
\end{align}
$$

Notice that we had to use height times y-stride here because we've defined our pictures to be row-major. In column-major format, we'd use width times x-stride. (and in time-major format: I'll leave this up for you as exercise in basic understanding of storage formats).


What is fft input size for ofdm? From the Transmitter comes a block of size N+L (where N is the number of subcarriers and L is the cyclic prefix length). Then the block goes through a Channel which adds delay and AWGN, so this process can be described as a convolution of a transmitted block with a channel (the channel length H).

The question: The output of the convolution increases the size of transmitted block (the output from a channel is of lenth [N+L+H]) so how can FFT be performed (of course after removing cyclic prefix)? (because the number of subcarriers must be the size of IFFT/FFT)

You're doing the same length as in the transmitter.

You're ignoring the purpose of the cyclic prefix: it's both to allow the "echoes" from the previous symbol to die off (i.e. it contains H). That's why your CP is always a bit larger than the length of channel you're expecting.


Filter Implementation I am rather new to the world of signal processing, and am struggling to understand a fundamental concept: How are filters actually implemented?

I have read a significant portion of this online book, and scrounged the internet, finding snippets of useful information here and there, but I cannot quite tie it all together.

In brief, I have a vibration profile (~200,000 samples) in the time domain, which I would like to analyze using a particular frequency weighting. This involves a four-step filter, using a (Butterworth) high-pass and low-pass, followed by an a-v transition and an upward step. The analog equations are shown below.

Yet I am still struggling on how to implement them. I have seen mentions of the bilinear transform, used in programs in MATLAB/Python, but such implementations seem to omit the 's' or 'p' variable from the filters shown, typically creating the numerator 'B' and the denominator 'A' from the coefficients of the analog filters.

Other sources suggest taking the Fourier of the filter equations and the signal, multiplying, then applying the inverse Fourier to bring it back into the time domain. Yet I cannot discern how to take the Fourier of an equation.

My same problem occurs with using convolution in the time domain, as I have not understood how to convolve an equation yet.

I am sure I have come near the answer, and I know this is a basic principle, but I cannot seem to comprehend it, and my Googling is off the mark. Any help on this would be appreciated. General advice is also welcome, as I have a lot to learn.

** If important, I am currently working in Python, using the scipy.signal package, with access to MATLAB for testing and evaluation purposes.

Currently I am looking at substituting j2πf for p, then simply running the data through the equations and taking the inverse Fourier to bring it back into the time domain.

Alternatively,  my eyes have just noticed what was staring me in the face: the second equation given for each filter is in terms of 'f', which would suggest I am able to simply run the frequency-domain data through this equation before taking the inverse Fourier to recover the time-domain. If this is a correct realization, this question can be closed (or I can 'answer' my own, now self-evident, question).

There are a few different ways to proceed. They involve doing some algebra but it is more an issue of tedium. The notation of your analog filters is a little non standard, essentially your $p$ is $s$ in most books.

Traditionally, phase wasn't considered important in hearing so standards tended to be specified by magnitude, so there is some ambiguity with respect to the digital implementation's phase response. You have 2 basic simple choices, a digital IIR filter or FIR filter. FIR is likely to be easier but easy and better are typically not the same. FIR can be direct or FFT based.  IIR has a number of choices such as biquads or direct. Floating point math has its own easy/better issues.

One path way is to focus on the terms inside the magnitude brackets, i.e. the complex analog transfer functions. The other path is to focus on the analog magnitude.

So, given the analog transfer functions, you can either use a table of Laplace transforms and deduce the impulse response of each component you listed, or use the bilinear transform to express each component that is in $s$ in terms of $z$.

If you derive the impulse responses, you can sample points of the impulse response and provided a generally exponential decay, use those samples as coefficients of a FIR filter of sufficient length. The FIR filters are then applied sequentially. One nice property of ideal LTI filters, the order that you apply the filters are ideally interchangeable.

If you choose the bilinear transform and do the algebra to derive the discrete time numerators and denominators, you need to be aware that the bilinear transform squeezes the response as it approaches the Nyquist frequency. One needs to plot the discrete time response to asses how the corresponding analog response is distorted. If some modest tweaking results in an acceptable response then, the numerators and denominators need to be factored so that they are compatible with the IIR implementation, such as a cascade of biquads.

If you go with the magnitudes, the task becomes a direct sampling from the analog frequency domain to the discrete time periodic frequency domain. One can choose a linear phase response, and appropriately impose symmetry and define the frequency domain FIR filter. This might not sound like an analog implementation.

The essential challenge is mapping a frequency response in $-\infty < \Omega < \infty$ to $-\pi \le \omega \le \pi$

There are other approaches like state variables, so this is not the entire story.


How can I vectorize the computations for a first-order recursive filter? I have a simple single pole low pass filter (for parameter smoothing) that can be explained by the following formula:

$$
y[n] = (1-a) y[n-1] + a x[n]
$$

The architecture that I'm using has access to single-instruction, multiple-data (SIMD) instructions that can perform multiple vectorized calculations in parallel. I would like to take advantage of this capability, but I'm not sure how to do that for a recursive filter like this. The problem is that every computation needs a previous result.

Assuming that you perform vector operations $M$ elements at a time, you can unroll the difference equation by a factor of $M$ pretty easily for the simple single-pole filter. Assume that you have already calculated all outputs up to $y[n]$. Then, you can calculate the subsequent ones as follows:
$$
\begin{align}
y[n+1] &= (1-a)y[n] + ax[n+1] \\
y[n+2] &= (1-a)y[n+1] + ax[n+2] \\
&= (1-a)((1-a)y[n] + ax[n+1]) + ax[n+2] \\
&= (1-a)^2y[n] + a(1-a)x[n+1] + ax[n+2] \\
&\vdots
\end{align}
$$

In general, you can write $y[n+k]$ as:

$$
y[n+k] = (1-a)^ky[n] + \sum_{i=1}^k a(1-a)^{k-i}x[n+i]
$$

For each sample index $n+k$, this looks like an FIR filter with $k+1$ taps: one tap multiplies the last filter output $y[n]$, and the other $k$ taps multiply the filter inputs $x[n+1], \ldots , x[n+k]$. The nice thing is that the coefficients used for all of these taps can be precalculated, allowing you to unroll the recursive filter into $M$ $M+1$-tap parallel non-recursive filters (these filters compute the output samples $y[n+1], \ldots , y[n+M]$), updating the feedback term every $M$ output samples. So, given an initial condition $y[n]$ (which is assumed to be the last output calculated on the previous
vectorized iteration), you can calculate the next $M$ outputs in parallel.

There are some caveats to this approach:

If $M$ becomes large, then you end up multiplying a bunch of numbers together in order to get the effective FIR coefficients for the unrolled filters. Depending upon your number format and the value of $a$, this could have numerical precision implications.
Also, you don't get an $M$-fold speedup with this approach: you end up calculating $y[n+k]$ with what amounts to a $k$-tap FIR filter. Although you're calculating $M$ outputs in parallel, the fact that you have to do $k$ multiply-accumulate (MAC) operations instead of the simple first-order recursive implementation diminishes some of the benefit to vectorization. The non-vectorized approach uses 2 MACs per output, so you need $2M$ operations to calculate $M$ outputs. The vectorized scheme computes the $M$ outputs at once, requiring $M+1$ MACs in the process. So, the reduction in operations can be expressed as a function of $M$ as:

$$
    R = \frac{M+1}{2M} = \frac{1}{2}\left(1 + \frac{1}{M}\right)
    $$

So, even with $M$ very large, you can get a maximum of a 50% reduction in the amount of computations using this method. For the common values of $M=4$ and $M=8$, the reduction is 37.5% and 43.75%, respectively. However, in addition to pure reduction in the number of operations, you may gain additional performance benefits due to the different memory access pattern employed by the unrolled approach; for the simpler implementation, you may run into delays due to the dependency of each output sample $y[n]$ on the previous sample $y[n-1]$. This effect is obviously very platform-dependent.


algorithm for feature extraction from audio signals I am trying to build a model for speaker identification, and I understand that the first step is to extract the features from the audio signals that are in my database. I have done quite a bit of research and cant find how to do this extraction and to which features. How do I select which features to extract to? what would be the algorithm? Can someone please point me to either an algorithm that I could implement or a sample matlab code?

For speaker recognition, the features you should probably start with are MFCC.

There are libraries offering MFCC extraction modules, such as YAAFE, aubio (C/C++), the MIR toolbox or Dan Ellis' implementation (Matlab) - and of course speech recognition frameworks (Sphinx, HTK) provides MFCC extraction tools.


Has there been any research on time-windowed spectral estimation using more than one window function at the same time? In the Welch's method, a signal is divided into overlapping segments, multiplied with a window function, and transformed to frequency domain via fourier transform.

Why does the method use only one window function? Suppose I would use (say) ten window functions to produce ten spectral estimates. By averaging the estimates, my SNR would be higher than by using just one of the windows?

Congratulations! Sounds like you are about to reinvent the multitaper method, in which a number of orthogonal window functions are used to get independent spectral estimates that are then averaged (with weights). A suitable window function family are the Slepian sequences:


What is the effect of a larger window size on "temporal resolution"? I guess my question boils down to what "temporal resolution" means?
I'm taking a signal processing class right now and we're learning about DFT and windowing at the moment. We've learned that when we increase the window size, $N$, we have greater "frequency resolution"

From the image we see that the more we increase, $N$, the more we get the actual $\Omega$ we want and the less smear there is.
Isn't this all there is to resolution? We're getting the $\Omega$ that we actually want? Is there something else important that I'm being really oblivious to?
My professor mentions that a larger window, $w[n]$, provides less "temporal resolution"... What does that mean? And how is it different than having good frequency resolution?

You have probably heard of the uncertainty principles that arise in physics, with the most famous of all being the Heisenberg's uncertainty principle. But the uncertainty principle applies to signals as well. Just as Heisenberg's uncertainty is described by the following inequality:
$$
\Delta_x\, \Delta_p \ge \dfrac{h}{4\pi}
$$
The above imposes a limit in the level of resolution in the measurement of a particle's position and momentum. In fact their relationship is inverse meaning that as the particle's position measurements precision increases, then its momentum measurement precision decreases or in other words the uncertainty in its momentum measurement increases.
Now in signal processing the goal is to extract information about the signal's structure (if any) usually in the frequency or time domain (it could also be in the spatial domain).
But as you correctly described and analyzed in your graphs the more precise we get in the time domain the more uncertain our estimation of its frequency spectrum becomes. You can see it intuitively in the ripple. The smallest temporal window $N=16$ results in the biggest ripple at higher frequencies, taking also a sizeable proportion in the signal's power spectrum, but as $N$ increases so does the ripple increasing our certainty about the signal's frequency spectrum structure.
$$
\sigma_t\, \sigma_f \ge \dfrac{1}{4\pi}
$$
where $\sigma_t$ and $\sigma_f$ are the standard deviations of the time and frequency domain's estimations, respectively.
This minimum threshold is called the Gabor limit, and the principle is called Gabor's uncertainty.
Note:
  What is meant by "a larger window provides less temporal resolution" is that some of the properties of the signal in the time domain may be less prominent whereas in the frequency domain we would be more certain about its properties and structure.
From a practical perspective this limit is also important because in addition to other factors it could help us choose an appropriate window for our particular application where it could be desirable to extract features from both domains, sometimes we care only for the time domain properties, sometimes for the frequency domain properties and other times for both (also sometimes implications can be made about properties in one domain from properties measured in the other domain).
In conclusion, the more a signal spreads in the time domain, the less it spreads in the frequency domain and vice versa. In analog signals, a classic example of this phenomenon is the a rectangular pulse in the time domain, which Fourier Transform is the sinc function (exhibits spectral leakage).


Matlab BOC (binary offset carrier) I am trying to implement a binary offset carrier (BOC) direct sequence spread spectrum scheme in MATLAB. In that scheme, a code $c(t)$ (a binary vector which repeats after a certain number of zeros and ones) should be calculated by an other square (sub) carrier of the following form sign(sin(1*pi*fs*t) where fs is the subcarrier frequency and finally should be calculated by a carrier frequency having the usual form of a cosine.

I have generated the code needed (a vector of 1023 elements). These 1023 chips (vector’s elements) should be repeated every 1 ms in order for the chipping frequency to be achieved (10Mchips/s). After, that vector (code) should be multiplied by a square subcarrier of 1023 Mhz and finally by the ‘actual’ carrier (a co sinusoidal signal) of 4092 Mhz. All these should be demonstrated via plots.

My first problem is that in order to demonstrate the spectrum of the codes, spectral nulls should be represented at frequencies multiples of 1.023 Mhz. The only way that I could somehow think how to do this is to use massive vectors where are not really supported by my laptop's memory (RAM).

I say "somehow", because the sampling rate should change (for example wherever I have one ‘1’ I should add some more and the same applies to zeros). This procedure ruins the whole purpose of the code because the code has some particular autocorrelation properties which I don’t know if will continue to exist after.

Moreover, if that should reach the 10Mchips how will this could be done by using the particular vector. While I understand that it is easily done by using a random vector of 10*10^6 elements. Should I repeat the 1023 chips 1000 times?

First off- I know you are doing you're best (which is why I am answering) but your question is very hard to read.  Regardless of your English skills, it would be easier to read if you would break it up into paragraphs and throw in some pictures to help show what you're talking about.

  These 1023 chips (vector’s elements) should be repeated every 1 ms in order for the chipping frequency to be achieved (10Mchips/s).

1023 chips repeated every 1 ms means you have $\frac{1023}{.001 s} = $1 Mchips/s, not 10 Mchips/s.

  My first problem is that in the first place in order to demonstrate the spectrum of the codes,spectrul nulls should be represented at frequencies multiples of 1.023 Mhz.The only way that I could somehow think how to do this is to use massive vectors where are not really supported by my laptops memory (RAM).

You should not need a "massive" vector to demonstrate this.  Increasing the number of samples that you FFT just gets you better frequency resolution- it doesn't change the frequency range.  I think that you should be able to see the nulls just fine with several repetitions of the spreading code.

  I say somehow because the sampling rate should change (for example wherever I have one ‘1’ I should add some more and the same applies to zeros).

I have no idea why you say that the sampling rate should change, unless you're talking about going from the unspread rate to the spread rate.


Guess the frequency response of a filter with little information I have a processing unit that does a nice EQ filtering on soundfiles.
The only thing I can do is : I give the original .wav file, and the processing unit gives me back the filtered .wav. I have no other information than this.

I'd like to know the frequency response graph of this processing unit. How to do it ?

1) feed the unit with a "flat spectrum white noise", and see the spectrum of the filtered file? It works but it is not very accurate
2) feed the unit with sinewaves soundfile with frequency = 20, 30, 40, 50, etc... 20 000 Hz.
Then I examine the response for each file... and I can draw a frequency response graph after 10 hours of work :)
3) feed the unit with a frequency sweep ?  and then try to see what's the unit's reponse to this sweep ? How to decode the response into a nice frequency response graph  ?
4) I imagine that there is another better solution?

Thank you.

If you know that the system is linear and time invariant, the easiest method (assuming that you have no noise added in the process) is to let the system act on an impulse function. The Fourier transform of the output is the frequency response of the system.


recovery plane from 4 point imaging using calibration camera I have a camera and its K matrix (calibration matrix) also I have image of plane, I know the real points of the 4 corners and thier correspondence pixel. I know how to compute the H matrix if z=0 (H is homography matrix between Image and the real plane).
And Now I try to get the real point of the plane (3D point) with the rotation matrix and the transltion vector
I follow this paper :Calibrating an Overhead Video Camera by Raul Rojas in section 3 - 3.3.
My code is:

ImagePointsScreen=[16,8,1;505,55,1;505,248,1;44,301,1;];

screenImage=imread( 'screen.jpg');
RealPointsMirror=[0,0,1;9,0,1;9,6,1;0,6,1]; %Mirror
RealPointsScreen=[0,0,1;47.5,0,1;47.5,20,1;0,20,1];%Screen
imagesc(screenImage);
hold on
for i=1:4
    drawBubble(ImagePointsScreen(i,1),ImagePointsScreen(i,2),1,'g',int2str(i),'r')
end

Points3DScreen=Get3DpointSurface(RealPointsScreen,ImagePointsScreen,'Screen');

figure
hold on
plot3(Points3DScreen(:,1),Points3DScreen(:,2),Points3DScreen(:,3));
for i=1:4
    drawBubble(Points3DScreen(i,1),Points3DScreen(i,2),1,'g',int2str(i),'r')
end

function [ Points3D ] = Get3DpointSurface( RealPoints,ImagePoints,name)
M=zeros(8,9);

for i=1:4

M((i*2)-1,1:3)=-RealPoints(i,:);
M((i*2)-1,7:9)=RealPoints(i,:)*ImagePoints(i,1);
M(i*2,4:6)=-RealPoints(i,:);
M(i*2,7:9)=RealPoints(i,:)*ImagePoints(i,2);

end

[U S V] = svd(M);
X = V(:,end);
H(1,:)=X(1:3,1)';
H(2,:)=X(4:6,1)';
H(3,:)=X(7:9,1)';
K=[680.561906875074,0,360.536967117290;0,682.250270165388,249.568615725655;0,0,1;];

newRO=pinv(K)*H;
h1=newRO(1:3,1);
h2=newRO(1:3,2);

scaleFactor=(norm(h1)+norm(h2))/2;
newRO=newRO./scaleFactor;
r1=newRO(1:3,1);
r2=newRO(1:3,2);
r3=cross(r1,r2);
r3=r3/norm(r3);

R=[r1,r2,r3];

RInv=pinv(R);
O=-RInv*newRO(1:3,3);
M=K*[R,-R*O];
for i=1:4
   res=pinv(M)* [ImagePoints(i,1),ImagePoints(i,2),1]';
   res=res';
   res=res*(1/res(1,4));
   Points3D(i,:)=res';

end
Points3D(i+1,:)=Points3D(1,:);  %just add the first point to the end of the array for draw square

end

function drawBubble (xcenter, ycenter, radius, color, string,colorNum)
% draw circle with input xcenter, ycenter, radius, color and text
angles = linspace(0, 2*pi, 40);
bubx = xcenter + (radius * cos(angles));
buby = ycenter + (radius * sin(angles));
fill(bubx, buby, color, 'EdgeColor', color)
text(xcenter, ycenter, string, 'Color', colorNum, 'FontSize', 13, ...
       'FontWeight', 'normal', 'HorizontalAlignment', 'center')

end

My problem is that I don't get good results

1.The point 1 is at (0,0,0) and this is not the real location

2.the points are upside down

What I am doing worng?

I found the answer in the paper: Calibrating an Overhead Video Camera by Raul Rojas in section 3 - 3.3.

for the start:
 H=K^-1*H

1)  Given four points in the image and their known coordinates in the world, the matrix H can be recovered, up to a scaling factor . We know that the first two columns of the rotation matrix R must be the first two columns of the transformation matrix. Let us denote by h1, h2, and h3 the three columns of the matrix H.

2)  Due to the scaling factor we then have that xr1 = h1 and xr2 = h2
Since
|r1| = 1,
then

x = |h1|/|r1| = |h1| and

x = |h2|/|r2| = |h2|.

3) We can thus compute the factor and eliminate it from the recovered matrix H. We just set

H'= H/x

4) In this way we recover the first two columns of the rotation matrix R. The third column of R can be found remembering that any column in a rotation matrix is the cross product of the other two columns (times the appropriate plus or minus sign). In particular

r3 = r1 × r2

5) Therefore, we can recover from H the rotation matrix R. We can also recover the translation vector (the position of the camera in field coordinates). Just remember that

h'3 = −R^t

6) Therefore the position vector of the camera pin-hole t is given by

t = −R^-1 h3


Signal processing for undersampled signal I have got a set of pressure data from a boiler which has been sampled at 45Hz and I need to extract a frequency component of 250Hz from the signal. Using a standard FFT method doesn't display anything apart from frequency components between 0 and 10Hz.

I have read similar questions posted on the forum with no success in solving my problem.

%Matlab code

filename = 'boilerp.xlsx';
columnA = xlsread(filename,'A:A');

Fs =45;

L=length(columnA);

nfft = 2^nextpow2(L);

XT = fft(columnA, nfft)/L;

f= Fs/2*linspace(0,1,nfft/2+1);

plot(f, 2*abs(XT(1:nfft/2+1)))

Many thanks,

Khel

Combining @JRE's comment and @Nikita's answer: if no anti-aliasing filter was applied to the signal, and no other frequencies alias to 20 Hz, then the 250 Hz component you're interested in will appear at 20 Hz in your sampled signal. See https://en.wikipedia.org/wiki/Aliasing#Sampling_sinusoidal_functions for more information on how this works.


About the BPSK signal waveform My basic knowledge is too bad. I am confused about the waveform of BPSK signals. To create a BPSK signal, I generated BPSK symbols with ideal rectangle pulse shaping, which is given by:

Then, the signal after carrier modulation is shown below, which is consistent with the waveform on the textbook.

But, the root raise cosine (RRC) filter is often used in practical. By using the same BPSK symbols above, the waveform after pulse shaping is given by:

After carrier modulation, the waveform becomes like this:

This strange image makes me wonder, is this the correct waveform? Because there is no such waveform in the textbook，I am really not sure about that, Thanks guys. By the way, I have tried that this modulated signal can be demodulated after coherent demodulation and matched filter.

ones and zeros:

Modulated symbols after pulse-shaping:

After carrier modulation:

If you can demodulate then it's probably correct. Perhaps it's just the fact that the carrier is not very fast that it looks weird. Can you increase the symbol period and the carrier frequency. Otherwise, my guess is that you are not convolving correctly. Just double-check if your RC has the correct width (or zero crossings).
For reference I generated this:

Pulse-shaped BPSK with rect() through convolution:

Modulated by multiplying with the carrier tone:

Now with the raised-cosine pulse shape:

Note that by sampling at the correct times (each symbol period), we can recover the BPSK sequence perfectly.

Now modulated:


Convolution Input/output Units How come are the units of convolution for LTI system not consistent?

$$y(t)= \int_{-\infty}^{\infty} x(\tau)h(t-\tau)d\tau$$

$y(t)$ unit is volt, whereas the RHS has volt.time unit!!!!

impulse responses, $h(t)$, of continuous-time systems (those with output that is the same species of animal of the input) have dimensions 1/time.

consider a wire which has impulse response

$$ h(t) = \delta(t) $$

the area of the $\delta(t)$ is equal to the dimensionless $1$.  the width of the nascent delta is in dimension of time, so the height must be in dimension of 1/time.

the impulse response for a simple RC low-pass filter is

$$ h(t) = \tfrac{1}{RC} e^{-^{\, t}/_{RC}} \ u(t) $$

the dimension of $h(t)$ is in 1/time because $RC$ is in units of time.


Guitar String Replucks I'm analyzing guitar string plucks and sustains.  I'm having good success with auto correlation using FFT's.  Now I'd like to detect plucks while the string is still vibrating.  Since I already am periodically taking FFT's to find the pitch, I thought maybe I can take advantage of the fft results and look for changes that might indicate a repluck.  Do I simply add up the bins for some kind of power measurement?

Any ideas?

A pluck might produce significantly more broadband noise than a free string.  An FFT of such noise would show more relative energy outside of all the FFT result bins that are related to (F0 or overtones of) a single pitch.

Also, a free string has a more predictable decay rate in any FFT magnitudes related to the pitch across successive offset FFT frames windowing the sound data.  A sudden change from this decay rate (a stop or a big increase) might infer a disturbance, such as a pluck.


Add harmonics to a FFT decomposed signal i would like to add even harmonics from the data of FFT of a signal.

How do i do that ? do i have to search for higher magnitude frequencies then add the $2f$ then $4f$ etc ?

I understand how to do it from a sin wave for example but not from a noisier signal.

For each added harmonic, for each bin number below N/2, you could divide by the harmonic ratio, then use that resulting number to interpolate a point below in the FFT, unless the division results in an integer, in which case you could just use the value from a lower FFT bin.  Use a high-quality interpolator, such as a windowed Sinc or other polyphase upsampling filter kernel.  For a strictly real result from an IFFT, make sure the values above N/2 in your harmonic frequency synthesis remain conjugate symmetric with those below N/2.


Compare phase and magnitude spectrum results of 2 Images I have read that from Fourier transform we obtain magnitude and phase spectrum. The magnitude spectrum tells you how strong are the harmonics in a image and the phase spectrum tells where this harmonic lies in space.

I have used Matlab to compute the magnitude and phase spectrum of the 2 grayscale images of babies using fftn function, but the all spectrum are difficult to understand.

Can anyone please explain the features given in 1st paragraph by comparing the spectrums of both babies ?

1.phase spectrum

and

2.magnitude spectrum

(If you have results other than in Matlab, explaining with your method is also okay)

image1=imread('D:\baby1.jpg');%Read colour image
image1=rgb2gray(image1);%covert colour image into Gray scale

image2=imread('D:\baby2.jpg');%Read colour image
image2=rgb2gray(image2);%covert colour image into Gray scale

figure,
subplot(1,3,1);
imshow(image1);
title('Gray scale Image of Baby1');

%Apply Fourier transform on Gray scale image
fft1=fftn(image1);
fft1=log(1+fftshift(fft1));%Use log function for scaling

%Find the magnitude spectrum
magnitude1=abs(fft1);
subplot(1,3,2);
imshow(magnitude1,[]);
title('Amplitude Spectrum 1');

%Find the phase spectrum
phase1=angle(fft1);
subplot(1,3,3);
imshow(phase1,[]);
title('Phase Spectrum 1');

figure,
subplot(1,3,1);
imshow(image2);
title('Gray scale Image of Baby2');

%Apply Fourier transform on Gray scale image
fft2=fftn(image2);
fft2=log(1+fftshift(fft2));%Use log function for scaling

%Find the magnitude spectrum
magnitude2=abs(fft2);
subplot(1,3,2);
imshow(magnitude2,[]);
title('Amplitude Spectrum 2');

%Find the phase spectrum
phase2=angle(fft2);
subplot(1,3,3);
imshow(phase2,[]);
title('Phase Spectrum 2');

When you create the transform, an embedded assumption is that you extend your image to infinity by duplicating it vertical and horizontal direction.
Each point in the transform provides information about a 1D wave (frequency that has certain direction based on the point location and phase (how it is moved along the image).
Try to draw a 1D sine wave over a 2D array and then run the 2D FFT to see how the point is created - the higher the frequency the point further from the origin. Starting at a peak will result a $90^0$ phase...


Segmentation of warped text lines I am trying to segment warped text lines found in photos taken from a page of a book (from the camera of a smartphone). This is meant as a preprocessing step before running OCR on them.

I have found two relevant papers on this:

https://sites.google.com/a/iupr.com/bukhari/home/papers/Bukhari-Coupled-Snakelets-IJDAR12.pdf?attredirects=0

https://www.researchgate.net/profile/Rafael_Lins/publication/221000335_An_Efficient_Algorithm_for_Segmenting_Warped_Text-Lines_in_Document_Images/links/0fcfd5063a49cc779c000000.pdf?origin=publication_detail

I read the papers but they're too hard to implement, can anyone help?

Example input no.1:

Example input no.2:

Your question is very broad and I won't go read the papers and implement it for you. If you could show specific problems, we might help you better (or the guys at stackoverflow for programming questions).

Meanwhile, you can have a look at these two resources, which treat a similar problem, including some code:

http://www.pyimagesearch.com/2014/09/01/build-kick-ass-mobile-document-scanner-just-5-minutes/

https://mzucker.github.io/2016/08/15/page-dewarping.html


What is the output of IFFT operations, continuous or discrete time? Assume we have a sequence of symbols
$$S= \{0, 0, 1+j, -1-j, 0 ,0 , 1+j,-1+j\}$$

and that this sequence is a frequency domain sequence and will be input to an IFFT operation.

Is the output a continous time or discrete time signal?

Thanks

The FFT is just a fast algorithm for evaluating the discrete Fourier transform (DFT). The domains for the DFT are discrete in both time and frequency. So, the result of an inverse DFT using the inputs you showed above will be a discrete-time signal.


Power Spectrum Analysis In order to find a peak or max response in the fourier domain as in the spatial domain, I have been studying bartlett's method, welch's method, and the blackman-tukey method for more accurate power spectrum analysis after correlating two functions together in the fourier domain as depicted in this document.

But I had a few questions after reading, which I could not find in any book nor online

1) For the Bartlett method (page 11 of the pdf) what determines the value of k, or is it arbitraty. And following that, does the accuracy of the power spectrum analysis go up if k is higher, and vice versa

2) Same question as 1) but for the welch method (page 12 of the pdf)

3) For the blackman-tukey method (page 14 of the pdf), will i get the same result if I window the signal in the spatial domain, but then do the autocorrelation in the fourier domain?

Hopefully you understand why and in what situations those methods are used.

The integer 'K' in the Bartlett and Welch methods is not arbitrary.  The value of 'K' is a "trade off" that you must make. In both methods, the greater is 'K" the lower is the variance in your final DFT results. That's a good thing. What the PSD's author doesn't tell you, though, is that the greater is 'K" the larger is your final DFT result's frequency-domain sample-to-sample spacing (measured in Hz). And that's usually a bad thing. So you must decide what value of 'K' is acceptable to you. I'm not qualified to comment on the Blackman-Tukey method.


Change duration/tempo of a .Wav file without introducing pops or pitch shift? I am attempting to develop a similar function as Audacity's Effects>Change Tempo effect.

Or more specifically: I am in need of a way to create a function in my own project that can modify the length of a .wav file without changing it's pitch or introducing pops/leakage.

I am writing a standalone application (in C#) that already has some functions to create and save .wav files. But I just realized that I must add some processing to allow me to adjust the duration/length of a .wav so I can sync it to be exactly as long as an associated animation (which can not be altered).

I am very interested in knowing if there are any examples that exist which demonstrate concepts on how to change a .wav file's tempo (without noticeable distortion)?

I should specify that I would not need to adjust the length by much more than ±5%.

What I have tried already is to remove slices of data at regular intervals for a given "slice width", which kind of worked but introduced many pops as the sine waves were no longer smooth/continuous.

I have also tried simply adjusting the playback speed at the player by a multiplier to adjust to a desired duration, but this obviously cause the pitch to change (and somewhat noticeably too).

This is something that I am very anxious to understand as my project now seems to be in trouble unless I can somehow synchronize the audio. Could someone please help me figure out if I can write my own function to process this?

I would suggest the following steps:

identify voices and unvoiced sounds: for this you can use zero-crossing rate and intensity as a basic approach or use more complex methods: cepstrum, harmonicity, others.
the unvoiced sounds can be simply cut, usually without concerns
for the voiced sounds apply PSOLA or WSOLA. For PSOLA there are many implementation and it might be easier. So first you must identify the glotal pulses to be able to identify glotal periods. Then you can "cut" full periods (in a overlap and add sense). Always cut periods in the middle of the phonemes and never near the transition zones (borders)


filter and resample or resample and smooth? Currently doing some signal analysis in python for a major project in my physics degree which is due really soon. I need some help!

Say I have two signals, f(t) and g(t) which are recorded over the same amount of time.

f(t) and g(t) both feature different amounts of noise.

f(t) and g(t) are recorded at different sampling rates so I have two different quantities of measurements.

For analysis I imagine I will need:

a) Interpolation (or resampling?) both signals to get the same sampling rate.

b) Smoothing/filtering to reduce noise.

a then b or b then a?

I would like to ideally have two signals with equal sampling frequencies. I would also like all values of f(t) and g(t) to occupy the same points on the t (x) axis. In addition to this I would like to smooth the data as much as possible with a minimal loss of information.

Could you describe what you want to find in your data -- peaks, trends, model ?
and what kind of noise there is -- $\infty$ possibilities ?
A plot would be good.
Without some more info, I can only say, with others, "try it and see".

As you may know,
a good way to smooth and resample (interpolate) in Python is with
scipy.ndimage.map_coordinates,
along the lines

import numpy as np
from scipy.ndimage import map_coordinates

def smooth_resample( y, newlen ):
    """ resample y to newlen, with B-spline smoothing """
    n = len(y)
    newgrid = np.linspace( 0, n - 1, newlen )
        # e.g. n = 3, newlen = 5 -> [0 .5 1 1.5 2]
    return map_coordinates( y, newgrid, mode="nearest", prefilter=False )

map_coordinates also works for 2d images, 3d voxels ... hence the name "ndimage".
For more, stackoverflow has quite a few Q+A s that use it.

(Cubic-spline interpolation uses the nearest 4 points to each x in newgrid,
whereas dsp smoothing filters are usually much longer.
For smoothing that preserves peaks,
see Savitzky-Golay filters, e.g.
Numerical Recipes p. 771 ff.)


What Are Less Computationally Demanding Alternatives to the Viterbi Decoder? What are less computationally demanding alternatives to the Viterbi Decoder?

Ideally what I would like is a list of the most commonly used approximate methods, along with brief pros and cons.

An straightforward alternative to a Viterbi decoder for decoding convolutional codes is a Sequential Decoder.   The algorithm is simpler, but the performance isn't as good.   For example, the common constraint lengths for systems using Viterbi decoding are typically k = 7 or k = 9.  A system decoded with a sequential decoder typically has a constraint length in the thirties or so.   The complexity of a Viterbi decoder would be prohibitive at such a long constraint length, and the performance of a Sequential Decoder is very poor with the short constraint lengths used with a Viterbi decoder.

Another alternative would be to use an APP or MAP decoder, but those are more complex than a Viterbi decoder (since they're essentially two Viterbi decoders running together in opposite directions).  But if you just want to say you're not using a "Viterbi" decoder that might qualify.


Applying multiple IIR and FIR filters at once I would like to appy multiple filters to a signal at once.
Let's assume I got two frequency bands I want to extract. I compute a Butterworth Filter for 300-500Hz and one for 700-800Hz. I would like to combine those filters into a single filter and apply it to the signal. The resulting signal should contain only frequencies in range 300-500Hz and 700-800Hz.

I assume that if the filters are FIR (finite impulse response) filters I just need to add them. Using that single filter, I apply a convolution to the signal. Is that correct?

I am pretty lost if it gets to IIR (Infinite impulse response) filters which are - as I understand - divisions of polynoms in a radial space I do not understand. How can I combine two IIR filters?

Kind regards,
Thomy800

When someone says “I want to do x” , it helps to know why?
Instead of going into why you want to do your “x”, it might be instructive to go into why it is probably not a good idea.

There are some filters where there are  mathematically equivalent recursive and non recursive forms, but they don’t have the same robustness to rounding errors.
Combining recursive and non recursive forms into a single filter is possible. The numerator in a ratio of polynomials that typically make up a transfer function can be thought of as a non recursive component , so one can end up with a combined filter that is recursive.  You can also truncate the impulse response to end up with an overall non recursive filter. In either case it really can’t be both.
Polynomials are delicate. Small perturbations in coefficients can make big changes in roots, particularly when those roots are poles near the unit circle. Analytically, it can be stable but in practice it isn’t. In most cases, we factor transfer function is structures like biquads. Lots of cascades little filters are almost always better behaved.
Humans delight in their cleverness, but very few things are universally acknowledged as being cleaver. In practice, if this filter is “product” someone else will come along and maintain your code. In most cases, that person ( which might be you after a few years) will not share your joy. They are more likely to curse you and if you are very cleaver, they will curse your mother and father. One should chose clarity over cleverness.
You really need to do something a number of different ways in this sort of framework to make tradeoffs. The term “exhaustive enumeration” applies to this
problem. The different ways a band pass can be implanted is really not hard to compare most of the
salient parameters. You have to appreciate how much harder you make
your filter optimization
When a parent gives their child a bicycle, it will take a number of years before they enter X-Games. Given that you asked the question, the ultimate answer, is yes but will require great skill.

.

“Whatever doesn’t kill you makes you stronger” I encourage you to follow up on your ideas. Try to make it work. Endeavor to make me wrong,


A priori spectral efficiency of OFDM scheme Of course I can compute the spectral efficiency as the ratio of the measured bit-rate (of data bits, no header, no preable, nor redundancy bits included) and the measured spectral bandwidth. This  would be an a posteriori approach.

However, I wonder if there exists an analytical expression for a priori determination of the spectral efficiency of an OFDM communication system.

The system uses preable, header, forward error correction (in and out), cyclic redundancy check, cyclic prefix and taper window. And the modulation by subcarrier can also be changed, as well as the number of subcarriers, and which of them are nulled or pilots.

Is there some kind of formula that includes all of these parameters and returns the effective bit rate and spectral efficiency?

Taking a look to the source code of the liquid-dsp project, and to the ofdmflexframegen.c file in particular, I found a partial answer to my question.

The spectral efficiency is:
$$\beta = \frac{R_b}{W}~~~~~~~~~~~~~~~~~~~~~~~~~~~(1)$$
where $R_b$ is the bit rate and $W$ is the OFDM bandwidth.

The bit rate is:
$$R_b = \frac{B}{T}~~~~~~~~~~~~~~~~~~~~~~~~~~~(2)$$
where $B$ is the data bits (of the decoded payload) and $T$ is the total transmission time.

Then, putting (2) into (1):
$$
\beta = \frac{B}{TW}~~~~~~~~~~~~~~~~~~~~~~~~~~~(3)
$$
and the product $TW = N$ is the number of samples of the transmission.

As for $N$, we have:

$$
N = (M + c)(N_{preamble} + N_{header} + N_{payload})~~~~~~~~~~~~~~~~~~~~~~~~~~~(4)
$$
where $M$ the total number of subcarriers, $c$ the length of the cyclic prefix.
$N_{preamble}$, $N_{header}$ and $N_{payload}$ are the number of OFDM symbols of the preamble, header and payload, respectively. Each OFDM symbol is comprised of $M + c$ samples.

Finally, the spectral efficiency for such a system would be:
$$
\beta = \frac{B}{(M + c)(N_{preamble} + N_{header} + N_{payload})}~~~~~~~~~~~~~~~~~~~~~~~~~~~(5)
$$

However, practical measurements of the spectral efficiency of my communication system are a little higher of what this formula predicts. It is like if the measured bandwidth was actually smaller than the predicted one. The difference increases with the payload size. I would appreciate any extra idea in this regard.

EDIT 1: Taking a deeper look to my system, it is possible that the difference between the spectral efficiency predicted by Eq. 5 and what I am measuring is due to the presence of null subcarriers. Null subcarriers at both sides of the OFDM spectrum are not measured, but they are included in Eq. 5. I wonder if I may just substract the number of null subcarrier from the term $M+c$ in Eq. 5. Does that make any sense?


Possible to break out complex WiFi signal? Is it possible to break out the i and q of a complex WiFi signal?  If I have a laptop with a standard WiFi card in it, could I break out the i and q in analog, or upon sampling in an FPGA via digital means?

It must be possible since it is being done on a WiFi receiver, right?

The in-phase and quadrature components are easy to break out. You just need to demodulate using in-phase and quadrature local oscillators. See the diagram below.

The received signal $r(t,k)$ is mixed with a $\cos$ and a $-\sin$ wave of the appropriate frequency, and then low-pass filtered.

The actual IEEE 802.11 modulation schemes used are a little more complicated than indicated, but not overly so.


correcting mirror ghosting effect I have an application in which I have to process images taken through a mirror. The camera and mirror positions are fixed, and the object is always in the same position with respect to the mirror. The problem is that I get a ghosting effect because the image is also reflected by the transparent surface of the mirror.

I will try to get a first surface mirror for my application, but I wonder if this ghosting effect can be corrected by software. Any hint would be appreciated (a paper, algorithm, anything). I believe this should be possible, since the camera and mirror are both fixed, and the object is always at the same distance from the mirror, but I don't know where to start looking.

Thank you for any hint/advice you can provide.

We can assume
$g=i(x,y)*f(x,y)+n(x,y)$, where $g$ is the above images, $i$ is the original image , $f$ is a degradation function and $n$ is image noise. If you find $f$ and deconvolve it with $g$, you'll get $i$ plus some noise. Note, simple deconvolution would not work, because the noise ($n$) in the process is amplified. The Wiener filter can used to restore images that have been degraded by a known degradation function.Here I am not sure what is degradation function.
I suggest you read the following material as well:
Lecture notes of Prof. William Hoff


FFT padding to power of 2 I have a vector v1 of size N1, I want to pad it in a certain way until it reaches a size power of 2 (v2 of size N2). I need that the first N1 coefficients of the FFT output to stay the same as when I would apply FFT(v1), is that possible? If so, how to?

It is possible but it is usually less complex to do the FFT of size n1.
To be honest given the FFT of a signal of size N1 padded to N2, the most efficient way that comes to my mind to find the coefficients of an FFT of the same signal without padding is to do the inverse FFT of the signal, removing the padding and do the FFT of size N1 with those clever libraries such as fftw.
Why?
Let $p(f) = \sum_{k=0}^{N_1-1} x_i \exp(2j \pi f k)$, the FFT of size $N_1$ is $p(n_1/N_1)$ but the FFT of size $N_2$ (with the padded signal) will give you $p(n_2/N_2)$.
So let's see for which $n_1$ we find solutions with integer $n_2$
$$\frac{n_1}{N_1} = \frac{n_2}{N_2}$$
it turns out that $n_2 = n_1 N_2 / N_1$, let $N_2 = 2^{a_2}$ and $N_1 = 2^{a_1} b$ for an odd $b$. Then we have $n_2 = n_1 2^{a_2 - a_1} / b$, i.e. only the coefficients at indices multiples of $b$.
In a more general way we could say that only the coefficients with indices multiples of $N_1 / \gcd(N_1, N_2)$ can be recovered.


Structuring Kalman filter for tracking problem where only position is known I'm new to Kalman filters and my extensive web search about them has helped me understand the majority of it (or so I think). However I still need some light shed on my problem formulation.

I have a set of points $(x,y)$ provided by some other part of my program that correspond to the coordinates of a target whose position is being estimated. Let's assume the coordinates have have error with known variances $\sigma_x$ and $\sigma_y$, zero mean and are normally distributed. I have no other information (for example, regarding velocity or acceleration). The path of the target is completely random as it modeled by a walking person and it can remain long periods of time just standing still. The samples $(x,y)$ are obtained with a well defined and constant time interval $\Delta t$.

Is a Kalman filter the best approach to smooth the path in my case? It seems to me that, although the person can walk freely and apparently randomly, their movements still have a high correlation in a small scale, so a Kalman filter could be a choice.
There is an example in Wikipedia on Kalman filtering with a truck whose velocity is not known and where its unknown acceleration is modeled by the noise vector $\boldsymbol{w_k}$. Is this the best approach to solve my problem?
Also, how would the system be formulated in the case of two dimensions $(x,y)$?
Is there any harm in feeding the same coordinates $(x,y)$ to the filter for long periods of time? (Meaning the person is standing still)
How does the Kalman filter differ from, for example, the Savitzy-Golay approach to data smoothing?

Thank you in advance!

Just a slight nit $\sigma^2$ is a variance $\sigma$ is a standard deviation.

Let me answer your last question first.  In a KF, one has measurements and states, and they are not usually the same.  In a Savasky Golay filter, one is estimating the derivatives of the measurement, not the state.  In your case the measurements are trivially related to the states.

I'm guessing but I think you want something like
$$
\left| \begin{array}{c} x \\ y \\v_x \\ v_y \end{array}\right| =
\left| \begin{array}{cc}
1 & 0 & \Delta_t & 0 \\
0 & 1 & 0 &\Delta_t \\
0 &  0 & 0 & 0 \\
0 & 0 & 0 & 0 \end{array} \right| +
\left| \begin{array}{c} 0 \\ 0 \\ w_1 \\ w_2 \end{array} \right|
$$
for the state model  $\Delta_t$ is your time step and $w_1, w_2$ iid Normal noise.

If this model is your problem,  yes the KF is what is wanted.

What you will need to figure out, is what is process noise and what is measurement noise.

A KF needs measurement noise to work, process noise is optional, but you would need a different state model if you want to just use measurement noise.


Particular Solution to Difference Equation I have a system given by
$$y[n] - \frac{1}{4} y[n-1] - \frac{1}{8} y[n-2] =3x[n] $$
I want to solve for $y[n]$ for $x[n]=(\frac{1}{2})^nu[n]$.
The complementary solution evaluates to $[k_1(\frac{1}{2})^n+k_2 (\frac{-1}{4})^n]u[n]$.
But when I evaluate the particular solution consider $y_p[n]=k(\frac{1}{2})^n u[n]$ I get an absurd answer as $0=3$. And boundary conditions are not given either. What initial conditions am I to consider? I want to solve this difference equation without using Z-transform or Fourier transform.
Edit:
I took the particular solution of the form $$y_p[n]=kn(\frac{1}{2})^nu[n]$$ since the input is of the same form as the roots of the characteristic equation.
and get k as equal to 2.
and evaluate the coefficients of complementary solution as
$$y[n]=\frac{8}{3} (\frac{1}{2})^n u[n]+\frac{1}{3} (\frac{-1}{4})^n u[n]+2n (\frac{1}{2})^n u[n]$$
but the solution given is $$\frac{1}{3} (\frac{1}{4})^n u[n] + 4(n+1)(\frac{1}{2})^{(n+1)} u[n+1] + \frac{2}{3} (\frac{1}{2})^n u[n]$$
I don't where the anomaly in my evaluation is? Or is the given solution wrong?
I'm new to the concept of difference equations.
I just want to know if my understanding is right.
I wish to solve this difference equation without using z transform.

The solution you came up with is the correct homogeneous solution (i.e., when $x[n]=0$). Thing is, this is a non-homogeneous difference equation, and its solutions are of the form
$$y[n]=y_h[n]+y_p[n]=k_1\Big(\frac{1}{2}\Big)^nu[n]+ k_2\Big(\frac{-1}{4}\Big)^nu[n]+y_p[n]$$
where $y_p[n]$ is the particular solution that solves for the $x[n]$ term, and you are missing the corresponding particular solution in your evaluation. One way to find $y_p[n]$ is to guess the form of $y_p[n]$ with unknown coefficients given the form of $x[n]$, and use the method of undetermined coefficients to find those coefficients.
For $x[n]=r^nu[n]$, you would try $y_p[n]=ar^nu[n]$ at first. But, in this case $x[n]$ corresponds to a mode of the system. Therefore, the particular solution cannot be of the form $a(\tfrac{1}{2})^nu[n]$, but it has to be of the form $an(\tfrac{1}{2})^nu[n]$. Therefore, the solution of the equation is
$$y[n]=k_1\Big(\frac{1}{2}\Big)^nu[n]+ k_2\Big(\frac{-1}{4}\Big)^nu[n]+an\Big(\frac{1}{2}\Big)^nu[n],$$
and you find $a$ by evaluating this solution in the equation.


Causality, Linearity, and Time Invariance for Systems Described by Linear Constant Coefficient Differential Equations I am currently using Signals and Systems by Alan Oppenheim as a reference to learn about LTI systems.

Before introducing systems represented by linear constant coefficient differential equations, the author first considers general systems of the form $x(t) \rightarrow y(t)$. When discussing linearity of systems, he introduces the "zero in-zero out" property of linear systems (not necessarily time invariant or causal), that is, given a linear system, if the input $x(t)$ is zero for all $t$, then the output $y(t)$ must be zero for all $t$. In addition, a linear system is causal if and only if $y(t) = 0$ for all $t < t_0$ if $x(t) = 0$ for all $t < t_0$, a condition known as initial rest, from which I concluded that initial rest for linear systems implies causality but not time invariance.

Now given a linear constant coefficient differential equation $$Ly(t) = x(t),$$ where $L$ is a linear differential operator with constant coefficients, the author emphasizes that although the differential equation is linear, the system that it represents may or may not be linear depending on the initial conditions. I find this confusing because given two input output pairs $x_1(t) \rightarrow y_1(t)$ and $x_2(t) \rightarrow y_2(t)$, we have $L(ay_1(t) + by_2(t)) = ax_1(t) + bx_2(t)$ for any constants $a, b \in \mathbb C$ by virtue of the linearity of the operator $L$. It seems that $ax_1(t) + bx_2(t) \rightarrow ay_1(t) + by_2(t)$ regardless of the initial conditions, why is the system not necessarily linear in that case? What am I missing here?

In addition, the author states that for the system to be linear, the initial conditions have to be zero, and for it to be causal and time invariant in addition to being linear, its initial conditions must obey initial rest. But isn't initial rest in general a condition required for causality, not time invariance? Besides, doesn't this imply that a differential equation can only represent a causal system if the input is zero before some time? Given an input that is nonzero for all time, how do we impose causality via initial conditions that correspond to initial rest? I hope my question makes sense. Thank you in advance.

Note that the system is not uniquely characterized by the linear differential operator. You need auxiliary conditions, and dependent on these conditions, the system may be linear, causal and time-invariant or not. So if you just use the operator $L$ to show that the system is linear then you're actually not showing anything because $L$ is a linear operator by definition.

The output of a system described by a linear constant coefficient differential equation can be split up into two contributions: the zero-state response (ZSR) and the zero-input response (ZIR). The ZSR is the response of the system with zero initial conditions, and, consequently, the ZSR is fully determined by the input signal. The ZIR response is the response of the system caused by non-zero initial conditions, and it is independent of the input signal. Since non-zero initial conditions cause an output term that does not depend on the input signal, scaling or shifting the input signal will not affect that term. Consequently, a system with non-zero initial conditions can neither be linear nor time-invariant, because if it were, the output should scale and shift according to the scaling and shifting of the input signal.

Causality simply means that the system's response at a given time does not depend on the future of the input signal. Clearly, causality is a property of the system, so it cannot depend on the input signal. A causal linear time-invariant (LTI) system is characterized by an impulse response $h(t)$ that satisfies $h(t)=0$ for $t<0$, and it is irrelevant whether the input signal is zero before a certain time or not.

Also take a look at this related answer.


How to create a distortion compensation filter I am creating a channelizer that consists of a complex mixer, CIC decimator, and a FIR compensation/decimation filter.  The final FIR filter may, if it matters, be implemented as more than one filter.

My question is, how do I design a filter such that it compensates for the very non-flat frequency response of the CIC filter?  Do you create the frequency response that you want by calculating the reciprocal of the CIC's response and then run it through an inverse FFT to get the impulse response?

As you can see, though my particular problem pertains to CIC filters, my question is really about how you create any kind of distortion compensation filter.

Thanks for your time.

The following addresses the CIC compensation and is not a general "distortion" technique.  But it is a straight-forward method to "compensate" non-varying "distortion".  If the frequency response is known the inverse of the frequency can be used to compensate.  Examples like the CIC filter, where a poor filter might be used because of reduced complexity, are compensated later the signal chain.  In this example the frequency response is know and the inverse can be used.  Note, with multi-rate filters you only want to use the "usable" spectrum after decimation.

Generally, to compensate a CIC filter the inverse of the CIC filters response
can be used to generate the compensation filter.  The CIC has
a response of (see reference [r2] or [r3])

$$
  H(\omega) = \left| \frac{sin(\omega D/2)}{sin(\omega M/2)} \right|^N
$$

Where D is the differentiate delay, M is the decimation rate, and N is the
filter order (number of cascaded filters).  The inverse can be specified as

$$
  H(\omega) = \left| \frac{sin(\omega M/2)}{sin(\omega D/2)} \right|^N
$$

Once we have the
frequency response of the compensation filter, we can simply choose the length
of FIR filter that we desire.  The length of the FIR is application
specific.  Obviously the longer the FIR filter the better compensation.

The following are plots of this straight forward compensation.

The following is the Python code to create the frequency responses and plots.

import numpy as np
from numpy import sin, abs, pi
import pylab

D = 1; M = 7; N = 3

Hfunc = lambda w : abs( (sin((w*M)/2)) / (sin((w*D)/2.)) )**N
HfuncC = lambda w : abs( (sin((w*D)/2.)) / (sin((w*M)/2.)) )**N

w = np.arange(1024) * pi/1024

G = (M*D)**N
H = np.array(map(Hfunc, w))
Hc = np.array(map(HfuncC, w))
# only use the inverse (compensation) roughly to the first null.
Hc[int(1024*pi/M/2):] = 1e-8
plot(w, 20*log10(H/G))
plot(w, 20*log10(Hc*G))
grid('on')

See [r1] for other approaches and $sinc^{-1}$ approximation.

[r1] Altera, "Understanding CIC compensation filters"

[r2] R. Lyons, "Understanding Digital Signal Processing",
    2nd ed., Prentice Hall, Upper Saddle River, New Jersey, 2004

[r3] R. Lyons, "Understanding Cascaded Integrator Comb Filters"


Derivation of Circular Mean Square Error I would like to understand how Eq. (36) in [2] was derived:

The rationale behind the definition of circular sample mean in Eq. (37) is clear, but there is no motiviation for the CMSE definition in Eq. (36).

[2] Lovell, Brian C., and Robert C. Williamson. "The statistical performance of some instantaneous frequency estimators." IEEE Transactions on Signal Processing 40.7 (1992): 1708-1723.

The first term in $e^2_p$ is just the variation of the estimates around the true value.

The second term is due to the problem that happens at the end of the periodic region. Sometimes, the noise is enough to move the estimate from $-\pi+\alpha$ to $\pi - \beta$.

For example, the orange plot in the figure below represents the variation of the estimates for a true value around zero.  This variation is captured by the first term.

The blue plot in the figure below represents the variation of the estimates for the true value around the end of the period.  The trouble is that some estimates end up all the way at the other end of the periodic region.  This causes a bias to appear.  The second term in (36) is aimed at capturing that bias.


If noise is your signal, what is your noise? Consider the following contrived situation. Imagine a Gaussian white noise process $x[t]$, with bandwidth $Δf$, with PSD equal to some quantity $A$ which you would like to measure.

So the way to measure this seems to measure the variance of the process $x[t]$, which by Parseval's theorem will be $AΔf$.

So you measure points with some frequency $f_m$, probably $2Δf$. And at a lower rate, say $f_v$, you compute the variance of the preceding block of $f_m/f_v$ points and take that to be a measurement. What will the noise/variance be in this measurement of the variance of $x[t]$? How can I approach this question?

There is a derivation over on the math.SE site that might be useful.

Provided your data is Gaussian with variance $\sigma^2$.

If you use the unbiased estimator of the population variance:
$$
S^2_{N-1} = \frac{1}{N-1} \sum_{n=0}^{N-1} (x[n]- \bar{x})^2
$$
where $N$ is the number of points in your batch and $\bar{x}$ is the sample mean, then the variance of this estimate will be:
$$
\frac{2 \sigma^4}{N-1}
$$


Sampling: How should I sketch a signal spectrum? As revision, I am going through Signals and Systems, and in Section 7.1the authors are sketching representative spectra for x(t) and xp(t)as follows

However when I search for "How to draw a signal spectrum" I get results (here and here), where signal spectra are shown like this:

For example,how should I go about sketching the the spectrum of signal x(t), where x(t) = 4cos(8000πt)?

Here are my two answers:

Which one is the correct approach? Any suggestions would be appreciated.

Figure 7.4 just concerns itself with an example signal, which happens to have triangular spectrum.

your $x(t)=a\cos(2\pi f t)$ signal is just that: a single cosine tone, and thus has the "shape" of two dirac impulses with integral area $\frac a2$ each.


By using array antenna and a MUSIC algorithm, how many targets could be detected? By using an array antenna which has N receivers, How many targets could be detected and classified when we use the MUSIC(MUltiple SIgnal Classification) algorithm?

not a foolish question at all.

for a uniform linear array.

there are a number of different answers in the literature that do not agree.

the first issue is that direction finding and beamforming are different and have different upper bounds on the number of signals that can be resolved.

for beamforming $N-1$ sources under the narrow band assumption is typically the upper bound. you might assume polarization might double that number.

again under the narrow band assumption, the rank of the cross spectrum matrix is an upper bound for direction finding.
One might reasonably assume $N$ for $N$ sensors but one can construct a virtual array from the negative coarray.

The virtual array is based subtracting outputs of pairs of real sensors.  One can form $N^2$ such pairs if $x_k-x_j$ is distinct from $x_j-x_k$.  If one were to believe this, the upper bound at a single frequency is $N^2$ sources.

The $N^2$ result is given in the peer reviewed literature by an eminent IEEE Fellow but I can’t say it is compatible with my intuition.


Is $O (N \log N)$ FFT speed the fastest we can ever attain? I am wondering about whether or not there is a theoretical limit as to the speed at which we can compute a DFT. We all know that the FFT executes in $O (N \log N)$ time. However, is this a lower bound of some sort?

Might there be yet faster algorithms that have yet to be discovered which we do not currently know about that are faster, or is this a holy boundary that cannot be crossed?

Edit: Apparently this is a known unsolved problem in CS, (thanks to @John).

So, what are the details of the problem? Where would one start to try to make improvements?

While the general problem remains an open problem in theoretical computer science, there are (randomized) algorithms known as sparse fast fourier transforms which can approximate DFTs to a certain approximation factor when there are only a few DFT coefficients you are concerned with.

You can find some details in these papers and related papers (downloadable here):

Hassanieh, Haitham, et al. "Nearly optimal sparse Fourier transform." Proceedings of the 44th symposium on Theory of Computing. ACM, 2012.
Hassanieh, Haitham, et al. "Simple and practical algorithm for sparse Fourier transform." Proceedings of the Twenty-Third Annual ACM-SIAM Symposium on Discrete Algorithms. SIAM, 2012.


Fast/Slow Fading: How are doppler shift and coherence time related? In the literature I've read that fast fading happens if the bit duration is larger than the coherence time of the channel.
The coherence time of a channel is proportional to $\frac{1}{f_{Doppler}}$. So coherence time and doppler spread are related somehow. Why do frequency shifts determine the time span for which the channel coefficients are constant, e.g. is there any physical meaning of this relation? Is there any intuitive explanation?

Here's an intuitive explanation.
Imagine a communication system (wireless) in a static environment (maybe picture the inside of a closed warehouse). There will be an amount of fading, but it is static: it never changes, because the channel doesn't change.
For the channel to change, the multipath scenario must change, and that means that objects must move around (it can be the reflectors/scatterers, the antennas, or both). If they move relatively slowly, you have "slow" fading; if they move relatively quickly, the fading becomes "fast".
Now, whenever you mix waves and moving objects, you get Doppler shifts. Furthermore, the faster the movement, the higher the Doppler. In consequence, you can measure how fast the channel is changing by either measuring how fast the objects are moving, or by how large the Doppler is.
Since the Doppler shift is easy to quantify, it is the preferred way to measure how fast a fading channel is changing.


How does the world's simplest sigma-delta DAC work? I have run into a simple delta-sigma (supposedly, where is the delta part?) DAC implementation utilizing an accumulator and then using the overflow bit as a modulated 1-bit output.
From https://www.fpga4fun.com/PWM_DAC_2.html:
module PWM(clk, PWM_in, PWM_out);
input clk;
input [7:0] PWM_in;
output PWM_out;

reg [8:0] PWM_accumulator;
always @(posedge clk) PWM_accumulator <= PWM_accumulator[7:0] + PWM_in;

assign PWM_out = PWM_accumulator[8];
endmodule

I have implemented this and tested on an FPGA and compared it to a delta-sigma DAC from Xilinx app note: https://www.xilinx.com/support/documentation/application_notes/xapp154.pdf
module ds_modulator (
        output DACout,
        input [7:0] DACin,
        input Clk,
        input Resetn
);

reg DACout;
reg [9:0] DeltaAdder;
reg [9:0] SigmaAdder;
reg [9:0] SigmaLatch;
reg [9:0] DeltaB;

assign DeltaB = {SigmaLatch[9], SigmaLatch[9], 8'b0};
assign DeltaAdder = DACin + DeltaB;
assign SigmaAdder = DeltaAdder + SigmaLatch;

always @(posedge Clk, negedge Resetn) begin
    if (!Resetn) begin
            SigmaLatch <= 10'b1111111111;
            DACout <= 1'b0;
    end else begin
            SigmaLatch <= SigmaAdder;
            DACout <= SigmaLatch[9];
    end
end

endmodule

I am looking for as simple as possible delta-sigma DAC containing minimum logic while still maintaining reasonable performance. The very simple accumulator-overflow solution seems to work great. I understand it conceptually--high values will cause more frequent overflows and low values will cause less frequent overflows. Are there any disadvantages? I am using 120 MHz sampling frequency and generating sine waves between 100 kHz and 500 kHz through direct digital synthesis.
I am interested in this simple approach. Can someone shine light on why it works so well while being so simple? Also, if there is an official name for this approach or if there is some analysis somewhere online, it would be greatly appreciated if you can reference it in your reply.
I also found the simple accumulator-overflow DAC reference here (page 71): https://github.com/hamsternz/IntroToSpartanFPGABook/blob/master/IntroToSpartanFPGABook.pdf
Thank you.

This is just a first-order sigma delta DAC with the input as an unsigned binary number. The classical variant with a signed digital number is shown below where the MSB as the sign bit would subtract (Delta) from the input. In this signed variant, consider the input being a scaled digital number ranging from +/- Full Scale, and the feedback would either be +Full Scale or -Full Scale; in this case the accumulator can count up or down with a ramp rate based on the difference between the actual input value and +/-Full Scale with the output toggling when the accumulator output crosses midrange.
In comparison for the OP's case, the input is unsigned and the accumulator grows only in the positive direction with a ramp rate based on the difference between the input and +Full Scale. When the accumulator overflows, it is reset to 0 plus the overflow, which is essentially subtracting Full Scale (which is the Delta in this case).

This has all the drawbacks of any first order Sigma Delta: lower noise supression than a higher order architecture (-9 dB/octave) requiring a higher oversampling ratio to achieve a desired equivalent number of bits, and suffers from pattern noise leading to higher spurious content if the input is not changing.


How to smooth a signal? The easiest way to smooth a signal is by moving window average.

A more advanced way is to use a Savitzky-Golay filter.
From wikipedia:

  The main advantage of this approach is that it tends to preserve
  features of the distribution such as relative maxima, minima and
  width, which are usually 'flattened' by other adjacent averaging
  techniques (like moving averages, for example).

There is also a whole range of window functions.
As I understand this: any finite filter wil cause spectral leakage,
but the moving window average is the worst. E.g. a Gaussian window is better in this respect.

Is flattening and spectral leakage the same?

When should one use Savitzky-Golay and when should one use Gauss, Hann, Hamming etc?

Thanks in advance for any answers!

When using a moving average filter, or rectangular convolution in the time domain, the time domain smoothing and so-called spectral leakage are nearly opposing effects.  The smoothing is due to removing high frequency content, and the leakage is the name given for portions of this frequency content that end up not being removed.  Since sharp peaks and transients can be partially composed of high frequency content, smoothing by this kind of convolution can diminish them.

Time-domain convolution windows other than rectangular may remove less of the very nearby frequencies, but better remove more of the far away frequency content (relative to some transition frequency, which is related approximately to the window width).

An adaptive polynomial regression will remove different spectral characteristics from a signal at different points in time, which could either expose or hide "important" elements of the signal, depending on the signal, and what you are looking for.


Sine wave amplitude modulation - Analysis using Audacity I am trying amplitude modulation in C++ on a sine wave which is a 1000 Hz sine wave sampled at 48000 Hz. The original sine wave is like this (In Audacity) -

When I basically double the amplitude for each sample, it works except for the flipping up of the extremities -

When I multiple each sample with 4, this is what I get -

Why is this the case, when the equations that I am using in the algorithm don't have anything of this sort? What if I want to remove this 'error'?

P. S. When I import the .raw data, I am using Signed 32-bit PCM encoding, Big-endian Byte order and single channel, with zero offset and 48kHz sample rate.

Your samples have a finite bitwidth. For example, assume your samples have 16 bit width; if you increase their amplitude beyond what these 16 bit can represent, you end up with an overflow.

In case your audio samples being signed integers, that overflow typically first manifests as a sign inversion – exactly what you're seeing here. (That's because only the lower 15 bits in your 16 bit signed integer are actually the number – the uppest bit is the sign; you shift a 1 there, the result becomes a negative number in 2-complement represenation.)

This isn't a limitation of audio processing only – in fact, integer limits is something basically (nearly) every programmer needs to be (slightly) aware of. Your number variables can't hold arbitrarily large values!


Unit impulse response of a cascade interconnection of three discrete-time systems I am nearly at the end of finishing a problem in my textbook but I couldn't understand something in the answer;

I did everything to the point I found the overall response of the system in terms of $h_{1}[n]$ which is $$h[n]=h_{1}[n]+2h_{1}[n-1]+h_{1}[n-2]$$ Now while comparing each $n$ value with the given graph for overall response, the answer is saying that for $n=0$;  $$h[0]=h_{1}[0]$$ but I didn't get why it isn't like this: $$h[0]=h_{1}[0]+2h_{1}[-1]+h_{1}[-2]$$ How do we know $h_{1}[n]$ is zero for $n$ values smaller than zero?

HINT:

Consider the cascade interconnection of three causal LTI systems, [...]"

What does that imply for the impulse response $h_1[n]$?
What is more, since the total impulse response is causal, and you know from its definition that $h_2[n]$ is causal, $h_1[n]$ must also be causal.


How to detect a spiky signal embedded in a periodic signal Say I've got a signal that I know is composed of periodic components. However, I'm worried that there are "spiky" events that are embedded in this signal, and potentially occurring periodically. I can't just look at the power spectrum of the signal, because the spiky event will show up across many frequencies in the power spectrum.

What kinds of methods are out there to tease apart components of a signal that can reasonably be decomposed into sine waves, from those that cannot (like a spike)?

You can use linear prediction: look at the residual of an AR model. If the signal is a sum of sinusoidal components - and provided the order of the model is appropriately chosen - the prediction error will be small. Unexpected spikes will directly show up in the error signal.


Apply Singular Value Decomposition on Gabor filters According to my question in this link, I applied (with a big help from lennon310) the separability of Gabor filters of any orientation.
The only disadvantage of this method is the convolution in the complex domain. So we can treat  the Gabor ﬁlter G(x,y) as a two dimensional matrix that has a low rank, since the ﬁlter is fairly smooth at lower frequencies. Singular value decomposition can then be used in order to decompose the ﬁlter into a linear sum of real separable ﬁlters.

Here, ui, vi are the columns of the orthogonal matrices U,V. s_i is the singular value. The convolution of the ﬁlter G(x, y) with the image I(x, y) can now be approximated using:

Please I need a help about how can I edit my code in this link in order to apply that in MATLAB. Any help will be very appreciated.

I won't recommend you apply SVD in your Gabor filter, since it does not bring much benefit but increase the computing load. If you implement SVD, you may not separate the filter at first, instead the SVD is performed on the 2D Gabor filter:

                for i = -filtSizeL:filtSizeR
                    for j = -filtSizeL:filtSizeR

                        if ( sqrt(i^2+j^2)>filtSize/2 )
                            E = 0;
                        else
                            x = i*cos(theta) - j*sin(theta);
                            y = i*sin(theta) + j*cos(theta);
                            E = exp(-(x^2+G^2*y^2)/(2*sigmaq))*cos(2*pi*x/lambda(k));
                        end
                        f(j+center,i+center) = E;
                    end
                end

               %% SVD %%%%
               [u,s,v]=SVD(f);

With u, s, and v, you implement the convolution with your image:

    convv = zeros(size(image_double));
    for i = 1:filtSizeR+filtSizeL+1
        convv1=imfilter(image_double*s(i,i),  u(i,:),'conv');
        convv2=imfilter(double(convv1),v(:,i)','conv');
        convv = convv + convv2;
    end

    figure
    imagesc(imag(convv));

Yet I don't think this method is as good as your separable f and g method.