Best Approach on How To Sample and Process I'm looking for some advice, I am hoping one of the many forum experts could help me with some advice on the sampling and processing of signals. I am using an STM32F4 processor with 12 bit ADC and DSP core. Presently, I can use either CMSIS DSP or Siglib libraries. Here's a point description of the problem space:

I can sample at whatever rate is required: in front of the ADC is a 10th order analog amplifier and LPF, Fc = 55Hz. Amplifier can be set to provide ideal ADC range of 3Vpp

Sinusoidal incoming signals are in the range of 10Hz to 40Hz that can be set for any duration in time. Ideally, in the range of 200ms to 300ms (can be changed)

In between these tones, there is dead space of (ideally) 50ms to 100ms (can be changed)

There would be a series of these tones with a total time of about 4, or possibly 5 seconds (can be processed after last signal arrives, ideally real-time)

I need the best amplitude and frequency accuracy possible, ideally about 2dB and 0.5Hz. The environment is noisy, but incoming signals can be kept 20dB over noise.

I am not looking for a code handout. Just a best approach algorithm and ideally some specifics such as sampling rate (is decimation required), buffer size, to use overlap-add/sliding window, signal processing chain - FIR required (how many taps), FFT method/size?

It may not be complex to many of you, but I haven't finished all the theory I need to try an implementation without spending a long time,(I have all the books but there's a lot of reading yet!) and the problem requires me to try for a solution now.

Thank you for your help!
Roger

As I understand, you want to sample a noisy tone signal at 20dB SNR and in 10Hz-40Hz range with best possible amplitude and frequency accuracy using 12 bit ADC of STM32F4 ARM+DSP chip.

To get the most frequency accuracy you need more samples to be buffered and large FFt sizes for a given samplig rate Fs. This is for frequencies on a given range.

Using high sampling rate (oversampling) may be advantageous for increased SNRQ. See Sigma-Delta converters and Noise shapers for details.

Size of FIR filters: As again it is a direct compromise between computing power and processing quality. STM32F4 is quite powerful chip at 168 Mhz operating frequency and 32 bit FPU.

For FFT, you can write your own code if you have time. I can provide you a butterfly implementation. Or you better use it from some free library.

You should select long fft sizes for better frequency resolution.

please have a look at this document: http://support.ircam.fr/docs/AudioSculpt/3.0/co/Window%20Size.html


What is the difference in between the different STFT implementations? (Overlap-shift with zero-padding, overlap-add, overlap-save, etc) I'm trying to understand STFT overlapping, why segments are concatenated and what are the consequences of this.

My implementation (from other questions and repositories found) of the STFT in Python is:

def stft(x, L, overlap, window): # Linear Spectrum (LS) [V]
opoints = int(L*overlap)
hopsize = L - opoints
X = np.array([(2 * np.abs(np.fft.fft(window*x[i:i+L])[:L//2])) / np.sum(window)
             for i in range(0,len(x)-L,hopsize)])
return X

where, x is the time-domain signal sampled at 44.1kHz, L is the FFT window size (4096), overlap is the overlapping factor [0, 1), and window is the window funtion used (Hanning). The absolute value of half the complex dft is scaled by 2/sum(win) so that the units of the linear spectrum are the same as for the input signal (V).

I used to understand the frequency resolution is fs/L and the width of the time slices or time resolution is L/fs, so that the number of time windows is N/L. However, if using overlapping (overlap!=0), the number of time windows increases ~(N/hopsize)-1.

Why are overlapped segments concatenated instead of added and what are the consequences of this regarding the time axis? (See picture below)
Is compressing the concatenated segments into the original duration of the signal correct? (i.e 3 segments of 1 second into 2 seconds, see picture below)

EDIT1:

Now, I do understand:

Overlap-shift and compression of the time axis of the overlapped segment into the original duration is one way to make a smoother spectrograph image
Overlap is needed in order to avoid lossy signals caused by the window

However, from @user31990 comment I still need to clarify:

What's the difference between the different methods/implementations?

Implementations I've found/understood:

Overlap and concatenate side-by-side the overlapped blocks (increasing the number of time windows and compressing the x time axis to the original duration)
Overlap adding the segments of the block that are overlapped together (the number of time windows doesn't increase = N/L)
Overlap save (zero-pad the beginning of the block and discard the transformed segment, the number of time windows doesn't increase)
Overlap shift with zero-pad

I don't totally understand the overlap shift with zero-pad, what is it doing?

My implementation is an example of the overlap and concatenate method.

Are This and Google's pytfd stft implementations examples of the overlap shift with zero-pad method?

EDIT2:

I was confused about the methods, question is answered/solved.

Zero-padding is a method that can be used: 1. To obtain an FFT of size N of a window whose size is smaller 2. To increment the size of the FFT, interpolating values in the frequency axis and obtaining smoother data.

Overlap-add and overlap-save are methods for the synthesis of an IFFT output in order to recover the original signal x.

Compressing the time axis of the overlapped windows into the original duration is one way to make a smoother spectrograph image.

Also, bell-shaped windows are lossy near the edges (after re-quantization), but a 50% overlap of Von Hann windows recovers all that information in each overlap.


Phasor notation of a complex number with apparently non-constant radius I'm asked to find the phasor of this signal:

$$ x_1[n] = 3\sin\left(\frac \pi7 n\right) + 4j\cos\left(\frac \pi7 n\right),\quad (0\le  n \le 20) \,.$$

The phasor is defined as:

$$ G z_{0}^{n}= Ae^{j\varphi}\cdot r^{n}\cdot e^{j\theta*n} = Ar^{n}[\cos(\theta n+\varphi)+ j\sin(\theta n+\varphi)] $$

But I really don't know how to make the conversion! I mean, one of them has a radius of 3, the other is 4. How can I do this?

Hint: a phasor is defined by its magnitude and phase. The magnitude can be obtained by computing

$$
    \sqrt{\left(\mathcal{Re}\left\{x_1[n]\right\}\right)^2
      + \left(\mathcal{Im}\left\{x_1[n]\right\}\right)^2}
$$

And the phase can be obtained with

$$
    \tan^{-1}{\frac{\mathcal{Im}\left\{x_1[n]\right\}}{\mathcal{Re}\left\{x_1[n]\right\}}}
$$


Kalman filter for position and velocity: introducing speed estimates Thanks to everyone who posted comments/answers to my query yesterday (Implementing a Kalman filter for position, velocity, acceleration ).  I've been looking at what was recommended, and in particular at both (a) the wikipedia example on one dimensional position and velocity and also another website that considers a similar thing.

Update 26-Apr-2013: the original question here contained some errors, related to the fact that I hadn't properly understood the the wikipedia example on one dimensional position and velocity.  With my improved understanding of what's going on, I've now redrafted the question and focused it more tightly.

Both examples that I refer to in the introductory paragraph above assume that it's only position that's measured.  However, neither example has any kind of calculation $(x_k-x_{k-1})/dt$ for speed.  For example, the Wikipedia example specifies the ${\bf H}$ matrix as ${\bf H} = [1\ \ \  0]$, which means that only position is input.  Focussing on the Wikipedia example, the state vector ${\bf x}_k$ of the Kalman filter contains position $x_k$and speed $\dot{x}_{k}$, i.e.

$$
\begin{align*}
\mathbf{x}_{k}  & =\left(\begin{array}[c]{c}x_{k}\\
\dot{x}_{k}\end{array}
\right)
\end{align*}
$$

Suppose the measurement of position at time $k$ is $\hat{x}_k$.  Then if the position and speed at time $k-1$ were $x_{k-1}$ and $\dot{x}_{k-1}$, and if $a$ is a constant acceleration that applies in the time interval $k-1$ to $k$, from the measurement of $\hat{x}$ it's possible to deduce a value for $a$ using the formula

$$
\hat{x}_k = x_{k-1} + \dot{x}_{k-1} dt + \frac{1}{2} a dt^2
$$

This implies that at time $k$, a measurement $\hat{\dot{x}}_k$ of the speed is given by

$$
\hat{\dot{x}}_k = \dot{x}_{k-1} + a dt = 2 \frac{\hat{x}_k - {x}_{k-1}}{dt} - \dot{x}_{k-1}
$$

All the quantities on the right hand side of that equation (i.e. $\hat{x}_k$, $x_{k-1}$ and $\dot{x}_{k-1}$) are normally distributed random variables with known means and standard deviations, so the $\bf R$ matrix for the measurement vector

$$
\begin{align*}
\mathbf{\hat{x}}_{k}  & =\left(\begin{array}[c]{c}\hat{x}_{k}\\
\hat{\dot{x}}_{k}\end{array}
\right)
\end{align*}
$$

can be calculated.  Is this a valid way of introducing speed estimates into the process?

Is this a valid way of introducing speed estimates into the process?

If you choose your state appropriately, then the speed estimates come "for free". See the derivation of the signal model below (for the simple 1-D case we've been looking at).

Signal Model, Take 2

So, we really need to agree on a signal model before we can move this forward.  From your edit, it looks like your model of the position, $x_k$, is:

$$
\begin{array}
xx_{k+1} &=& x_{k} + \dot{x}_{k} \Delta t + \frac{1}{2} a (\Delta t)^2\\
\dot{x}_{k+1} &=& \dot{x}_{k} + a \Delta t
\end{array}
$$

If our state is as before:
$$
\begin{align*}
\mathbf{x}_{k}  & =\left(\begin{array}[c]{c}x_{k}\\
\dot{x}_{k}\end{array}
\right)
\end{align*}
$$
then the state update equation is just:
$$
\mathbf{x}_{k+1} =
\left(\begin{array}[c]{c} 1\ \ \Delta t\\
0\ \ 1\end{array}
\right)
\mathbf{x}_{k} +
\left(\begin{array}[c]{c} \frac{(\Delta t)^2}{2} \\
\Delta t \end{array}
\right)
a_k
$$
where now our $a_k$ is the normally distributed acceleration.

That gives different $\mathbf{G}$ matrix from the previous version, but the $\mathbf{F}$ and $\mathbf{H}$ matrices should be the same.

If I implement this in scilab (sorry, no access to matlab), it looks like:

// Signal Model
DeltaT = 0.1;
F = [1 DeltaT; 0 1];
G = [DeltaT^2/2; DeltaT];
H = [1 0];

x0 = [0;0];
sigma_a = 0.1;

Q = sigma_a^2;
R = 0.1;

N = 1000;

a = rand(1,N,"normal")*sigma_a;

x_truth(:,1) = x0;
for t=1:N,
    x_truth(:,t+1) = F*x_truth(:,t) + G*a(t);
    y(t) = H*x_truth(:,t) + rand(1,1,"normal")*sqrt(R);
end

Then, I can apply the Kalman filter equations to this $y$ (the noisy measurements).

// Kalman Filter
p0 = 100*eye(2,2);

xx(:,1) = x0;
pp = p0;
pp_norm(1) = norm(pp);
for t=1:N,
    [x1,p1,x,p] = kalm(y(t),xx(:,t),pp,F,G,H,Q,R);
    xx(:,t+1) = x1;
    pp = p1;
    pp_norm(t+1) = norm(pp);
end

So we have our noisy measurements $y$, and we've applied the Kalman filter to them and used the same signal model to generate $y$ as we do to apply the Kalman filter (a pretty big assumption, sometimes!).

Then following plots show the result.

Plot 1: $y$ and $x_k$ versus time.

Plot 2: A zoomed view of the first few samples:

Plot 3: Something you never get in real life, the true position vs the state estimate of the position.

Plot 4: Something you also never get in real life, the true velocity vs the state estimate of the velocity.

Plot 5: The norm of the state covariance matrix (something you should always monitor in real life!).  Note that it very quickly goes from its initial very large value to something very small, so I've only shown the first few samples.

Plot 6: Plots of the error between the true position and velocity and their estimates.

If you study the case where the position measurements are exact, then you find that the Kalman udpate equations produce exact results for BOTH position and speed.  Mathematically it's straightforward to see why.  Using the same notation as the wikipedia article, exact measurements mean that $\mathbf{z}_{k+1}=x_{k+1}$. If you assume that the initial position and speed are known so that $\mathbf{P}_k=0$, then $\mathbf{P}_{k+1}^{-}=\mathbf{Q}$ and the Kalman gain matrix $\mathbf{K}_{k+1}$ is given by

$$
\mathbf{K}_{k+1} = \left(\begin{array}[c]{c}1\\
2/dt\end{array}
\right)
$$

This means that the Kalman update procedure produces

$$
\begin{align*}
\mathbf{\hat{x}}_{k+1} & = \mathbf{F}_{k+1}\mathbf{x}_k + \mathbf{K}_{k+1}\left(\mathbf{z}_{k+1} - \mathbf{H}_{k+1} \mathbf{F}_{k+1}\mathbf{x}_k\right)\\
& = \left(\begin{array}[c]{c}x_k + \dot{x}_k dt\\
\dot{x}_k\end{array}
\right) +
\left(\begin{array}[c]{c}1\\
2/dt\end{array}
\right)
\left(x_{k+1} - \left( x_k + \dot{x}_k dt\right)
\right)\\
& = \left(\begin{array}[c]{c}x_{k+1}\\
 2 \left(x_{k+1} - x_k \right) /dt - \dot{x}_k\end{array}
\right)
\end{align*}
$$

As you can see, the value for the speed is given by exactly the formula you were proposing to use for the speed estimate.  So although you couldn't see any kind of calculation $(x_k-x_{k-1})/dt$ for speed, in fact it is hidden in there after all.


Find the time lag for a wave suppose that we have a wave which has some peaks in some indexes like 2,12,25,... and otherwise zero.

I want to extract the time lag or lags between peaks of a signal. could anyone help me ?

The time lag is equal to the index difference times the sample period.  The sample period is:
$$
T_s = \frac{1}{f_s} = \frac{1}{16000} = 62.5\mu s
$$
Thus, the time lag from the first peak to the second is $(12 - 2) * T_s = 625\mu s$.  The time lag from the second peak to the third is $(25 - 12) * T_s = 812.5\mu s$.


Scaling down the modal value when thresholding I have a piece of code that I am translating and I am trying to understand the reasoning behind it.  The application here is thresholding medical (cellular) images.

The function in question takes a histogram and returns the lower threshold value.  The bit of code that I am translating mutates the histogram in a way that I do not fully understand.

It first finds the mode of the histogram.  It then finds the signal which appears most often which is not the mode (i.e., the signal which appears most often behind the mode). I'll just call that the "secondary signal".

If the modal value appears more than twice the number of times as secondary signal, the modal value is scaled down to 1.5 times the modal value.  For example:

mode = find_mode(hist);
second_max = 0;
for(int i = 0; i < hist.length; ++i) {
    if(hist[i] > second_max && i != mode) {
        second_max = hist[i];
    }
}

max = hist[mode];
if(max > secondary_max * 2 && secondary_max > 0) {
    hist[mode] = secondary_max * 1.5;
}

My problem is the "magic" behind that 1.5 number as well as the logic behind scaling only when the modal value is > 2 times the secondary signal.

I realize that this results in a higher threshold value and I assume that it is to compensate for noise (move the floor and next signal closer together), but I just don't know if these constants are well understood or if they were just chosen via experimentation on some unknown data set.

BTW, the developer is unreachable and I'm the only one here who knows anything about image processing... but I'm no expert, hence me coming here.

Most likely this is chosen by experimental results. Seems like for the algorithm to scceed it needs a value that is lower than the modal value. But not significantly lower (this is where the 1.5 scaling factor comes in).

I suggest trying the mid-range value and its variants to see how that works. You might be able to skip the scaling procedure.
mid-range = (max value + min value)/2


What's the difference between the Gabor and Morlet wavelets? The Gabor wavelet is a kind of the Gaussian modulated sinusoidal wave (source)

  Gabor wavelets are formed from two components, a complex sinusoidal carrier and a Gaussian envelope.  (source)

and

  In fact, the wavelet shown in Figure 2a (called the Morlet wavelet) is nothing more than a Sine wave (green curve in Figure 2b) multiplied by a Gaussian envelope (red curve). (source)

Are these just different names for the same thing?

Update:

Not to be confused with the "Gabor transform", which seems to be just another name for "STFT with a Gaussian window".  There's also the Gabor atom, which I guess is the same as the Gabor wavelet?

Since asking this on math.SE, I've also found terms like "Gabor/Morlet wavelet" and "Gabor-Morlet transform", implying that they are the same thing.

Also this has been asked before: Gabor transform/wavelet vs. Morlet wavelet but the answers are not clear to me.

The Gabor wavelet is basically the same thing.  It's apparently another name for the Modified Morlet wavelet. Quoting from Wavelets and Signal Processing:

  [The Modified Morlet wavelet] does not satisfy the admissibility
  condition but is nonetheless commonly used.  Sometimes this wavelet is
  called the "Gabor wavelet," but that term is improper because Gabor
  had nothing to do with wavelets.  He was ... one of the founders of
  time-frequency analysis.

That book is a collection of papers, and that paper ("The Wavelet Transform and Time-Frequency Analysis") is by Leon Cohen (of time-frequency distribution "Cohen class" fame), so I think it's reasonably authoritative.

At the very least, it sounds like the confusion is just a naming disagreement.  According to A Friendly Guide to Wavelets (pg. 114), Gabor was the first person to propose the use of Gaussian windows for time-frequency localizations, so his name tends to get attached anytime they're involved.


Derive velocity and acceleration from position I generate random spike data that represent the ouput from a rotary encoder. Here the output of my algorithm:

The first plot is position, then velocity and then acceleration. Even though I downsample my position vector and pass it to a average filter of 20 points, the derivative are very sensitive as you can see. I check the following question:
have position, want to calculate velocity and acceleration. Few option seem to be available , but I don't get the answer , the Savitzky-Golay Filters  is just a smoothing function, I don't get how he get the velocity, also what are the orther alternatives.

There are several ways to calculate the derivative of a stochastic signal. I will show you two methods.

First of all, adjust the output of a second order transfer function to your data: $${\omega_n^2}/({s^2+2\zeta\omega_n s+\omega_n^2})$$ Once you are satisfied with the fit, put an additional $s$ in the numerator and you will get the derivative of your signal: $${\omega_n^2\times s}/({s^2+2\zeta\omega_n s+\omega_n^2})$$ You can use a bilinear transformation and convert this solution into a discrete one, to implement in a software easily, if you like.

Another solution is to use the spectrum and calculate the fft of your signal. Then filter it vanishing the small values in the higher frequencies. Now, take the result and multiply it by $j\omega$. This is taken from the Fourier Transform Properties Table to calculate the derivative of a signal. Finally, redo the signal with ifft.

I added here a code to exemplify these aforementioned algorithms. The output picture is added as well at the end of the post. I hope it helps.

% \brief: this code exemplifies two techniques suitable for
%         calculating the first-order derivative of an
%         stochastic signal.
%
% \algorithm 1: LTI transfer function integrated with Tustin.
% \algorithm 2: FFT derivative property.
%
% \author:     Luciano Augusto Kruk
% \web:        www.kruk.eng.br

% some constants:
fs = 1000; % [Hz] sample rate
T  = 1/fs;
t  = 0:T:0.5;
nT = length(t);
% signal + colored noise:
[B,A] = butter(10,250*T*2);
s     = sin(2*pi*20*t)+filter(B,A,0.1*randn(1,nT));

figure(1); clf;
subplot(3,1,1);
plot(t,s);
title('signal s(t)');
grid on;

% first technique: second order transfer function
qsi  = 0.7;
wn   = 2*pi*60;
wnT2 = (wn*T)^2;
b    = 2*T*wn*wn;
a1   = 4 + (4*qsi*wn*T) + wnT2;
a2   = (2*wnT2)-8;
a3   = wnT2 - (4*qsi*wn*T) + 4;
dsdt = zeros(size(s));
for i = 3:nT
    dsdt(i) = (1/a1) * ...
        ((-a2*dsdt(i-1)) - (a3*dsdt(i-2)) + (b * (s(i) - s(i-2))));
end
subplot(3,1,2);
plot(t,dsdt, t,2*pi*20*cos(2*pi*20*t))
title('algorithm 1');
ylabel('ds/dt');
grid on;
set(gca, 'ylim', 150*[-1 1])
legend('algorithm', 'real')

% second technique: fft/ifft properties
s_fft = fft(s);
idx   = 1:(nT/2);
f_fft = (idx-1) * (fs/2) * (1/length(idx));
idx   = round((nT/2) + ((-nT/2.5):(nT/2.5)));
s_fft(idx) = 0;
w     = 2*pi*(((-nT/2):((nT/2)-1))/nT)*fs;
dsdt  = ifft(s_fft .* sqrt(-1) .* ifftshift(w));

subplot(3,1,3);
plot(t,dsdt, t,2*pi*20*cos(2*pi*20*t))
title('algorithm 2');
ylabel('ds/dt');
grid on;
set(gca, 'ylim', 150*[-1 1])
legend('algorithm', 'real')


Frequency analysis : relation between FFT size and sampling I am currently trying to perform an experiment given an audio signal. I am therefore sampling the audio frequency range (0 to 22kHz) with a 48k sampling rate.

When performing a FFT on my signal, I had the following question :

Is it equivalent to have a very small FFT size (for example 256), producing a 172 Hz resolution on my signal, in comparison with having a large FFT size (in a way more bars..) and performing a mean on my resulting frequencies in order to bring the number of bars down?

The two are not equivalent. Taking a larger FFT means you are using more time domain samples (unless you zero pad), in which case you have introduced new information compared to a smaller FFT size


The magnitude of taking discrete cosine transform of an image two times is similar to the original Consider the following:

im = double(imread('lena.bmp'));
subplot(1,3,1), imshow(im,[]), title('original');
t1 = dct2(im);
subplot(1,3,2), imshow(log(abs(t1)+1),[]), title('DCT transform');
t2 = dct2(t1);
subplot(1,3,3), imshow(t2,[]), title('DCT(DCT) transform');

the output is shown below:

Could anyone explain why the result of the second order DCT is similar to the original image?

Thanks in advance!

P.S. The same thing happens with every 2N'th DCT.

It's beacuse the inverse discrete fourier transform (DCT) is almost identical to the forward DCT. So taking twice the transform will be similar to the original signal. In fact if you provide which DCT type (DCT-I, DCT-II etc) you have used, one can show the effect more explicitly.


Fountain codes and LDPC codes I am having a small clarification with the difference between LDPC and fountain codes.

In LDPC codes, each parity bit is depenedent on numerous data bits. Isn't that similar to fountain codes as the encoded blocks are dependent on each other(I'm thinking of luby transform codes)?

Also, my professor was explaining that it doesn't matter if some blocks are missed, as long as it receives a subset of the encoded blocks. Fountain codes can decode the message because it just keeps listening to receive more information as there is enough redundancy built in.

My doubt is that, if it just keeps listening, when does the transmitter transmit the next set of messages if it's always sending redundancy to make up for what could be possibly deleted? Especially since there is no feedback.
Or is it that the message sizes are so large and that the encoding is done so that the redundancy is built throughout the message so that even if certain blocks are lost, we can still decode the information?
Lastly, fountain code applications are lossy connection, etc. In which case, why is this only applied to erasure channels?
Because if it is good at combating lossy connection, then it should also be able to effectively combat deletion codes right?

LDPC codes are a subset of linear channel codes. In terms of functionality, they are similar to other linear channel codes. However, they are constructed using a sparse bipartite graph.

Fountain codes are rateless codes. These type of codes are more useful in broadcast/multicast settings. If we use fixed-rate codes in such settings, we need to make sure that each and every packet is correctly received at all receivers. This can become inefficient since each broadcast can potentially be re-iterated  (depending on the re-transmission protocol) until all receivers receive all packets. Consider for example a case with 10 receivers, all have correctly received the packet #2 except one. This specific packet will be broadcast to all until the last receiver gets the packet. However, with a rateless code, the packet index is not a matter anymore (i.e. whether it is packet #2 or not). This is because a group of packets are selected at first and new encoded packets are generated randomly (hence, all randomly encoded packets look just like a new piece of information about the group and the index does not matter anymore). Any linearly independent subset of valid packets is enough to decode (acquire) the original group of packets.

Also I can see you have difficulty understanding acknowledge policy. We can have many different types of such policy. So based on the Ack protocol, the transmitter stops and listens for the ack from the receivers to find out their status in terms of the number or the index of missing/collected packets.


Understanding FFT of an Image I am having trouble identifying the frequency components of an image,
Here I simply generated, 256x256 binary-image.
a = [zeros(256,128) ones(256,128)];
imshow(a);

Taking the image FFT and shifting the Zero-frequency to the center,
f = fft2(a)
shft = fftshift(f);
imshow(log(shft)) %stretching

will result,

How do I identify how many frequency components are there and what are those frequencies?
Also if I calculate FFT of any 1-D sequence and plot that sequence there are certain values on negative, what does that mean? How could there be negative frequencies?

My first suggestion is that you understand FFT in 1 dimension before trying to interpret results in 2D.

The Discrete Fourier Transform (FFT is an implementation of DFT) is a complex transform: it transforms between 2 vectors complex vectors of size N.

So in the 1D case, you will get not only negative values, but complex values in general.

The same applies in 2D. The imshow function is probably taking the real part of the complex matrix (it is not clear in the imshow documentation).


Why does toning up a piano make it sound similar to that a string instrument I have a .wav file of a piano being played at 220 hertz. When I tone this up to 880 hertz, it begins to sound "stringy". However, if I were to use the piano plugin in fl studio and play it at 880 hertz, the piano would still sound normal.
My question is why a piano starts sounding stringy at higher frequencies(although I do know that a piano is in a way, a string instrument) as well as how to I can tone up a piano audio file so that it still sound like a piano(I am hoping to find a way to programmatically tone up a piano audio file as opposed to a program)
Link to file of a piano being played at 220 hertz(from fl studio): https://drive.google.com/file/d/1bUzBVXYrKS2Mzcg4VFEoBjIXqgGz1sFk/view?usp=sharing
Link to file of the above file being sped up to 880 hertz:
https://drive.google.com/file/d/1H16aydDErdBrrkDtHScfgITqrmMdi1I7/view?usp=sharing
Link to file of a piano being played a 880 hertz(from fl studio):
https://drive.google.com/file/d/1sEN4nvYS6QJeU0gtOruKvQGvDoZpzTx-/view?usp=sharing

The human ear, in collaboration with the brain, recognises what it is that is being perceived (as well as where it is coming from) within the first few ms of the sound waveform.
This "what is perceived" (is it a piano? is it an organ? is it a poggo-stick?) is called Timbre.
This modification of the attack phase is sometimes used deliberately.
Here, it is not the delay that is effecting the modification. Notice his hand on the volume knob when the guitar now starts sounding like a violin. It is clearer from 01:43 onwards in the "Without echo" part of the video. Here, he is killing the guitar's distorted "twang" as the string is excited and "brings it in" just a little after that so that you get more of the tone and less of the "twang". It then sounds like a violin. You might also notice on a piano, if you simply play the waveform after the attack, it sounds more like a decaying pipe-organ rather than a piano.
By "speeding up" the sound, you are distorting those crucial first few ms that the ear/brain are using to recognise the instrument. For example, the instrument now has a sharper attack than a piano. The reason it begins sounding "stringy" is another byproduct of the speeding up as the combination of sinusoids (and their relative phase) might now be tending towards a sawtooth wave that is typical of stringed instruments.
If you are looking to pitch shift but keep the duration of the sample the same, you need to be looking at techniques that attempt to preserve both the amplitude and phase spectrum of a sound.
Librosa implements this with pitch_shift() but other implementations are also available.
Hope this helps.


Why does a synchrosqueezed wavelet transform show oscillating behavior? This question came up in the context of the ssqueezepy library. As a basic experiment I did compute the synchrosqueezed wavelet transform of three basic signals:

A sine of 440 Hz.
A sine of 880 Hz.
A signal that mixes (1) and (2).

The result looks like this:

Full reproduction code: here
Basically the transform manages to perfectly localize signals (1) and (2), but for the mixed signal there is a surprising oscillation pattern. Zooming in a bit:

In general the issue depends on the choice of mu of the underlying wavelet. In this particular example, increasing mu slightly can help to make the transform of the mixed signal non-oscillating.
Nonetheless I'm wondering why synchrosqueezing leads to this oscillating interpretation of the signal in the first place?

This was interesting to figure out. The key lies in the phase transform, and how CWT interacts with own derivative upon insufficient component separation. Relevant are, and I'll be answering, the following:

What causes the wavy pattern?
Is the wavy pattern truly a sine, or a lookalike?
Why is only one band wavy, not both?
What can be done about it?

Answer split in two to separate meat from supplement, but much of insight's in latter; sections labeled in recommended reading order.
TL;DR: skip to Conclusion. Answer code.

0. How does the phase transform work?
A helpful requisite for following sections, covered here.

1. What causes the wavy pattern?
Begin by visualizing the phase transform, w:

The lower frequency's CWT and SSQ are nearly that of a pure tone, but the higher frequency is all wiggly. Let us plot one such wiggly row:

This sure looks sinusoidal. But is it? Stay tuned.
First: is the wavy SSQ explained by the w? Absolutely: w is wavy across timesteps, meaning for every column, we reassign CWT rows to a different SSQ row, in a wavy/sinusoidal pattern (see [complement] on reassignment). Our answer thus lies in understanding why w is wavy to begin with. Second answer begins here.

4. Inseparability: pre-transform perspective
Recall, the CWT and dCWT operations both take place in the frequency domain; we reach the lowest level of analysis by inspecting it directly. With xh = fft(pad(x)), psih = freq-domain wavelet, and dpsih = freq-domain wavelet with the derivative (* 1j * xi / dt), we plot both, overlapped and rescaled to fit in same graph, as well as the result of the multiplication - for our same row, and one on higher scale:

First, note for our original row scales[80] = 5.08, we indeed have a sum of two pure cisoids (analytic wavelet -> all negative frequencies are zero, and Morlet has pure-real DFT; cosine input -> pure real DFT; thus a DFT coefficient multiplies its complex exponential basis without cancellation), one at higher frequency having much greater amplitude.
Next, compare the cisoid amplitudes for psih vs dpsih; higher frequency dominates more for dpsih than for psih for both example scales. Coincidence? Well, they both do * xh, thus we expect the dpsih wavelet itself to skew toward right. Plotting,

this is confirmed. Does this happen for every scale? Yes: * xi = * linspace(0, pi, N) - taking the derivative involves multiplying by a line, thus skewing the entire wavelet toward higher frequencies, also breaking even symmetry (though not a lot).
All puzzle pieces are in. The skewness creates a difference in relative amplitudes at the two frequencies w.r.t. the non-derivative. The general form for dWx / Wx is then (still assuming purely real fft(Wx)):
$$
\frac{dW_x}{W_x} = \frac{a\exp{(j(\phi_1 + \pi/2))} + b\exp{(j(\phi_2 + \pi/2))}}{c\exp{(j\phi_1)} + d\exp{(j\phi_2)}}, \tag{4} \\
\phi_i = \omega_i \frac{n}{N}, \ \ n=[0, 1, ..., N - 1]
$$
($+\pi/2$ per the derivative). Let's first inspect interesting cases numerically, then derive an exact relation as proof (second answer).

8. What can be done about it?
Strictly speaking, nothing, unless the wavelet decays completely to zero eventually (untrue for many wavelets); as long as there's more than one frequency in the entire DFT of the input, there will be some form of interference. However, practically, such interference can be made negligible by:

Ensuring bands are sufficiently separated
Reducing wavelet frequency width (for Morlet, increase mu)

Let's inspect Morlet mu=5 vs mu=20:

mu=20 is clearly narrower; to get a better idea, single out a few rows:

Overlapped at same center frequency (~same row):

We thus predict mu=20 will show no wiggliness we can see. In fact, mu=20 is overkill, we'll settle for mu=10:

Can we also guess how to make both bands wiggly for small mu? Note, increasing scales lowers the wavelet's center frequency (shifts it left), but also reduces its frequency width; thus our wiggle odds improve toward lower scales (or higher frequencies). This turned out easier said than done; after some experimentation, I've had limited success, but maybe it can be done better:

9. Conclusion

What causes the wavy pattern? The frequency-domain wavelet and its derivative intercept the input bands at slightly different ratios per latter being skewed. The resulting phase transform is (imag. part of) a ratio of sums of complex sinusoids, which is of form $(\cos(\phi) + k_1)/(\cos(\phi) + k_2)$. The reassignment pattern dictated by w can then have a sine-like pattern across rows, which manifests as a sine-looking synchrosqueezed CWT.
Is the wavy pattern truly a sine, or a lookalike? Latter; it can never be exactly sine, but close to. The pattern is characterized by $(\cos(f_1 - f_2) + k_1)/(\cos(f_1 - f_2) + k_2)$, and has a period of $|f_1 - f_2|^{-1}$, where $f_1,f_2=$ input frequencies.
Why is only one band wavy, not both? With increasing scale (and decreasing center frequency), the frequency-domain wavelet moves to the left while shrinking in width; the resulting cisoid proportions work out such that dominantly sine reassignment pattern is toward the higher input frequency, where the wavelet is wider, but the lower one isn't devoid of such a pattern either.
What can be done about it? Separate input frequencies more, or reduce wavelet frequency width across all scales (for Morlet wavelet, raise mu).


Cascade filter realization equivalence? Given

$$H(z) = \frac{11 +4.6z^{-1} -26z^{-2}-3.75z^{-3}}{1-z^{-1}-8.75z^{-2}}$$

I'd like to know whether this realization:

Is equivalent to this one?

In short, is the direct form realization equivalent to the canonic one?

As MattL already mentioned, they are the same block diagrams as Direct Form II (canonical) implementation of a third order (on the FIR side) LCCDE.

However, the adders in the first block diagram are two input adders, while the the ones in the second block diagram are multiple input adders. Theoretically there is no difference.

In a popular DSP text, only the former block diagram is displayed, with sum (sigma) symbol replaced with a + sign.


How to reduce latency under mean filter for high noise? We have a position sensor that under some conditions receives some high frequency noise.  We can eliminate that very well with a simple mean filtering.

Unfortunately this causes too much lag when there is not any noise (at the beginning of the graph at s=90 to 95). Thus we have been experimenting with Savitzky-Golay filtering:

However, the results are not that great. The filtering under heavy noise is worse than mean filtering (Or the latency reduction is not that good). What else can we do?
Some more background:

The problem is an online problem, so we cannot retrospectively analyze the data and for instance shift the mean by half a period.
The 'regular' movement is performed by a human moving an object, the 'noise' is from a high frequency vibration.

Although this is an old entry and probably the author of the question has already solved the issue, I leave here my solution, just in case, it could help someone else.
As explained in the question, you need a kind of real-time "smoother" for the trajectory. You are going to have trade-offs: speed, latency, energy efficiency, etc.
I propose to use Alpha-Beta filters, Moving Average Filters or some kind of Adaptative Average filter. They should be faster, easy to implement and with less latency than Kalman filter or similar. I have made a couple of tests, trying to reproduce the type of data used in the question, and testing those filters with quick(and nasty) own implementations. The results are shown below in the next figures:

Cheers.
P.D. Probably, the vibration noise commented in the original question is not absolutely uncorrelated. The data I have generated is using random noise.


Converting PSD to dB scale I am taking the PSD of the output of an RTL-SDR with Welch's method.  The output is an array like: [1.34225796e-07   1.37964799e-07   1.24558405e-07, ....].  I'd like to convert this array to dB, so that I can find the point 3dB down from the center (which in this case is the first element of the array, because I am not using fftshift.)

How can I do this?  I don't understand the units given in the PSD in this case, nor how to convert them to dB (one problem follows the other).  I understand that the RTL-SDR output is not wrt some known unit (i.e. it's not dBm, but dBFS.)

A dB scale is simply a logarithmic scale.  Bels are base 10, and decibels are base $ 10^{1/10} $.  Because $log(a/b) = log(a) - log(b)$, a decibel scale converts ratios to intervals.  3 dB corresponds roughly to a factor of 2.  So, what you are really asking when you say looking for a dB value that is 3 dB down is equivalent to asking for where the raw value is half of the original raw value.

Added:  The units of the raw values don't matter since you are taking a ratio.  Simply take $ dB = 10 \cdot \log( raw ) $ to convert to a dB scale.

One thing you have to be careful of when taking the log of PSD values is that the PSD may have a zero value which makes the log undefined, but can be thought of as negative infinity.

Hope this helps,

Ced


Design of FIR Filters with Arbitrary Magnitude and Phase Responses I found many information in this thesis "Algorithms for the Constrained Design of Digital Filters with Arbitrary Magnitude and Phase Responses",i want to understand it.

This question is a continuation of a previous question(Complex Least Squares Approximation).

in Matt L's answer

$$H(e^{j\omega})=\sum_{n=0}^{N-1}h[n]e^{-jn\omega}=\mathbf{c}^H(\omega)\cdot \mathbf{h}\tag{1}$$

I want to approximate (1) by the following desired frequency response:

$$D(\omega)=M(\omega)e^{jP(\omega)}\tag{2}$$

My problem with this specifications :

Frequency domain specification is not symmetrical (linear phase filter with not symmetric impulse response),with order of filter and constant group delay, constant weighting function .

1-Before to start the design filter, must be know the problem is linear or no,does this have relation by 'h'?

In the case of FIR filters, we formulate the design problem using the complex least squares .The passband error function is

$$E(\mathbf{h})=\int_{\omega_1}^{\omega_2}| H(e^{j\omega})-D(e^{j\omega})|^2d\omega \tag{3}$$

How can i solve the integral in this case?or what is the steps which allow me to formulate an equation by least squares?

In complex leastsquares there is continuous and discrete approximation .However, often these integrals cannot be solved analytically.In
this case one either has to resort to numerical integration or formulate the problem a priori as a discrete approximation problem. What is the better for (3)?use the continuous or discrete approximation?

For the sake of clarity, let me point out that you do not approximate $(1)$ by $(2)$, but the other way around: you approximate a desired frequency response, i.e., your specification, by the filter's frequency response, and this is done by choosing the coefficients in the vector $\mathbf{h}$ such that some error measure (e.g., least squares) is minimized.

Since the frequency response $(1)$ is linear in $\mathbf{h}$, the problem is a linear least squares problem, which can be solved by solving a system of linear equations. This is true for the continuous error measure $(3)$ as well as for the discrete error measure

$$E(\mathbf{h})=\sum_{k=1}^K|H(e^{j\omega_k})-D(e^{j\omega_k})|^2\tag{1}$$

where $K$ is the number of discrete frequency points $\omega_k$. In general we require $K\gg N$, where $N$ is the number of FIR filter coefficients in the vector $\mathbf{h}$.

In the general case it is more convenient to use the discrete error measure $(1)$. There are only a few special cases in which the integral of the continuous error measure (Eq. $(3)$ in your question) can be solved analytically, otherwise you have to resort to numerical approximations anyway. And even if an analytical solution is available, the difference when minimizing the discrete error measure is negligible if the number of frequency points $K$ is sufficiently large.

Minimizing the discrete squared error $(1)$ is equivalent to solving the following overdetermined linear system of equations in a least squares sense:

$$H(e^{j\omega_k})\stackrel{!}{=}D(e^{j\omega_k}),\qquad k=1,2,\ldots,K\tag{2}$$

Formulation $(2)$ is convenient because many software packages such as Matlab have syntactically elegant and computationally efficient methods to solve such systems.

I've implemented a complex FIR filter design based on solving $(2)$ in the (extremely simple) Matlab/Octave function cfirls.m.


Is there a separate word for the time domain when working with digital signals? Typically the phrase "time domain" implies an analog signal x(t), but I've also seen it referring to a digital signal x[n].  Is there any phrase that refers to x[n], but not x(t)?

I use sampled-time domain or discrete-time domain.  I think I've heard other people use other terms, so they can chime in (discretely).


What is impulse response in simple AWGN channel? Assume a sender and receiver communicate through an AWGN channel. Let $x(t)$ be a transmitted signal and $y(t)$ be a received signal and $z(t)$ be the Gaussian noise. It is known that $y(t) = kx(t-\Delta t) + z(t)$ where $k$ is the channel gain (say due to attenuation) and $\Delta t$ is the delay. I was wondering what is the channel impulse response? I assume it should be $h(t) = k\delta(t - \Delta t) + z(t)$? Is it correct?

A linear system cannot add anything new to its input signal, so additive noise is modeled separately from any linear distortion of the signal. In your case with just scaling and delay, the channel impulse response is $h(t)=k\delta(t-\Delta t)$, and the complete channel is modeled by first filtering the input signal with an LTI system with impulse response $h(t)$ and then adding the noise.


Amplitude Response at greater than half the sampling frequency I am hoping to clear up some confusion I have. In a lab I am taking, we analyzed the amplitude response of a simple system. We found that as we increased the input signal frequency to greater than half the sampling frequency, the output signal began to flat line. The picture below show the corresponding Amplitude Response.

I am aware of the Nyquist theorem, and I thought that  the reason the output signal flatlined at frequencies greater than 8 kHz was because the sampling frequency was 16 kHz.

But during a second experiment of generating a sine wave, I discovered the concept of folding about the Nyquist frequency. Now in this case, increasing the generated signal frequency to greater than the Nyquist frequency caused the output signal to not decrease in amplitude, but to become symmetric about the Nyquist frequency.

E.g.

Generating a sine wave of 14 kHz produced an output of a sine wave of 2 kHz.

How are these two concepts related? It seems in the first example, going above the Nyquist frequency caused the output signal to flatline. In the other case, generating a signal above the Nyquist frequency didnt cause the output to flatline, but caused the output signal to only have a different output frequency.

This is a phenomenon called aliasing

https://en.wikipedia.org/wiki/Aliasing

Usually, there are anti-aliasing filters that will filter out signals with frequencies higher than fs/2.

However, in your case I think there isno anti-aliasing filter. Hence your 14 kHz-signal aliases to 2 kHz. Since f > fs/2, then apparent frequency = 16 kHz - 14 kHz = 2 kHz


Matlab SOS Filter in C++ I'm struggling to validate that my C++ implementation of a filter I designed in Matlab is correct.

My first issue is to make sure I'm implementing the sections in the right order. This is the output of the SOS matrix from fdatool:

0.5961   -1.1815    0.5961    1.0000   -1.9792    0.9884
0.5793   -1.1475    0.5793    1.0000   -1.9555    0.9649
0.6200   -1.2266    0.6200    1.0000   -1.9300    0.9398
0.5658   -1.1167    0.5658    1.0000   -1.9013    0.9119
0.4348   -0.8542    0.4348    1.0000   -1.8689    0.8806
0.3066   -0.5967    0.3066    1.0000   -1.8330    0.8459
0.1814   -0.3451    0.1814    1.0000   -1.7960    0.8104
0.0735   -0.1285    0.0735    1.0000   -1.7637    0.7794
0.0168   -0.0087    0.0168    1.0000   -1.7443    0.7609

I'm assuming that the sections should be processed from top to bottom, is that correct?

My second question is to do with scaling. I was scaling my filters using scale(Hd) directly, but I'm noticing that I get inconsistent results in my program. I have generated multiple filters that should output roughly the same data but the peak values are quite different. I noticed there is a ScaleValues property on the filter object, after scaling the values are mostly 1.0:

0.1047
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000
1.0000

I was all ready to assume these are applied to the output of each section, then I noticed there is 1 more scaling value than section. How are these values used?

Final question, given the options to scale, not to scale, various ways to scale, if I'm targeting a 64-bit CPU or maybe GPU doing my math in double-precision floats, is there an optimal way to configure scaling?

The ScaleValues property is easy to incorporate. Multiply the entire input signal by the first scale value, then the output of each section by subsequent values.


Indian Space Research Organization - Extent of spectral spread 1/3rd Octave Analysis I am a scientist at Indian Space Research Organization (ISRO) and will be representing my team. We work on Telemetry and Data Handling units for our Launch Vehicles.
We are preparing for our upcoming SSLV Launch (Small Satellite Launch Vehicle) and need some information to confirm if we are at par with the industry, or if we need to switch to better available options.
It would be really helpful if you share the name of the company which excel in software or equipment for 1/3rd octave analysis of acoustics signals, that you know of or have used and can vouch for best accuracy.
Please mention what is:

Order of band pass filter

Type of filter (Elliptic/butterworth etc.)

Order of Low pass filter (If one is being used in the design), type of low pass filter, its stop band attenuation and passband ripple.

Your help is duly appreciated.
Warm regards from our entire team at ISRO.

Rekhanshi Varma
Scientist/Engineer 'SC'
Vikram Sarabhai Space Center
Indian Space Research Organization (ISRO)
Office Ph: 04712564741 | Country Code: +91

That's not really a signal processing  question, but here are are some pointers.
3rd octave filters are well defined in ANSI S1.11-2004 and any decent commercial implementation should be fully compliant to the standard. The implementation is rather simple, so compliance isn't much of a problem. It's easy enough to do it yourself and their are ready-made packages in analysis products like Matlab https://www.mathworks.com/help/audio/ug/octave-band-and-fractional-octave-band-filters.html or measurement products like B&K https://www.bksv.com/en/instruments/handheld/sound-level-meter-apps/bz-7132-and-7131-frequency-analysis (no product endorsement intended, but these are reputable names).
One thing to keep in mind: In most real world acoustic problems, the difficult part is NOT the analysis but actually making the measurement and/or data capture. Where and how to place what type of microphones, when and how to capture with what signal, making sure that the signals are neither overdriven nor overly noisy, good signal to noise ratio at all frequencies of interest, good reproducibility, etc.. This requires fare more experience and expertise than the data analysis, so if you want to spend extra money, I recommend you spend it there.
Sound fields are very complicated 3-dimensional phenomena and a microphone only samples a single point in space. Moving a microphone by a few centimeters often makes much more of a difference than any analysis parameter.


What are the minimum knowledge required to begin the study of Digital Signal Processing? I am a student of computer science (not the university but technical school).
I am interested in the world of DSP, and I saw that it is closely connected with mathematics. I noticed for example that a lot is based upon the Fourier I'm slowly trying to learn.

I also want to buy the book "Understanding Digital Signal Processing", but I do not know what I would go against it, and if I can understand it. What I want to know is a short list of topics needed before we can deal with such a study that one should know, is listed for reference material.

Let me give you an example of what I mean (to be taken only indicative, since we do not know almost nothing of DSP)

Mathematics: Integral
Physics: Pressure
Science: Anatomy of the human ear, Anatomy of the human eye.

I know most of you are engineers with years of work experience and master behind it, but somewhere we must also begin ;)

I personally believe that the approach of studying something just in order to "feel ready" to study something else later on is not efficient and tends to overwhelm the student. (This of course only applies because I assume that you know what a computer scientist should know anyway: basic math, and - most importantly - how to think clearly and logically.) I could give you a whole list of mathematical topics which are important for doing DSP, but this would just paralyze you and would make you feel that there is so much you need to know that you won't even know where to start. Note that DSP has become such a huge topic in itself that there is hardly anyone who can claim that he knows even the basics of all its different sub-fields.

So what do I recommend to you after all? Get a good book on DSP and simply start studying it. If you come across problems that you cannot solve yourself, dig deeper by searching the internet, asking questions on this site, etc. In this way you will learn all necessary basics that you might be lacking. But the point is that only in this way will you find out what it actually is that you are lacking, so filling the gaps in your knowledge will be most efficient.

One good (and free!) book I would like to recommend to you is Introduction to Signal Processing by Orfanidis. It contains a lot of necessary basics and it is very practical. I think you will learn a lot in quite a short time. Also take a look at this online-book.

The answers to this question might also be relevant for you.


Interpreting the amplitude of signals in fourier transform I totally understand the concept of fourier transform, but one thing thats bothering me is the amplitude that we plot in the frequency domain. What does that amplitude of each frequency signifies? Is it the amplitude of that particular signal component in the input signal, as it is for the frequency? Is there any other meaning or any other way we can infer the amplitude?

Thank you for your time.

"I totally understand the concept of Fourier transform"

Lucky you if you really do. Some of us (me, in first place) don't (in totality). The Fourier transform (and its avatars) is a prototype for duality. Duality here means that you can represent a signal on some primal domain (time) onto a dual domain (here frequency). This transformation is meant to possess useful properties, to preserve signal information, and to add insight, like a somewhat easier interpretation.

Signals are traditionally represented as "amplitudes of something for each time". These amplitudes are not always absolute. But we can hope they are linearly related to the actual physical values. And that the linearity factor does not change over time. Under the above premises, sines and cosines, or complex exponentials (cisoids), are a very good way (perhaps the best) to model  the system. And in turn they can be used to represent the data, as a linear combinaison of shifted sines.

Since the sines are orthogonal, the squared magnitude of the coefficients weighting the sines is proportional to the energy of each corresponding sine frequency.

The proportional factor depends on how the coefficients are computed, often up to a scaling factor. The amplitude of a given frequency component can be directly computed with correlation. Note however that this coefficient should be taken with care, as outer conditions play a role: the finite horizon observation of a function/signal, amplitude discretization, noise, etc.


DTFT and Inverse DTFT Homework Problem I'm trying to solve this signals homework problem:

So for part a, since multiplication in the time domain is convolution in the frequency domain, I just used a DTFT table, found the DTFT for $\left(\frac12\right)^n$ and $\cos(\pi n/2)$, convolved them, and solved for $H(\Omega)$. I got the same answer they have in part a.

I don't quite understand their solution though, which leads me to believe there is a more intuitive way to think about this.

Anyways...for part b, you would think you'd use the DTFT of $\cos(n \pi/2)$, multiply that by $H(\Omega)$ and take that whole result and Inverse DTFT it back to the time domain, then solve for $y[n]$.

However, I can't get the math to work, and I can't seem to follow their solution.

Can anyone show me mathematically or intuitively how they get that final $\frac43\cos(n\pi/2)$?

Thanks,
Clint

The only thing you need to know is that an LTI system's response to a sinusoidal input signal is a sinusoidal output signal with the same frequency as the input signal, with its amplitude and phase modified by the magnitude and phase of the system's transfer function evaluated at the input frequency, or, in math:

$$x[n]=A\cos(\omega_0n)\Longrightarrow y[n]=A|H(\omega_0)|\cos(\omega_0n+\phi(\omega_0))\tag{1}$$

where $x[n]$ and $y[n]$ are the input and output signals, respectively, and $H(\omega)=|H(\omega)|e^{j\phi(\omega)}$ is the system's frequency response.

So if you got $H(\omega)$, you just need to evaluate it at the given input frequency $\omega_0=\pi/2$ and apply $(1)$. Splitting up the cosine into two complex exponentials, as suggested by the given solution, is not necessary (and not even easier I'd say).


FFT of white noise with different sampling rate I have a question on how I should interpret the white noise power level (noise floor) obtain from FFT for different hardware sampling rates. I realized if I sample the same noise at different rates (100 MHz, 10 MHz), even though I have down-converted (FIRfilter and decimate) them to the same, lower rate (1 MHz), I will obtain a different noise level if I FFT the data.

For example, with 1 sec worth of data, I obtain 100 million points sampled at 100 MHz. By down-converting it to 1 MHz and doing FFT, I obtain a noise floor of -80 dBm. If now I reduce my hardware sampling rate to 10 MHz and down-convert it to 1 MHz, the FFT result will give me -70 dBm of noise floor which is higher than when I sampled at 100 MHz.

May I know what the cause is?

My initial guess was the broad-band noise caused aliasing and they added up. Does it mean I should always filter out physically the noise above double my sampling frequency?

The variance of sampled white noise that is filtered by a brick-wall lowpass filter is reduced by a factor equal to the ratio of the filter cutoff frequency to 1/2 the sample-rate.  You start off with the same variance for the 100mhz and 10mhz cases, but the reduction factors are different due to the different sample-rates. Or to put it another way, the 100mhz-sampled signal has less power/Hz than the 10mhz sampled signal, so when you compare them over the same bandwidth, the filtered 100mhz signal has less power.


Short Time Fourier Transform analysis of newborn baby cry signal I am working on analysis of baby cry signal. fundamental frequency(Fo) of signal is already estimated using auto correlation and cepstrum analysis. cry signal is inputted in .wav format. now i want to do short time fourier transform analysis of same signal inorder to find variations in fundamental frequency over the time and plot its "pitch contour". i.e. Fo v/s t. its like finding Fo value for each small interval and plotting over time to see the variability. pls, can anybody help me with the matlab code?

Have a look at SAaC from LabROSA. You also might be interested in delta pitch without explicit pitch tracking.


Receiving frequency modulated signal I have two questions:

If the FM radio at the receiver mixes with a another frequency to bring down the carrier frequency, wouldn't that affect the message signal?
In today's digital world, why the FM radio still uses frequency modulation (FM) and not FSK?

If the FM radio at the receiver mixes with a another frequency to bring down the carrier frequency, wouldn't that affect the message signal?

This is the fundamental thing that all DSP bases on:

Any band-limited signal is 100% accurately representable by a baseband signal of the same bandwidth, and vice versa.

So, no, this doesn't affect the information at all!

In fact, this is the central result of Nyquist's works on the theory of continuous-time signals (with a few additions, it's equivalent to his sampling theorem).

  In today's digital world, why the FM radio still uses frequency modulation (FM) and not FSK?

Because of roughly 80 years of legacy devices! FSK is something that has a discrete alphabet of symbols, so it's only applicable to digital transmissions, really. Audio, as is, is an analog signal. So you'd first be converting your analog audio to digital, then modulating it, and then converting it back.

Fact-provenly, the spectral efficiency of digital modes is higher than that of analog FM – but you need digital receivers, which somehow still aren't that widespread.

That's a really long, and partly sad, story of badly introduced standards – there's multiple standards that strive to replace the analog radio broadcasts, and as far as I know, none has so far reached the proliferation that was once envision (note that this applies to terrestrial broadcasts only. For satellite radio, digital modes are overwhelmingly more popular):

DAB / DAB+: High-quality FM broadcast audio alternative; practically unused by anyone who isn't sitting in a high-end car with the "all-incluse" entertainment system. The difference between being able to receive DAB or not usually is something that the car dealerships earn a couple hundred Euros with. So not attractive at all, although German government has spent a lot of money on pushing that standard. Part of the problem was that early DAB modes had weak error correction, making the reception often worse than what people were used to via FM, receivers are still a lot more expensive, and DAB+, although it solves the quality issue, hasn't been adopted by everyone selling digital receivers. Coverage is mediocre, at best, compared to traditional FM worldwide.
DRM / DRM+: Digital Radio Mondiale, meant to replace the AM broadcast. Actually stands way better chances at seeing widespread adoption, because it solves actual problems, and the Indian Government alone has a 1MW transmitter to reach rural areas (read: most of India's area), so there's a lot of potential audience. Also possible candidate for broadcasting in Australia. There's about one receiver you can buy in Europe for that (source: been sitting in the same lab as the guys implementing the only FOSS transmitter in Europe).
add-on standards to "normal" FM: Just as stereo is an add-on to classical mono FM using more bandwidth, there are digital sideband and in-band modem standards (HD Radio etc) that provide additional functionality, like fine-grained weather info etc. Pretty much stillbirths, because where FM SNR is good enough, chances are high that people just get 2G/3G/4G coverage, offering actual advantages over just FM RDS.

By the way, none (to my knowledge) of these standards use FSK – if you're already doing digital radio in the 21st century, you'd typically go for something with higher spectral efficiency (most of these standards are OFDM things with PSK or QAM subcarriers).

So, to conclude, reasons for not adopting digital radio are:

literally billions of legacy devices that would need to be replaced
lack of advantage of digital modes ("FM sounds just fine, and I don't need more channels") or need ("99% of the population can get FM radio, and listener numbers are declining, why do we need new technologies")
actual technical disadvantages (e.g. DAB (not DAB+) was often worse than FM reception, and power consumption is still higher)
cost of new receivers high
bad time to market (broadcast standards take much longer to adapt than cellular standards, which now make broadcasting obsolete in metropolitan areas)
segmentation of market


Zero padding DFT intuition I'm trying to grasp some intuition about why zero-padding the time domain sequence $x[n]$ interpolates the frequency domain bins of the $DFT\{x[n]\} = X[k]$ and how does this relate to the $DTFT$ of $x[n]$.
Let's say we have $N$ point signal $x[n]$ and we pad it with 0's in the time domain as usual.
So, as we append more and more zeros to $x[n]$ and perform the $DFT$, we know that $DFT$ makes the assumption that the signal is periodic, hence the period tends to infinity, hence we get frequency samples (more and more and padding increaes) of the $N$ point $DTFT$ applied to the original $x[n]$ ?
I'd love some insights and intution please.
Thanks

Say a discrete signal $x(n)$ and its Z-transform $X(z)$, discrete time Fourier transform $X(e^{j\omega})$, and discrete Fourier transform $X(k)$. The relationships between DFT and ZT / DTFT are respectively:
$$
X(k) = X(z)|_{z=e^{j2\pi k/N}} \tag{1.1}
$$
$$
X(k) = X(e^{j\omega})|_{\omega=2\pi k/N} \tag{1.2}
$$
Since the original signal $x(n)$ can be perfectly reconstructed by IDFT, yes you can certainly obtain DTFT from DFT. However, you can't get DTFT by padding zeros as long as possible but by a technique called frequency domain interpolation. First is reconstruction of ZT from DFT:
$$
X(z) = \sum_{n=0}^{N-1}x(n)z^{-n}=\sum_{n=0}^{N-1}\Big[\frac{1}{N}\sum_{k=0}^{N-1}X(k)W_N^{-nk} \Big]z^{-n} \\ =  \frac{1}{N}\sum_{k=0}^{N-1}X(k)\Big[ \sum_{n=0}^{N-1}W_N^{-nk}z^{-n} \Big] \\ = \frac{1}{N}\sum_{k=0}^{N-1}X(k)\frac{1-W_N^{-Nk}z^{-N}}{1-W_N^{-k}z^{-1}} = \frac{1-z^{-N}}{N} \sum_{k=0}^{N-1}\frac{X(k)}{1-W_N^{-k}z^{-1}} \tag{2}
$$
where $W_N^{nk} = e^{j\frac{2\pi}{N}nk}$. Define the interpolation function
$$\varPhi_k(z) = \frac{1}{N}\frac{1-z^{-N}}{1-W_N^{-k}z^{-1}} = \frac{z^N-1}{Nz^{N-1}(z-W_N^{-k})} \tag{3}$$
and $X(z)$ can be written as
$$
X(z) = \sum_{k=0}^{N-1}X(k)\varPhi_k(z) \tag{4}
$$
Now let $z = e^{j\omega}$ we can derive that
$$
X(e^{j\omega}) = \sum_{k=0}^{N-1}X(k)\varPhi_k(e^{j\omega}) \tag{5}
$$
and the interpolation function
$$
\varPhi_k(e^{j\omega}) = \frac{1-e^{-j\omega N}}{N(1-W_N^{-k}e^{-j\omega})} = \frac{1-e^{-j\omega N}}{N(1-e^{-j(\omega-2\pi k/N)})}\\ = \frac{1}{N}\frac{\sin(N\omega/2)}{\sin\big[(\omega-2\pi k/N)/2\big]} e^{-j[\omega(N-1)/2 + k\pi/N]} \\
=\frac{1}{N}\frac{\sin(N\omega/2)}{\sin\big[(\omega-2\pi k/N)/2\big]} e^{jk\pi(N-1)/N} e^{-j(N-1)\omega/2} \tag{6}
$$
Let's write it in a more convinient form as
$$
\varPhi_k(e^{j\omega}) = \varPhi(\omega - 2\pi k/N) \tag{7}
$$
where $\varPhi(\omega)$ is the Fourier transform of a $N$-point rectangle window divided by $N$
$$
\varPhi(\omega) = \frac{1}{N} \frac{\sin(\omega N/2)}{\sin(\omega/2)} e^{-j(N-1)\omega/2} \tag{8}
$$
Finally, we get the equation of reconstructing DTFT by interpolating DFT
$$
X(e^{j\omega}) = \sum_{k=0}^{N-1}X(k)\varPhi(\omega-2\pi k/N) \tag{9}
$$


Irregularities in tone generated using FFT data I am implementing a simple program to generate tone in MATLAB using FFT/iFFT by zeroing all frequency bins except for the one frequency that I'd like to generate a tone in.

Sampling frequency is 16000 Hz.
I am using a 32 point FFT which gives me 16 frequency bins.

FFT data * Equalizer -> iFFT

If I want to generate a frequency of 1000 Hz, After FFT, using the Equalizer, I'd Zero all frequency bins except for frequency bin 2. This will generate a 1 kHz signal.

If I want to generate a frequency of 4000 Hz, After FFT, using the Equalizer, I'd Zero all frequency bins except for frequency bin 8. This will generate a 4 kHz signal.

I implemented this program.  Here are the results

For 1 kHz,

For 2500 Hz,

For 5500 Hz,

As you can see from the figure above, the 1 kHz is a clean signal. But as I try to generate higher frequencies, the signals are less smoother. The spectral figure shows that the correct frequency is generated, but i'm wondering if the irregularity is due to Harmonics or an other phenomena?

So how can I generate a more clean signal in the higher frequencies ? And if possible also the reasoning behind these 'irregularities'?

The closer a synthesized frequency gets to half the sample rate, the better the reconstruction filter ones needs to reproduce a pure waveform.  It looks like your plotting application has a poor reconstruction (upsampling to plot points) filter, causing the resulting plot to look irregular when the signal frequency gets much above about 1/8th the sample rate.  A graph using polynomial, spline, or Sinc interpolation for the plot points would likely look much cleaner.


What is discrete angular frequency? I am currently reading a paper about speech enhancement. The paper uses Spectral Subtraction Method. On that section, The paper stated that ω is the "discrete angular frequency index of the frames"

Question: What is the discrete angular frequency index? Is it part of Short Time Fourier Transform?

  Assume that the noisy speech $\boldsymbol{y[n]}$ can be expressed as $\boldsymbol{y[n] = s[n] + d[n]}$, where $\boldsymbol{s[n]}$ is the clean speech and $\boldsymbol{d[n]}$ is the additive noise. As the enhancement is carried out according to the frame, the above model can be expressed as
  $$y(n,k) = s(n,k) + d(n,k),\quad n=0,1,2,\ldots,(N-1);\quad k=1,2,\ldots,N\tag{1}$$
  Here $n$ is the discrete time index, $k$ is the frame number and $N$ is the length of the frame.
  $$y(\omega,k) = S(\omega,k) + D(\omega,k)\tag{2}$$
  Here $\omega$ is the discrete angular frequency index of the frames.

The paper is:

Navneet Upadhyay and Rahul Kumar Jaiswal: "Single Channel Speech Enhancement: Using Wiener Filtering with Recursive Noise Estimation", Procedia Computer Science, Volume 84, 2016, Pages 22-30, doi:10.1016/j.procs.2016.04.061.

$\omega$ is angular frequency in radians, $\omega =2\pi fT$, where $f$ is frequency (for example in Hz) and $T$ is sampling period (for example in seconds). This is revealed by Eq. 3:

  $$Y(\omega, k) = \textstyle\sum_{n=-\infty}^\infty y(n)w(k-n)e^{-j\omega n}\tag{3},$$

which represents discrete-time Fourier transform (DTFT) of windowed signal $y(n)w(k-n)$.

It's not common to call $\omega$ "discrete angular frequency index", which gives just 3 google hits.


What Does It Mean When Cross Correlation Between Two Images Occur at Negative Coordinates I am using Matlab normxcorr2 to calculate the cross correlation between several images, something like what has been done here: in MATLAB documentation. There are some of the examples that I get negative values for xoffSet and yoffSet, but I am not sure what it means when these values are negative.
Another issue, is practical since I need to give these xoffSet and yoffSet values as an input to coordinates of rectangle in an image, but since these are negative, the function throws an error. I googled it but no success in finding something that can explain such thing. I was thinking of just putting zero instead of negative values, but I'm not sure if that's a right approach.
I would appreciate if someone can explain about the meaning of such negative values and how I can substitute such negative values with appropriate quantities.

It means the best match to template happens outside the image.
For instance, let's say your template is 5 by 5.
And you got answer which is -1, -1. It means the part of the image which best matches you image is centered at [-1, -1] and you only have part of it in your image.
This is really an extreme case.
P. S.
If you share your data (2 Images) we'll be able to assist more.


Representing data in combination of QAM and BPSK I have a question regarding the data representation in different modulations. Assuming we are using MIMO system whith $4$ $T_x$ x $4$ $R_x$  which means four transmission antennas by four receiving antennas. My question, if three antennas are using QAM modulation, and the other antenna uses BPSK modulation. Is it possible the Receiving antennas can distinguish the data sent from each antennas ? I mean is it possible to detect the data sent from transmission antennas?

I'm asking that because in QAM modulation, data will be represented with real and imaginary parts, where in BPSK, data will be represented in real part only. So I see it logic that receiving antenna can detect the data perfectly. But I'm not sure.

Let's us take the simple case $2\times 2$ MIMO. Let $y_1$, and $y_2$ be the received signal at the first and second receive antennas, respectively. Then we have

$$y_1 = h_{11}x_1 + h_{21}x_2 + n_1$$

and

$$y_2 = h_{12}x_1 + h_{22}x_2 + n_2$$

where $h_{ij}$ is the channel coefficient from transmit antenna $i$ to receive antenna $j$, $x_i$ is the transmitted signal from transmit antenna $i$, and $n_i$ is the AWGN at receive antenna $i$.

We can write the above signals in matrix-vector form as

$$\mathbf{y} = \mathbf{H}\mathbf{x} + \mathbf{n}$$

Apparently, there is interference here. So, we need to do some sort of equalization. Let's use zero forcing (ZF) for simplicity. In this case you get

$$\mathbf{H}^{-1}\mathbf{y} = \mathbf{\tilde{x}}=\left[\begin{matrix} \tilde{x}_1\\\tilde{x}_2\end{matrix}\right]$$

where $$\tilde{x}_i = x_i + \tilde{n}_i$$

for $i = 1,\,2$.

Now, suppose that $x_1$ is 4-QAM, and $x_2$ is BPSK, both have the same symbol time. Will that be a problem? and why? They seem separable to me, and can be detected given that the receiver knows in advance that $x_1$ is drawn from a 4-QAM constellation, while $x_2$ from a BPSK constellation. Actually we do this in OFDM. It's called bit loading, where different subcarriers are modulated with different modulation sizes, and schemes when the CSI is know at the transmitter. Of course in OFDM the symbols are separable because of the orthogonality of the subcarriers, but we can somehow separate the signals in spatial multiplexing, although not perfectly.


How to derive a rate (e.g. water dripping) from an audio file I want to learn how to go about analyzing some audio files I have acquired.  Let's say the files contain some sort of repeating sound (e.g., water dripping), and I want to derive the pattern of the repeating sound.  I use R for programming and analysis and have briefly introduced myself to the tuneR and seewave packages, but I am open to any answers about how such an analysis should be conducted.

To make my question more concrete, let's say I have a 2 hour audio file of water dripping and I would like to be able to select some point or interval in the file and know the drip rate.  I will be pursuing further processing (e.g., finding times when the drip rate is most steady, fastest, slowest, etc.), but right now I am looking to learn how to get started with this task.

Specifics of my audio files:
.Wav files
Sampling rate: 44100 Hz
Bits: 16
Stereo

I briefly read about Fourier transforms, but due to my lack of knowledge and my rudimentary understanding of terminology, I am unclear if a Fourier transform would be appropriate starting place for this task.  I also briefly read about subcarriers, envelopes, and autocorrelation, but it quickly became clear to me that I could benefit greatly from some guidance before I start investing too much time learning about and introducing myself to the incorrect topics.

Beat detection is not what you are looking for.

Beat detection aims at looking for a periodic, sometimes implicit pulse in music signals. The events you are interested here are not periodic and every drip counts. What you are looking for is instead "Onset detection", to recognize each percussive event in the audio signal. There are ready to use software solution like aubio or sonic annotator to extract the timestamps at which percussive events occur in an audio file. From there, it is only a matter of running a sliding window over your output, count events within the window and divide by the duration of the window.


Why is the size of results from FFT half the size of the input, while that is not the case in image processing? When dealing with real data, Fourier Transform produces symmetric results and FFT algorithms discards half of it. An input with 1024 points will yield only 513 points in the output.

However, in image processing, when Fourier Transform is calculated, the size of the output image is same as that of the input image. Assume that the images are 2D. The (pixel) value of image function f(x, y) is still real. So, why do we still have an output image size same as input?

No.

Size of FFT in 1D is not half of the signal length. It's still the signal length (assuming FFT length was of the signal length). However for real data, due to conjugate symmetry, half of this is redundant and can be dicarded to reduce storage. Note however that, when it's required to process the FFT data, you have to construct the full length from the half and proceed.

This is also the reason why spectral analysis software display only positive frequencies for real data; i.e, negative frequencies will be a mirror copy (of the magnitude) of positive frequencies.

But you don't have to discard the half. You can just retain it.

For image processing, the symmetry of FFT for real input data again exist and if you wish you can also dicard half of the image FFT data. Whether this will be employed or not depends on the intentions of the package.


Time for a preamble (underwater acoustic communication) What is the maximum duration of a preamble for underwater acoustic communications? For example, can we have a preamble $>10 \textrm{ sec}$?

In communications, no matter whether it's radio, underwater acoustics or say powerline comms, there's multiple influences on the permissible duration of a symbol.

If you consider your preamble to be kind of a "super symbol" that you use for detection and equalization / phase correction, then you'd typically want to make sure that the channel coherency time is larger than that symbol's length; otherwise, your channel will change while you're still transmitting that preamble, and thus you will lose sensitivity detecting the preamble, or you'll learn about a channel state that simply is already "gone" the moment your preamble finishes.

However, if your receiver does e.g. OOK/power detection, then you can basically not go "too long" with your preamble – the longer it is, the more gain you'll have.

Now, I'm no expert in underwater communications, but all I know indicates that you might be dealing with a heavy multipath channel, and thus, frequency selectivity. Thus, your preamble detector might not only serve the purpose of detecting the start of a transmission, but also use the preamble to estimate the channel and feed the equalizer. In that case, as said above, a preamble that's longer than your channel coherency makes no sense.

Whether 10 s is longer than your channel is coherent is something I can't say. If you're doing communications in a container of liquid that doesn't change, then your coherency is probably "forever" and there's no upper limit. If you're doing communication across a stream with a lot of large debris rushing through quickly, then 10s might be far too long.


Time and frequency domain spreading in CDMA and MC-CDMA? For time domain CDMA transmitter, how the spreading sequences should look like? Is it necessary to be from $\{ -1,+1\}$ like PN sequences ? How will the codes be in frequency domain?
I am confused in spreading codes used for MC-CDMA and CDMA.
If I use a PN sequence $\in \{ -1,+1\}$ and then take FFT, it is MC-CDMA.
Now, If I directly use a DFT matrix as spreading codes, is it equivalent to MC-CDMA?

The spectrum of DS-CDMA without pulse shaping is equivalent to the sinc function. The Matlab code below illustrates the difference between DS and MC-CDMA. For MC you need the IFFT not the FFT function.

clear,clc

N=2;
% BPSK data
d=sign(randn(1,N));

% 7-chip m-sequence code
c=[1 1 1 -1 1 -1 -1];

% DS-CDMA
s1=kron(d,c);

% MC-CDMA
L=sqrt(length(c));
s2=L*ifft(d(1)*c);
s3=L*ifft(d(2)*c);
% etc, better to use a for loop for more symbols


Speech recognition library for detecting phrase in file I want to find all occurrences of a specific phrase in a set of audio files (it's a catch-phrase in a TV show). The phrase is pronounced by different people with different speed etc. What is the best approach to this problem?

Additional details: It's a one time thing. The end result should be a program that takes 30 hours of audio and produces timestamps of a particular short phrase. Obviously the training will be required. Problem is that while I'm a seasoned programmer I have pretty much no experience with speech recognition. I know that phoneme recognition engines exist will they get the job done? For now I need some pointers in the right direction.

The detection of a particular phrase in records sounded by different people may necessitate a complex speech recognition operation. Even if you train a good acoustic model, there is no guarantee that you will obtain a good recognition performance because you would also need some textual content which will constitute the language model.

Google speech recognition performs quite good and you can use its API instead of developing a system from scratch:
https://stackoverflow.com/questions/8830203/is-there-an-api-for-googles-speech-recognition-technology

You can send segments of voice data to the API and get the recognition results.


Removing values from FFT result same as filtering? I don't quite understand why the textbooks say it is impossible to implement an ideal low pass filter.

If I was to take the FFT of a discrete signal x[n], with Matlab's fft function I'd be returned with a complex sequence X[n] which represents the magnitude of each frequency from 0 to n/2 and is repeated. If I was to then say simply make all the lower frequency bins equal to zero and keep the remaining the same aren't I applying an ideal low pass filter. I am excluding all the frequencies I don't want and keeping the ones I do.

For example say I have the signal below

This is the fft, with the lower frequencies towards the ends.

I selectively remove the bins I do not want and keep the ones I do, note I don't actually have Matlab all the images are from a html page.

If you only have a finite number of samples and are using a finite length FFT, then you will end up with a finite number of FFT frequency result bins.  Each bin has a one-to-one relationship ONLY with frequencies that are exactly integer periodic in the FFT length.  Any other spectrum frequencies that are not exactly periodic in the FFT aperture will not have that one-to-one relationship, but instead will have their energy spattered across the entire FFT result (a sampled periodic Sinc or Dirichlet function).  Thus, zero-ing just the nearest neighbor frequency bin will not have the effect you think.

Zero-ing a bin is the same as adding an exactly periodic sinusoid of the same magnitude but of the opposite phase.  If what you are trying to filter out is spectrum of a very slightly different frequency from that of the bin center, then instead of zero-ing it out by exact cancellation, the result will be an added beat pattern.  Which is not the same as (FIR, IIR) filtering.

ADDED (in regards to ideal filtering) :

The basis vectors of a DFT are all orthogonal to one another.  Any other sinusoid (non-periodic-in-aperture) will be orthogonal to none of the basis vectors.  Which means it that that sinusoid will be represented, or have some energy end up in every DFT bin (because of not being orthogonal to any bin).

Therefore no subset of all the bins (a sequence of lower FFT/DFT bins for example) can exactly represent that signal.  Therefore zero-ing only a sequence of lower bins in an attempt to create a canceling signal will be imperfect, or non-ideal for that signal.  Thus a filter created by such a method cannot be ideal.

One might be able to get close for a signal near the center of a large sequence of filter bins.  But not for a signal near the edges of any sequence of bins, as a large portion of the signal's energy will be represented by those bins just outside the sharp edge of the filter rectangle.


White Gaussian Noise Spectrum and Power White Gaussian noise has constant power spectral density $N_0/2$. I know that the area under the power spectral density curve between two points gives the power of the signal between these two points.

If I want to know the power of a certain frequency in the signal (not in a range of frequencies), can we say that the power of each frequency in the signal is exactly $N_0/2$?

The total of power of additive white Gaussian noise is infinity, what does this mean? Is it reasonable to assume that the noise added to the signal have an infinite power?

If I want to know the power of a certain frequency in the signal (not in a range of frequencies), can we say that the power of each frequency in the signal is exactly No/2?

No, BUT: You mean the right thing, you just say it wrongly:
The Power Spectral Density is constant – a single frequency doesn't have any power; it has a "power per bandwidth"! To arrive at a power, you need to integrate the density over a non-zero mass of frequencies.
(This is kind of an important distinction to make – only infinitely long periodic signals, e.g. sine waves, have power at a single frequency; everything else has a "power distributed over frequencies".)

The total of power of additive white Gaussian noise is infinity, what does this mean? Is it reasonable to assume that the noise added to the signal have an infinite power?

Yes. Notice that you're never dealing with a truly white Gaussian noise in continuous-time systems (luckily for the universe, I might add); it's always approximately white for some bandwidth. Everything else is physically impossible – but rarely matters. Example: the thermal noise you can measure over a resistor is the classical example of white Gaussian noise in systems. However, it's not really white – the power density decreases at very high frequencies. But that's totally irrelevant to your observation – your measurement doesn't go into the terahertzes.
In time-discrete systems, things look different: for a sampled time-continuous stochastic signal (noise) to be white, it's sufficient that the original time-continuous signal had a constant PSD over a bandwidth. So, there's no physical problem in the time-continuous world. Since a discrete signal is just a sequence of numbers, there's no concern for "physicality" anyways.


Question about overlap-add @Jim Clay: No if the filter varies slowly, then the clicks go away. So I guess the filter needs to vary slowly to avoid audible clicking sounds.

I created a test in matlab where I'm comparing time-domain filtering with frequency-domain filtering. The two operations produce the same output, but I get clicks in the audio when I switch filter parameters periodically. How do I solve that?

Here's my script:

clc
close all
clear all

[x,fs] = wavread('audiosignal');
x = x';
numberOfFilterTaps = 256;
NFFT = numberOfFilterTaps * 2;
z1 = zeros(1,numberOfFilterTaps-1);
blockSize = 64;
numberOfSamples = length(x);
h1 = fir1(numberOfFilterTaps-1, 0.5);
h2 = fir1(numberOfFilterTaps-1, 0.9);
blockCounter = 0;
ola = zeros(1,NFFT);
enableThis = 1;
switchCounter = 0;

while 1
    idx1 = 1 + blockSize * blockCounter;
    idx2 = blockSize + blockSize * blockCounter;

    if (idx2 > numberOfSamples)
        break;
    end

    idx = idx1:idx2;

    sampleBlock = x(idx);

    % Switch filter periodically
    if (switchCounter == 0)
        switchCounter = 0;
        if (enableThis)
            h = h1;
            H = fft(h,NFFT);
            enableThis = 0;
        else
            enableThis = 1;
            h = h2;
            H = fft(h,NFFT);
        end
    end
    switchCounter = switchCounter +1;
    if (switchCounter > 60)
        switchCounter = 0;
    end

    % Method #1
    [y1(idx),z1] = filter(h,1,sampleBlock,z1);

    % Method #2
    tmp = real(ifft(fft(sampleBlock, NFFT) .* H));
    tmp = tmp + ola;
    ola = [tmp(blockSize+1:end) zeros(1,blockSize)];
    y2(idx) = tmp(1:blockSize);

    blockCounter = blockCounter + 1;
end

wavwrite(y1,fs,'output_y1');
wavwrite(y2,fs,'output_y2');

The main issue here is that you are implementing a time-variant filter. Both FIR filter and overlap algorithms are only valid if your filter is time-invariant.

If you need time variant filters, you need to solve two problems: continuity of the state variable and continuity of the output waveform. The easiest way to deal with this is to basically run both filters in parallel and then cross-fade between the outputs at switch time. The time constant of the cross fade determines what audibility you get, so this can be adjusted to the desired level.


Step response of a differentiating system Consider I can find the step response of a system with the following method:

We can easily find the step input of a system from its transfer function.  Given a system with input $x(t)$, output $y(t)$ and transfer function $H(s)$

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

the output with zero initial conditions (i.e., the zero state output) is simply given by

$$Y(s) = H(s)X(s)$$

so the unit step response, $Y_\gamma(s)$, is given by

$$ Y_\gamma(s) = \frac{1}{s} H(s)
$$

If my input is a "differentiating system" then what will be the step response ?

If your system is an ideal differentiator with input-output relation

$$y(t)=\frac{dx(t)}{dt}\tag{1}$$

then its transfer function is

$$H(s)=\frac{Y(s)}{X(s)}=s\tag{2}$$

From the equation in your question you obtain for its step response

$$A(s)=1\tag{3}$$

which in the time domain corresponds to a Dirac delta impulse:

$$a(t)=\delta(t)\tag{4}$$

This is also intuitively clear, because a step input has a derivative of zero everywhere except at $t=0$, where it is discontinuous.


Is there an adjective describing a filter with kernel that has zero mean? A linear filter with a kernel that has zero mean could be thought of as a "DC-rejecting" filter. Is there a better or more commonly used adjective for such a filter?

"Zero-Mean" is the word that's commonly used to describe signals and signals with a zero average. "This is a zero-mean filter."
If you really mean a filter that is specifically meant to cancel the DC component, a "DC blocker" is a name for that.


How to improve filter quality at low fs This my 1st post here.

I'm using this impulse invariant (MZT) based method to calculate coefficients for biquad digital filter (for RIAA and non-RIAA de-emph.):

double a0, a1, a2, b0, b1, b2;
double fs = 44100;
//timeconstants (case RIAA):
// frequency -> time conversion 1/(2*pi*fc) (= R*C)
//poles
double p1 = 3180e-6; // 1/(2*pi*50.05Hz)
double p2 = 75e-6;    // 2212Hz
//zeros
double z1 = 318e-6;  // 500.5Hz
double z2 = 0.0;       // 3.18e-6 for Neumann pole (50kHz)

double pole1= exp(-1.0/(fs*p1));
double pole2 = exp(-1.0/(fs*p2));
double zero1 = exp(-1.0/(fs*z1));
double zero2 = exp(-1.0/(fs*z2));

a0 = 1.0;   // = 1.0
a1 = -pole1 - pole2; // = -0.931176
a2 = pole1 * pole2; // = 0
b0 = 1.0; // = 1.0
b1 = -zero1 - zero2; // = -1.731986
b2 = zero1 * zero2; // = 0.733838

Known issue for MZT is magnitude error close Nyqvist at low sample rates - above filter results over -3dB error at 20kHz. I'm trying to find a way to fix this error. When using 96kHz sampling the error at 20kHz is around 0.5dB so the error (range 20Hz-20kHz) gets fixed by rising the sampling frequency high enough.

Not long ago Martin Visanec published an article "Matched Second Order Digital Filters" and he even kindly looked this case and tweaked his method suitable for above filter but the results he posted was not what I'm looking after (though, (without additional pole/zero pair) the result was less than ±0.3dB error in this above case).

Are there fitting methods one could consider to use for to improve the filter at low fs (for filtering system where coefficients are calculated in real time)?

Any thoughts and help with equations would be appreciated.

Juha

The impulse invariance method as well as the bilinear transform will give you relatively large errors close the Nyquist and there is nothing you can do about it. One thing you could do is use a filter design method that approximates the given analog response. I've tried this for the RIAA de-emphasis fitler with an (unpublished) algorithm I came up with some time ago, and the figure below shows the result:

The digital filter is a second-order IIR filter, the sampling frequency is $f_s=44100\text{ Hz}$, and the error is much less than 1dB over the whole frequency range up to Nyquist. The price you pay for it is 20 iterations of an approximation method, which solves a $5\times 5$ system of linear equations in each iteration.

EDIT:

If you add one more pole and one more zero (i.e. your discrete-time filter has 3 poles and 3 zeros), then the error goes down to less than 0.1 dB over the whole frequency range up to Nyquist:


GNU Radio loop bandwidth normalization From GNU Radio documentation, the control_loop block proposes loop bandwidth values in the range $[\frac{2\pi}{200}, \frac{2\pi}{100}]$ in radians per sample.

Some blocks, like the new timing symbol sync block recommends a value around $2\pi\cdot 0.040$. Given sample rate $R_{sample}$ and symbol rate $R_{symbol}$, the symbol rate can be expressed in radians per sample as

$$2\pi\frac{ R_{symbol}}{R_{sample}} = \frac{2\pi}{N}$$

where $N$ is the oversampling ratio is samples per symbol. From literature, the loop bandwidth around a few percentages of the symbol rate (say 3% for example) is recommended.

I have a problem in relating the recommended control_loop values to the values suggested in the literature. Suppose for example the signal is oversampled by a factor of $N = 100$ samples per symbol.

A loop bandwidth value of say $\frac{2\pi}{100}$ will be equivalent to 100% of the symbol rate, which doesn't seem to be logical. The same loop bandwidth will be equivalent to 3% for an oversampling factor of $N = 3$.

I think that the ratio between the loop bandwidth and the  symbol rate depends on $N$. For example, in the figure below, I don't think a loop bandwidth of $\frac{2\pi}{100}$ would mean the same thing to the FLL and the PFB ($N = 4\,\text{sps}$) as it will to the Costas loop ($N = 1\,\text{sps}$).

Why don't the recommended values in control_loop take N into account?
How can the loop bandwidth in GNU Radio synchronization be configured
as a percentage of the symbol rate?

I can only answer your second question:

"How can the loop bandwidth in GNU Radio synchronization be configured as a percentage of the symbol rate?"

The tracking loop in the symbol synchronizer block operates at the symbol rate, estimating timing error and making a correction once per symbol.  So the sample rate of the error signal from the TED is at approximately 1 sample/symbol.  (I say approximately, because technically the block is continually estimating symbol clock error and adjusting its estimate of the symbol clock period, but its objective, nominal operating rate is 1 error sample/symbol.)

This means that in terms of normalized digital radian frequency, $\omega T_s$, for the loop filter, the symbol frequency corresponds to $2\pi$ radians/symbol.

The loop bandwidth parameter of the symbol synchronizer block is expected in units of the normalized digital radian frequency.  So when one specifies $\dfrac{2\pi}{200}$ for $\omega_n T_s$, one is specifying an approximate one-sided loop filter bandwidth of $2\pi \cdot 0.5\%$ radians/symbol, a one-sided filter bandwidth that is $0.5\%$ of the symbol frequency.

To express the loop bandwidth input of the symbol synchronizer block as a percentage of the symbol frequency, just use

$$\mbox{(one-sided) loop bandwidth} = \omega_n T_s = 2\pi \cdot n\%$$


Finding delay where squared error between signal and data is minimal I have a set of points  which correspond to part of a known signal. What I want to do is estimate the delay between the set of points and the known signal.
In this case, I don't think cross correlation is appropriate, since the approach should work even if only a small number of points are available. I think this is an optimization problem that minimizes the error between the data points and the known signal with respect to the delay.

Is there any operator similar to cross correlation that achieves this, or does it need to be solved by going through each delay seperately?

Just doing something like:
$$
e_\delta = \sum_{n=0}^{N-1} \left| x[n-\delta] - x_\delta[n]\right|^2
$$
for several values of $\delta$ and choosing
$$
\hat{\delta} = \arg \min_{\delta} e_\delta
$$
(i.e. the value of delta that minimizes $e_\delta$) will probably be good enough.
Example signals below:

gives the following error plot

The error plot seems to give the right value of about 5.5 samples for this contrived example data.

Code Only Below
from numpy import log10, asarray, polyfit, ceil, arange, exp, sin, pi, log, random, sum, diff
import matplotlib.pyplot as plt
from scipy.signal import kaiserord, lfilter, firwin, freqz

T = 100
time_period = list(arange(1,T))

K1 = 10
K2 = -10
tau = 200
measurement = [K1 + K2*exp(-x/tau) for x in time_period]

fir_filter = [0,0,0,0,0,0.5,0.5,0,0,0,0]
channel_1 = lfilter(fir_filter,1, measurement)

max_delay = 20

def error_calculation(signal, delayed_signal, delay):
    error = [ (signal[time-delay]-delayed_signal[time])**2 for time in list(arange(delay,T-delay)) ]
    return sum(error)

plt.figure(0)
plt.figure(figsize=(20,20))
plt.plot(time_period, measurement)
plt.plot(time_period, channel_1)
error = []
for delay in arange(1,max_delay):
    error.append(error_calculation(measurement,channel_1,delay))

plt.figure(1)
plt.figure(figsize=(20,20))
plt.plot(arange(1,max_delay), error)


Determining linearity, causality, memory, and time invariability from a picture of the signal So, this is a picture of saccadic eye movement reflected by an EOG. How can I determine if the eye movement system is linear, causal, memoryless, and time invariant from this displayed waveform? Please explain.

Basically you can't. In basic information theory a system has an input and output and the terms you are asking about are defined in terms of the relationship between output and input. It's difficult to determine what the input to saccadic eye movements would be. It's mostly an unconscious process controlled by the neural periphery.

In order to answer your question you would first define an appropriate input signal and than look at a screen shot that shows both input and output simultaneously.

FWIW you could potentially argue that it's causal simply because it is a real physical system and all physcial systems are causal.


harmonics of a signal and FFT Consider that we have a discrete signal of finite length. How can we find the amplitude and phase corresponding to different harmonics of this signal in Matlab?

Thanks.

As Florent told: you must take fourier transform from your signal, by fft, then you will get your spectrum in normalized frequency, actual frequency is correspond to your sampling rate, then you can see amplitude and phase of each frequency(or harmonic) by these command: abs and angle.
For example, assume x is signal and T is time interval between your samples and

x=[10 -10 10 -10 10 -10 -10 10 10 -10 -10]
X=fft(x);%you can set frequency resolution by setting fft(x,resolution)
Fs=1000;
f=linspace(1,Fs,length(X));
plot(f,angle(X))
figure
plot(f,abs(X))