The minimal signal to noise ratio (SNR) for people to understand a speech in the noisy background This question seems easy to find an answer but I got lost in the end. I am trying to find the threshold for people to understand speech in a noisy environment. For example, assume the background noise is $\mathbf n$ and the speech signal is $\mathbf s$. I wanna know the threshold thr for the speech being recognizable by a human when:
$$
10\log\left(\frac{\lVert \mathbf s\rVert}{\lVert \mathbf n\rVert}\right) > \text{threshold}\implies\quad\text{so the speech $\mathbf s$ is recognizable by human}$$

The most relevant scientific fact I found is the absolute hearing threshold, such as the ISO 226.

However, I got confused to interpret this ISO 226 figure because the $y$-axis is in SPL (dB). Based on my understanding, the SPL is also defined based on the minimal hearing threshold, which is kind of weird to understand how much signal strength (of a speech) is actually necessary for people to recognize.

Can I say the signal strength of a $10\textrm{ kHz}$ signal needs to be $\approx 5\textrm{ dB}$ higher than the noise floor for being recognized by the human based on the red lines shown in this ISO 226 figure?
If it is true, how can I know the threshold for a speech signal covering multiple frequencies?

The Equi-loudness contours, chart provided at your link, is for hearing thresholds for isolated tones. Assuming speech is a sum of tones is a wrong assumption. But even if we go ahead with this assumption, it does not help to compute SNR for speech understanding under noise because of the following reasons.

a. Speech is a non-stationary signal implying that the spectral content evolves over time. In other words, the energy of the different frequencies in the spectrum from 0-8 kHz will change depending on the uttered phonemes (making up the speech). An example is shown in Fig.1 below.

Now, coming to your question on minimal SNR for people to understand speech, the notion of SNR is not well defined for speech as the power spectral density of speech is varying with time. Although you have provided a definition for SNR you can see that it just a ratio of full signal energy. Instead, we perceive speech as a function of time.

b. Speech perception in noise can be analysed from at least two aspects - intelligibility and quality. I assume by understanding speech you mean intelligibility, that is you can make out the phonemes uttered. Intelligibility in noise is dependent on the power spectral density of noise. Different acoustic noise have different spectrum, and affect the intelligibility to different extent. Like a tone interference may not effect speech intelligibility, but something like babble noise (example, in a party where many people are talking in the background) can easily impact the intelligibility.

c. Having recognized above nuances with respect to speech perception, the Equi-loudness contours specify the threshold of hearing for tone stimuli. The threshold of hearing is obtained in a quite environment (or absence of any other sound). The 0 dB is with respect to threshold of 1 kHz tone. Hence, the curves are of not much help to compute threshold in noise. However, these have been used in designing equilizers in sound playback systems.

d. In case you play two tones and ask the question of when is one sound masked by the other then this is termed as "masking threshold". Masking thresholds are used in MP3 codec to improve coding gains.


What Is the Bilateral Filter Category: LPF, HPF, BPF or BSF? I`m reading about Bilateral filter with images, I knew that it uses Gaussian filter (which is a LPF) as a domain filter plus a range filter and it is used to preserve edges(high freq component).

Can we consider Bilateral filter as a LPF ? .

I searched some articles:
1- Which approach is better for decomposing a image into high frequency and low frequency component?

The answer was that bilateral filter attenuates ""medium frequency "

2- How to Extract High Frequency and Low Frequency Component Using Bilateral Filter?

The answer was to consider bilateral filter as a LPF.

Neither. To me, filter classes using the notion of frequency bands (low-pass, high-pass, etc.) can be used  safely in the linear case. And the bilateral filter is nonlinear. Edges are not really high-frequency: they often have sharp variations across the edge, but slow variation along it. Hence linear directional filters often have a derivative part and and orthogonal smoothing part.
I would consider the bilateral filter as an edge-preserving smoother, a broad and somewhat imprecise class.
Suggested reading: Fast and Provably Accurate Bilateral Filtering, 2016

The bilateral filter is a non-linear filter that uses a range
filter along with a spatial filter to perform edge-preserving
smoothing of images. A direct computation of the bilateral filter
requires O(S) operations per pixel, where S is the size of the support
of the spatial filter. In this paper, we present a fast and provably
accurate algorithm for approximating the bilateral filter when the
range kernel is Gaussian. In particular, for box and Gaussian spatial
filters, the proposed algorithm can cut down the complexity to O(1)
per pixel for any arbitrary S. The algorithm has a simple
implementation involving N+1 spatial filterings, where N is the
approximation order. We give a detailed analysis of the filtering
accuracy that can be achieved by the proposed approximation in
relation to the target bilateral filter. This allows us to estimate
the order N required to obtain a given accuracy. We also present
comprehensive numerical results to demonstrate that the proposed
algorithm is competitive with the state-of-the-art methods in terms of
speed and accuracy.


Alternative to BIBO stability of a system In DSP textbooks a system is stable in the BIBO (Bounded-Input, Bounded Output) sense if and only if every bounded input sequence produces a bounded output sequence. After stating this definition stable always means BIBO-stable.

I'd like to know if there exists other form of stability for a system?

For linear systems, BIBO stability is the most useful and practical criterion. For systems described by rational transfer functions it coincides with the condition that all transfer function poles must be located in the left half of the $s$-plane (for continuous-time systems), or inside the unit circle of the $z$-plane (for discrete-time systems). For linear systems, I do not know of any other stability criterion that makes any sense and does not coincide with the BIBO-stability criterion. Since most basic DSP texts focus on linear (and especially time-invariant) systems, you will only find the notion of BIBO-stability there.

Things look a bit different for non-linear systems. There you have several reasonable stability criteria that do not all lead to the same basic criterion as is the case for linear systems. One important notion for non-linear systems is input-to-state stability, which basically means that for zero input, the system is stable about its zero state, and that well-behaved and bounded input signals produce a bounded state trajectory. This article reviews some of these concepts.

But if you're mainly interested in linear systems, BIBO-stability is all you need.


What is the difference between a controller and a compensator? I recently revisited some texts on control theory. In some textbooks the term compensator is used, seemingly interchangeably with controller (see e.g. The Control Handbook by Levine, chapter 14 on classical control system design methods).

This confuses me. In my field (aerospace engineering), we always talked about controllers. Since it's some time since my days at university I may have missed/forgotten some subtle differences...

Please help me to understand the difference:

What is the (technical) difference between a compensator and a controller?
When do you call it one way or the other?

As the book's title is control theory; a controller is the general broad name given to any piece of hardware (or sowftware) that aims to control some processes, in your aerodynamics field probably the flight process; by menas of receiving signals from sensors and sending commands to propellers - actuators, and in between is the control hardware or software...

A compensator is a specific control-box-tool to correct (or compensate) some lacking aspect of a system property under concern, such as phase, lag, magnitude etc... It is generally useful as part of a bigger control system, rather than itself being a controller except some simple applications.

Obviously there are many different control system examples in addition to a compensator...


Clustering algorithm I am looking for some suggestions regarding clustering.
Lets say
I have 5 subjects (1, 2, 3, 4 and 5)

1 says  2 and 3 are similar to 1
2 says 1 is similar 2
3 says 5 is similar to 3
4 says no one is similar to 4
5 says 3 is similar to 5

So we have a matrix like

1 - 2 3,
2 - 1,
3 - 5,
4 - 0,
5 - 3

I need to cluster saying

cluster1 - 1 and 2
cluster2 - 4
cluster3 - 3 and 5

any suggestions to do clustering based on this type of data?
most clustering are based on calculating the distance difference. But here I think that is not the case.
Thank you!

You are not really looking for a clustering algorithm in the common meaning of the word, but more of a Graph Theory concept.

  So we have a matrix like

That is not a matrix. That is an adjacency list that states nodes adjacent to others.

The "matrix" you might be after is the adjacency matrix. For a simple graph (a graph that does not contain nodes with self-connections), $G(V, E)$ where $V$ is the set of vectors and $E$ is the set of edges, the adjacency matrix $A$ tells you if there is a connection between two vertices in $V$. For example, $a_{2,3}$ would tell you if there is a connection between vertices $v_2, v_3 \in V$.

Applied to your problem, the adjacency matrix looks like this:

$$V = \{ 1,2,3,4,5\}$$
$$ A = \begin{bmatrix}0 1 1 0 0 \\ 1 0 0 0 0 \\ 0 0 0 0 1 \\ 0 0 0 0 0 \\ 0 0 1 0 0\end{bmatrix} $$

A single graph can be composed of a number of components. In your case, node 4 is completely unconnected and it could be considered as a component.

If your graph was undirected, the task of finding connected components would be very simple. But, your graph is directed. This is denoted by the lack of symmetry of $A$ around the main diagonal ($A_{1,3} \ne A_{3,1}$). In this case, you would be looking for strongly connected component algorithms.

A Python module you can use to compute both is Networkx. With functions for connected components and strongly connected components here and here respectively.

This concept of clustering is very interesting but (at first view) quite different than the traditional clustering of data points in some space based on "similarity". See for example this type of clustering and you will see why you need a concept of "distance". In your type of problem, this would be equivalent to having the nodes also report HOW SIMILAR they are to each other. Then, in addition to direction, you would also have a weight (similarity) associated with each edge. In that case, things get a little bit more complicated.

Hope this helps.


what is window size? i am working on simulation of an article "A Real-time ECG Data Compression Algorithm for a Digital Holter System". This algorithm has five steps "Difference=> ECG Signal=>Classification=> DCT=> Window Filter => Huffman" . it described window filter as "In order to implement a high compression ratio, the DCT window size, which is chosen to be proportional to the DCT result, is adjusted from 15% to 100%.". i do not understand this step and how can i do it? how can i change window size?
please help me
thanx

That paper is very poorly written. It's no wonder you have a hard time understanding it; I certainly did!

What I think is happening is that several of the DCT coefficients are being thrown out (windowed). So if the window is 100%, then all the DCT coefficients are kept; if the window is 30% then 70% of the coefficients are thrown out (zeroed).

See Figure 11 of that paper:

The left hand signal are the DCT coefficients. The right-hand graph is the applied window. The two are multiplied together.


Sample and reconstruct a real exponential (just one period) I have a function with an equation:

$$C = 1.6925\left( e^{-0.136t}-e^{-1.192t}\right) $$

Where $C$ is real and $t$ represents time in hours. Beneath is the representation of my function.

I am trying to find the best sampling step (based on Nyquist-Shannon theorem) to sample and reconstruct this function using a Fourier Transform.

After using a Fast Fourier Transform algorithm in MATLAB (fft), by setting:

x= fft(c);
xmag=abs(x); (amplitude estimation)

I find the below spectrum for the function for a $t$ between $0$ and $100$ sampled at $1000 \text{Hz}$  , with the highest amplitude at:

bin(1) = 1.1025e+04

and the second highest at bin(100001) with an amplitude of 9.9946e+03. the bin number 100001 represents also the length of xmag.

So my question is how would you interpret the spectrum? Specifically I would like to know how to determine the highest frequency at which I can sample this function and reconstruct it without aliasing?

First, since $t>0$ in your case, you can write your function as

\begin{equation}
C(t)=1.6925\left(\exp^{-0.136t}- \exp^{-1.192t} \right) u(t)
\end{equation}
where $u(t)$ is a unit-step function. Then, denote the FT of $C(t)$ as $C(F)$, which is given by
\begin{equation}
C(f)=1.6925\left(\frac{1}{0.136+2\pi jf}- \frac{1}{1.192+ 2\pi jf} \right)
\end{equation}

The spectrum $C(f)$ is plotted below. We can see that $C(f)$ is almost zero after 5 Hz. You can choose the sampling frequency according to the Nyquist criterion as $f_{\max} \geq 10$ Hz. You have taken a very high sampling rate (1000 Hz) that is 100 times more than the required rate. When I took 10 Hz to sample the signal above and took the FFT I got the results as follows  (in order time domain $C(t)$, frequency domain signal $C(f)$, frequency domain domain signal obtained via FFT, that is, FFT(C(t)). Note the inherent periodicity of the FFT spectrum)

And here is the MATLAB/Octave code for your reference..

clear all; clc;close all;

%% Time-domain signal
Ts=0.1; % Sampling time (obtained from the frequency domain plot)
t=0:Ts:60; % in sec
C=1.6925* (exp(-0.136*t)-exp(-1.192*t));
figure,plot(t,C);
title('Time domain signal');
xlabel('Time in sec')

%Frequency domain signal
f=0:.1:10;
F1=1.6925*((1./(0.136+(j*2*pi*f)))-(1./(1.192+(j*2*pi*f))));
figure,plot(f,abs(F1));
title('Frequency domain signal');
xlabel('Frequency in Hz')
% Via FFT
F2=fft(C);
L=length(F2);
figure,plot((0:L-1)*(1/(L*Ts)),abs(F2));
title('Frequency domain signal via FFT of C(t)');
xlabel('Frequency in Hz')


Understanding the attenuation of a signal I have been told that the signal in the graph below was obtained by passing a seismogram through a high-pass filter to obtain the high frequency-small amplitude waves common in tremors far from the epicenter of an earthquake.

The high-pass filter used to unmask the tremor was selected to avoid the 'attenuation' of the tremor signal.  What exactly is meant by 'attenuation' in this context?

Attenuation in this context refers to a reduction in amplitude of some signal component as it passes through a system (in your case, the system being the highpass filter. Linear filters are designed to pass certain frequency components while stopping (or attenuating) others. As the name implies, a highpass filter is designed to pass high-frequency content while attenuating low-frequency content.

As you noted, the tremor features that you're interested in are high in frequency and small in amplitude. Thus, they may be masked by low-frequency components in the original seismogram. A highpass filter can be used to accentuate the high-frequency tremor waves by attenuating any low-frequency waves that might also be present.


Speech Activity Detection- Adhoc Problem with multi source sounds(birds,trees,heavy air) and how to eliminate them completely Here i will first try to put my application agenda

  Typically i am sampling at 8Khz rate and filtering speech components using a BPF(digital)
  now i have got the 0-4Khz components in which essentially <250Hz are not speech so i will eliminate them also
  now i am trying to find out entropy of frame and trying to mark respective frames as voice or non voice if they belong to band 250<=f<=3750 Hz

so now the problem is while i was doing this at my realtime application side there also non human sounds like birds sound or a ringing of a bell which again fall in to the frequency band i am flagging as speech

so this causes a false alarm of speech even in case of bird sound and bell ringing how to really over come this problem

what else i can try to do for effectvie speech activity detection

This is not aliasing this is not any phenomenon, there are number of birds whose vocal fall in speech range even a frog vocal falls in speech range so essentially a training based speech detection is required to classify whether a given frame is speech or a non speech

one such implementation is shown here

Speech Discrimination based on multiscale spectro temporal modulations


Distort voice to hide talker identity as much as possible I need to develop an app which makes speech recognition. Because I need to be able to constantly improve the recognition capability, I would like to train the engine with real recordings and match the recognized text with the expected text.

To generate the expected text, I would need to listen to the recording, understand what the talker says and manually write the "expected text". To preserve the talker's identity for privacy reasons, I would make available for listening a distorted version of the original audio file, while using the original one for real recognition to improve the system.

From the signal processing point of view, what would be a good flow for this task?

To preserve the talker's identity for privacy reasons, I would make available for listening a distorted version of the original audio file, while using the original one for real recognition to improve the system.

This won't work.

The point is that the whole job of your speech recognition algorithm is to figure out which "components" in your audio signal are of significance, and what they mean, and to ignore things that are irrelevant to the textual content, namely identity of the speaker.

By wanting an algorithm that removes a dominant aspect of the signal but preserves the original speech content, you'd first need to have an algorithm able to tell the content-wise important components from the unimportant. That would be the major part of your speech recognition algorithm, which is exactly what you wanted to build in the first place.

Generally, this sounds like you want to do something that other companies literally spent several millions on with a small team and little experience. I'm afraid that's not really a promising approach.


Is magnitude information enough to reconstruct an audio signal I have used an MMSE STSA estimator to obtain the magnitude of an audio signal. The original signal is combined with white noise and I used an algorithm given in an old research paper by Ephraim and Malah to estimate the amplitude of the clean signal. Now I just have the magnitude information and I don't have the phase. The question is: is this enough to reconstruct an audio signal using ifft. If not can I use any arbitrary phase like the white noise phase? Please let me know

Many thanks in advance

Here is a MATLAB experiment for phase:

[x Fs Nbit]= wavread('test.wav',1024*100); % Get a piece of song ;)
x = x(:,1); % Get a single channel
Y = abs(fft(x)); % Take FFT of the entire piece at once!
y = real(ifft(Y)); % Throw away all-phase information and reconstruct y

y2 = zeros(1,length(x)); % NOW: we will process in BLOCKS!
for i=1:1024:100*1024 % get a block, take its fft, reconstruct a block from fft magnitude
 y2(i:i+1023) = real(ifft(abs(fft(x(i:i+1023)))));
end

figure,plot(x); % piece of a song
figure,plot(y); % obtained from IFFT of MAGNITUDE only of complete piece
figure,plot(y2);% obtained from IFFT of MAGNITUDE only of 1024 sample BLOCKS

sound(x,Fs,Nbit); % the original
sound(y,Fs,Nbit); % from single piece, this is garbage!
sound(y2,Fs,Nbit); % block based, robotic but inteligible!


Why can the :$\mathbf W^*_{MMSE} $ be written as $E[\mathbf s \mathbf y^*]E[\mathbf y \mathbf y^*]^{-1}$ I read this paper: https://arxiv.org/pdf/1305.2460.pdf

And there is a formula about MMSE:$\mathbf W^*_{MMSE}= E[\mathbf s \mathbf y^*]E[\mathbf y \mathbf y^*]^{-1}$

$\mathbf W _{BB}$:Baseband combiner

$\mathbf W _{RF}$:RF combiner

$\mathbf W $:The combination of $\mathbf W_{RF}$ and $\mathbf W _{BB}$

$\mathbf W^*$ : the hermitian  of $\mathbf W$

$\mathbf s$:the signal on the left

$\mathbf y$:the signal when $\mathbf s$ goes  through the channel

Why can the :$\mathbf W^*_{MMSE} $ be written as $E[\mathbf s \mathbf y^*]E[\mathbf y \mathbf y^*]^{-1}$.why should we multiply the these two signal ,$\mathbf s$ and $ \mathbf y$,and average them ?Does anyone have information about this formula ?

The mean squared error equalizer is one that minimizes the mean squred error

$$J(\mathbf{W})=\mathbb{E}\left[\|\mathbf{s}-\mathbf{W}\mathbf{y}\|^2\right]$$

which is the equalizer that minimizes the Euclidean distance between the transmitted signal $\mathbf{s}$ and the equalized received signal $\mathbf{Wy}$. To find the value of $\mathbf{W}$ the minimizes $J(\mathbf{W})$, take the derivative with respect to $\mathbf{W}$ and put it equal to $\mathbf{0}$, which will result in

$$-2\mathbb{E}\left[\left(\mathbf{s}-\mathbf{Wy}\right)\mathbf{y}^H\right]=\mathbf{0}$$
Solve the above equation with respect to $\mathbf{W}$ you get the MMSE equalizer.

Note: $(.)^H$ is the conjugate transpose, which is equivalent to $(.)^*$ in the paper.


Guitar Effects/Filters in real time using low cost DSP? I have an c5502 texas low cost DSP.
I can capture sound trough its P2 input and process every separated sample after sending it to the P2 output.
Im trying to implement some kind of effect in the sound, but this DSP has very low memory so i cant storage the signal to process it as a whole after.

There is any way of implementing effects like overdrive or even low pass filters in real time? I have tried somethings and it seems they arent enough to change the sound by modifying only sample by sample (using past results and samples)...

I'm sure you can implement a simple fuzz or overdrive effect on that board. Apply asymmetrical clipping to you your data, and filter the result with a biquad IIR filter. You'll need a lowpass filter to smooth out some of the nasty high frequencies after clipping (try a cut-off frequency $f_c\in [4..7]$kHz), and you might want a parametric EQ to boost some mid frequencies. You can compute the filter coefficients using the Audio EQ Cookbook. You can probably do much more than that but this should get you started.


How to learn image processing / computer vision in a practical way? I am new to the field of image processing / computer vision and I find it very hard to get my hands on the subject.
I read a lot about the basics/theory but I feel in no way able to use this knowledge in a practical way.

Do you have any tips on how to translate my theoretical knowledge into practical knowledge?
Are there any "cookbooks" for image processing that are not bound to a software/library? (I plan to use Halcon but only find resources for Photoshop or OpenCV.)

Thanks in advance!

Digital Image Processing. Rafael C. Gonzalez, Richard E. Woods. - it is the classical book on Image Processing, that is not bound to any library.

To translate theoretical knowledge into practical knowledge you can start with the book Mastering OpenCV with Practical Computer Vision Projects, which can be downloaded here.


FFT of BPSK signal interpretation I am plotting spectrum of BPSK signal just series of -1 and 1, and I can see random peaks (looks to me to be random) at various locations but I am really confused if I do not multiply my signal by carrier shouldn't the peak of bpsk modulated signal be at DC? If not then how to identify the exact location that peaks should be ?
Here is what I have done:

bits=randi([0 1],1,1000);
bpsk_mod=2*bits-1;
fft_bpsk=fft(bpsk_mod);
len=length(bpsk_mod);
f=[-len/2:len/2-1];
plot(f,abs(fftshift(fft(bpsk_mod)))

A random sequence of +1 and -1 is not really a BPSK signal; it's just a sequence of bits. As such, its FFT has not much meaning, at least not in the sense you're thinking of.

To create an actual BPSK signal, first denote your sequence of $N-1$ bits by $a_k$, where $k=0,\ldots,N-1$ is an integer. Second, choose a bit rate $R_p=1/T_p$. Then, choose a pulse shape, $p(t)$. The pulse shape is typically rectangular, a sinc, a raised cosine, or similar. The important thing is that the pulse is orthogonal to $T_p$-shifted versions of itself; that is, $$\int_{-\infty}^\infty p(t-mT_p)p(t-nT_p)=0$$ for all integers $m\neq n$. Now you're ready to create your BPSK signal, which is given by $$s(t)=\sum_{k=0}^{N-1} a_kp(t-kT_p).$$

Now you can find the FT of $s(t)$ (or the DFT of a properly sampled version of $s(t)$), and you can find the bandwidth and other spectral properties of your BPSK signal.


How to compute the cross correlation of two fuctions I am trying to implement Matlab's cross correlation function but I can't seem to replicate the results obtained when using Matlab's built in function "xcorr".

Documentation for xcorr can be found HERE.

Looking at the Matlab documentation it seems as though the cross correlation is implemented as

$$
  R_{xy}[m] =  \left.\begin{cases} \sum_{n=0}^{N - m - 1} x[n + m] y[n]
& m \geq 0\\ \\ R_{yx}[-m] & m < 0\\ \end{cases} \right\}
$$

Or equivalently as

$$
  R_{xy}[m] =  \left.\begin{cases} \sum_{n=1}^{N - m} x[n + m] y[n]
& m \geq 0\\ \\ \sum_{n=1}^{N + m} y[n - m] x[n] & m < 0\\ \end{cases} \right\}
$$

if x and y have indexes starting at 1.The result will then be C, where

$$
  C(m) = R_{xy}(m - N) \hspace{1cm} \text{for} \hspace{1cm} m = 1, 2,...,2N - 1.
$$

The implementation of this formula in Matlab is shown below:

function [result] = myXCorr(A , B)
%Implementation of the Cross Coorelation Function

N = size(A,2);
M = size(B,2);
result = zeros(1, N + M - 1 );
len = size(result,2);

for m = 1 : len
   arg = (m - N);

   if(arg < 0)
       negativeCondition = 1;
       limit = N + arg;
   else
       negativeCondition = 0;
       limit = N - arg;
   end

   for n = 1:limit
           if(negativeCondition == 0)
            result(m) = result(m) + A(arg + n) + B(n);
           else
            result(m) = result(m) + A(n) + B(n - arg);
           end
   end
end

end

This implementation does not produce the same results as Matlab and I can't seem to find out  what the error is.The plots of the autocorrelation's are shown below.

A = 1:1:10;
B = 1:1:10;

myResult = myXCorr(A,B);
mResult = xcorr(A,B);

lag = -9:1:9;

plot(lag, myResult);
plot(lag, mResult);

Autocorrelation using my function:

Autocorrelation using "xcorr":

Does anyone know what 'm doing wrong here? Any help will be appreciated.

In the loop where you compute your result vector, you should be multiplying the A and B values instead of adding them up.


What do I lose when I set overlap rate to a really high value in STFT If I don't care about computational cost, is there any cons when I put overlap rate really high in STFT?

If you don’t care about computational costs, you can start a window at each sample (e.g. 100% - 1 sample overlap).  It’s well into diminishing returns, but phase vocoder estimation methods work slightly better with greater overlap.


Fourier transform of a tilted line function Assume a line function (line segment to ensure integration):

$$y = a\cdot x + b$$

What is the Fourier transform of the line segment? Intuitively, it may give a sinc function, but actually not.

The Fourier transform of polynomial functions exists in the domain of generalized functions, or distributions. The (non-unitary) Fourier transform of

$$f(t)=t^n\tag{1}$$

is given by (cf. entry 308 in this table)

$$F(\omega)=2\pi j^n\delta^{(n)}(\omega)\tag{2}$$

where $\delta^{(n)}$ is the $n^{th}$ generalized derivative of the Dirac delta impulse.

Consequently, the Fourier transform of

$$f(t)=at+b\tag{3}$$

is given by

$$F(\omega)=2\pi \left[ j\cdot a\, \delta'(\omega)+b\,\delta(\omega)\right]\tag{4}$$

If you want the Fourier transform of a line segment, i.e., of a line multiplied by a rectangular function, you need to convolve $(4)$ with an appropriately chosen sinc function.


The sound of modulated audio signal As I was learned in signal processing, modulation took the signal spectrum to higher frequency range, so I made an experiment for audio file. The audio file "handel.ogg" plays a small part of music "hallelujah".
I was expected that I would hear higher frequency or higher tone of the modulated sound than the unmodulated one when I played the modulated sound, but it didn't.

Why the thing that frequencies in spectrum get higher doesn't mean that the tone or the frequency of signal will get higher? Can someone tell me what's difference, thank you.

Here is my code:

clear;
clc;

% handel.ogg is a sample audio file in matlab.
[x, fs]=audioread('handel.ogg');

% carrier frequency
fc = 100;

% modulated signal
time = (1:numel(x))/fs;
x_modulated = x.*exp(1j*2*pi*fc*time');

% play the music
sound(real(x_modulated), fs)

Eventhough your expectation is sensible on the surface, once you think deeper into audio signals, you'll see that modulation of audio in that way will not produce the expected sonic result. Furthermore your implementation is also wrong.

Consider the audio spectrum; it has two conjugate symmetric lobes into the positive and negative frequencies each about $W$ kHz denoting the bandwidth of your audio. And lower frequency band between $-20$ Hz and $20$ Hz is empty.

First note that you should keep this symmetry in the frequency domain if you want to listen to the modulated bare audio. Your multiplication of a single complex exponential completely destroys the spectral symmetry. [edit: as @DanBoschen reminded, you take the real part of the multiplication so that spectral symmetry is preserved!]. Nevertheless, it's cheaper to multiply with a cosine than a complex exponential, so you better multiply with a cosine() waveform instead of multiplying with a single complex exponential and taking real part of it.

Second, even if  the audio is multipleid with a cosine waveform, a shift of $100$ Hz simply creates a huge overlap in the spectrum resulting in aliasing distortion; i.e., shifting of the spectrum to right and left by $100$ Hz creates right and left lobes to cosite on each other, creating terrible aliasing. To prevent this, you can, in principle, apply separate shifts of left and right single side lobes independently and then sum the result.

Third, even if you do this perfectly, you will not get what you expect as the harmonicity of the audio spectrum is also distorted. The fundamental frequencies and their corresponding higher harmonics do not align on proper places. You should have scaled the frequency axis instead, but you have only linearly shifted the spectrum. To numerically see the effect, consider three harmonic frequencies in the original spectrum at, say, fundemantal at 120, and hamronics 240, and 480 Hzs. When you shift them by, say, 1499 Hz, the new frequencies will be 1619, 1739, and 1979 Hzs. And this new set is not harmonic anymore; i.e., they won't sound natural...

If you want to listen to the same audio at a higher pitch, then you can play with the output DAC sampling frequency. It will strecth the spectrum correctly keeping the harmonicity intact and prevent aliasing also, but it will change the duration of the playback...


Is the Dirac delta (impulse) signal a power signal or an energy signal? I'm a beginner so sorry if this question is very fundamental.  Dirac impulse has finite area i.e = 1. But I've heard that $|\delta(t)|^2$ is undefined. So area under $|\delta(t)|^2$ is also undefined and signal doesn't exist in all time $t$ so it cant be a power signal. So my guess is Neither Power nor Energy signal. Am I right?

[Added a reference on Schwartz's impossibility theorem for products of distribution]

The continuous Dirac delta  $\delta$ is not considered a true function or signal, but a distribution. From its wikipedia page:

  The delta function can also be defined in the sense of distributions
  exactly as above in the one-dimensional case.[25] However, despite
  widespread use in engineering contexts, (2) should be manipulated with
  care, since the product of distributions can only be defined under
  quite narrow circumstances.

It can be defined such that, for any function $f$ satisfying some important properties, and for $a\in \mathbb{R}$:

$$ \int f(t)\delta(t-a)dt=f(a).$$

From the best of my knowledge, those important properties are not satisfied by $\delta$, so one cannot directly replaced $f$ by $\delta$ and get a meaningful result. As far as know, the product of two Dirac distributions is not well defined, unless one talks about $n$-dimensional versions, or so-called "formal" manipulations as used in physics for instance, or more complicated maths.
A short account of Simplified production of Dirac delta function identities is provided by Nicholas Wheeler. If one wants to dig deeper, I'd suggest The Colombeau theory of generalized functions, by Ta Ngoc Tri, 2005:

  Soon after the introduction of his own theory, L. Schwartz published a
  paper in which he showed an impossibility result (see [Sch54]) about
  the product of two arbitrary distributions.

One result is  Schwartz impossibility result. It (somehow) says that if one wants to encompass the derivative of continuously differentiable functions whle keeping  Leibniz's rule of derivation, then one get $\delta^2(|x|)=0$.

However, from an informal point of view, sometimes used in DSP (and in physics), this "product" is, as far as I know, neither energy nor power.
From a logical point of view though, if it does not exist, one could affect this "product" a lot of properties...

Some related posts:

What is product of delta function with itself?
What is the square root of the Dirac Delta Function?
Products and compositions with the Dirac delta function, C. K. Raju, Journal of Physics A: Mathematical and General, Volume 15, Number 2


Rayleigh Fading Simulation I have a question about simulating Rayleigh Fading channel.

I know theory behind the Rayleigh channel and it is sum of two Gaussian random variable where the amplitude is Rayleigh distributed and the phase is uniform. Also if the symbol time is greater than the delay spread then it is flat fading and symbol time is less than delay spread then it is frequency selective channel.

However I got confused after reviewing different posts in this site and other websites in the internet about actually simulating the flat fading and frequency selective channel and I will appreciate any help.

Lets assume I have one tap filter then

N= 1000;
h=(1/sqrt(2))*(randn(N,1)+1i*randn(N,1));

will generate Rayleigh fading channel.

M=4
N= 1000;
h=(1/sqrt(2))*(randn(N,M)+1i*randn(N,M));

Question 1

Does this create a frequency selective channel solely because I have more taps?

Question 2

In non-of these definitions we are taking into account the delay and magnitude of each path and this means a little strange way to create Rayleigh fading channel. Assume for example I have A=[1 0.5 0.5 1] for gain
and T=[2 3  4 5] for delay. Then how may I incorporate this into my simulation? Above definition does not seem to be working.

I will appreciate any help given regarding the correct way to simulate the channel.

Your approach will not model the actual channel but can be used to create the received waveform as having passed through a Rayleigh fading channel. I describe considerations to actually create the received waveform using your approach in either flat and frequency selective conditions below, which should give insight into why it would not be a channel simulation.

To simulate a channel

The Rayleigh channel can be simulated with an FIR filter with a significant number of taps that match the sample rate of the waveform and enough taps that the span of the filter exceeds the delay spread of the channel.  As for "correct ways" of simulating a Rayleigh fading channel itself as a FIR filter, please refer to "Jakes' Channel Model" which uses a "sum of sinusoids" approach to modeling the Rayleigh channel. More information on Jakes' model and other similar approaches is available here.

Details on how to proceed with the "Dictionary Approach"

The dictionary would represent the effect of a "flat fading" Rayleigh channel as a "1 tap filter" as you described if your dictionary was generated at your symbol rate and you multiply (not filter!) each transmitted symbol with a sample from your dictionary. In this fashion each symbol would have an independent amplitude and phase from the Rayleigh channel and that fade would be constant over a symbol duration consistent with flat fading. However the fade would be completely independent from one symbol to the next, so in this case would poorly represent an actual channel.

One solution to that problem is to perform a moving average on the generated dictionary to create a new dictionary with some coherence from symbol to symbol. Much better would be to create an upsampled dictionary by $M$ and then do a moving average over $M$ samples. The higher $M$ is, the more accurately you will represent a true Rayleigh channel.  In this approach the transmitted symbols would also need to be interpolated to match the rate in the dictionary. Your transmitted symbols would be interpolated for pulse shaping and other purposes in your modelling, so the final dictionary using this last approach described can be decimated at the output of the moving average filter to match the number of samples required in your waveform for all other reasons to accurately represent your model.

This approach will also give you the "knob" you are looking for to change between frequency selective and flat fading; what is described above by creating a symbol length moving average over $M$ samples, with $M$ samples per symbol, is at the threshold position between flat fading and frequency selective fading.  If you decrease the span of the moving average to be less than the symbol rate, it will become frequency selective fading, as you will have uncorrelated fading samples within your symbol period which is frequency selective fading.

So in summary to create a received waveform that has passed through a Raylegih channel do the following:

1) Create your transmitted waveform (from $N$ complex baseband symbols) that is sufficiently upsampled to model your channel bandwidth of interest (usually for pulse shaping and spectral regrowth or whatever else you may be modelling; this will typically be between 2 and 10 samples per symbol). I will use $K$ to be the total number of samples, and lower case $m$ to be the upsample factor, so then the total number of samples will be $K = Nm$.

2) Create a vector of $KM$ samples of independent complex Gaussian values using the process you described, where $M$ is an additional smoothing factor.

3) Do a moving average over $L$ samples on the vector created in step 2 above.  Choose $L$ based on modeling either flat fading or frequency selective fading where the threshold boundary between the two cases is when the moving average span is equal to the symbol duration ($L \ge mM$ is flat fading and $L < mM$ is frequency selective fading).

4) Decimate the $KM$ samples a the output of the moving average smoothing filter to match the sample rate of your $K$ interpolated transmitted waveform to create the final dictionary. (Decimate by $M$).

5) Finally multiply each of the $K$ samples in the transmitted waveform by each of the $K$ samples in the final dictionary as a sample by sample multiplication to create the received waveform as having gone through a Rayleigh fading channel.

The higher $M$ is, the more accurate the channel result would be but I would think a factor between 6 to 10 would both be easy to manage and be sufficient for creating the coherence from sample to sample. For more frequency selectivity, an even higher $M$ may be needed using this approach. Ultimately with high frequency selectivity each sample is decorrelated and the use of the moving average is unnecessary (the degree of frequency selectivity would be controlled by the interpolation factor of your waveform). I do like the use of the moving average as it gives you the "knob" to control the degree of fading between flat and frequency selective.

To be most accurate (although a lot more complicated), you could optionally vary L over the course of the moving average process within the range consistent that is consistent with your desired fading environment (flat or frequency selective)...this changes the coherence itself over time which would more accurately represent channel conditions, but that may be a second order effect; meaning not critical to your model).


Calculating FIR filter coefficients with the windowing method and Convolution algorithm Please advice how to calculate FIR filter coefficients with the windowing method using the Rectangular Window function and convolution method? What utility or FIR design application can perform this task?

draw your desired frequency response.  uniformly sample that frequency response $\frac{N}{2}$ times from DC up to just below Nyquist.  reflect those $\frac{N}{2}-1$ samples (not the DC values) as complex conjugate in the upper half of the frequency response above Nyquist.  do the inverse FFT.  if it's rectangular window, you're done.


What is the DFT of DFT of discrete signal What is the discrete fourier transform of the discrete fourier transform of any discrete time signal. Is result same signal? How?

let

$$\begin{align}
 X[k] &= \mathcal{DFT} \Big\{ x[n] \Big\} \\
      &\triangleq \sum\limits_{n=0}^{N-1} x[n] \, e^{-j2\pi nk/N}
\end{align} $$

and

$$ y[n] \triangleq X[n]  $$

(note the substitution of $n$ in for $k$.)  then

$$ Y[k] = \mathcal{DFT} \Big\{ y[n] \Big\} $$

then, if the DFT is defined the most common way (as above):

$$ Y[n] = N \cdot x[-n] $$

where periodicity is implied:  $x[n+N]=x[n]$ for all $n$.


How to get coefficients for sine/cosine function from complex FFT? I'm working on a control system that measures the movement of a vibrating robot arm. Because there is some deadtime, I need to look into the future of the somewhat noisy signal. My idea was to use the frequencies in the sampled signal and produce a fourier function that could be used for extrapolation.

My question: I already have the FFT of the signal vector (containing 60-100 values e.g.) and can see the main frequencies in the amplitude spectrum. Now I want to have a function f(t) which fits to the signal, removes some noise, and can be used to predict the near future of the signal. How do I calculate the coefficients for the sine/cosine functions out of the complex FFT data?

This is a no brainer for anyone who's studied signal processing (Im wondering why its been so long since its been answered)

The Fourier Transform itself is a function that tells you what sinusoid do you need to add together to create this continuous signal. (Its like the Fourier series for more complex signals)

the FFT algorithm gives you COMPLEX coefficients (a+bj) (where j is the imaginary number... it looks like and i but backwards....)

calculate the magnitude of each coefficient A(n) = sqrt(a(n)^2 + b(n)^2)
calculate the phase of each coefficient phi(n) = inverse_tan(b(n)/a(n)); (remember to use atan2 so that it results in the correct quadrant)

now remember based on your sampling rate each coeff represents a "bin" of frequencies, however you should think of them as evenly distributed frequencies:

So if you have N bins each frequency_bin they represent is f = Fs/(2*N)

then your signal becomes

$$
f(x) = \Sigma^N_0 A(n)*cos(2*\pi*f(n) + \phi(n))
$$

also: Predicting function: Impossible. that would be asking for a non-causal system.


low pass filter in fft I have signal. I do Fast Fourier transform on it and now I have this output.

I want to delete the other except two that have high altitude. Can you help me ?

My f sample is 2500.

My complete time is 10s.

My time rate is 0.0004s.

I need to see the code behind your graph because the numbers you put down are not making sense with that figure. Your graph should cover X values from 1-12501 if it's just FFT indices or 0-1250 for actual frequencies.

If you're looking to filter the two frequencies with the greatest magnitude no matter what the input signal is, you can find their frequencies by searching for local maxima. If you have the Signal Processing Toolbox, you can use this code to find the peaks and their indices, and then convert those indices into an actual radial frequency.

[pks,locs] = findpeaks(mag_fft_out);   % this fft should be 0-fs/2 ONLY
[throw_away,I1] = max(pks);            % find biggest peak and return index in peak vector
n1 = locs(I1);                         % find location of peak index in FFT vector
omega1 = 2.*pi.*n1./N;                 % N is the length of your original sequence in samples (T./fs)
pks(I1) = [];                          % remove primary peak from peak vector
locs(I1) = [];                         % remove primary peak from peak location vector
[throw_away,I2] = max(pks);            % find (remaining/second) biggest peak and return index
n2 = locs(I2);                         % find location of peak index in FFT vector
omega2 = 2.*pi.*n2./N;

This will also work for finding the peaks once or you could figure it out manually (Mk I Eyeball + Math)

After that you need to design a filter for each peak frequency. I highly recommend you use the MATLAB 1-D Digital Filter function filter(b,a,x) and choose your coefficients based on what characteristics you need. In all likelihood, a simple 2nd order bandpass will work for both. Look up how to design a digital filter. The article on Wikipedia is a great start and even has diagrams showing how the coefficients work.

I'm not sure this is what you're asking for. You say "delete" but it's not quite that simple. I supposed you could set all other values to zero except the peaks, and then if you want it back as a signal you inverse FFT it, but I've never actually tried that.


Basis functions in DCT I'm new to Discrete Cosine Transform (DCT) and I have a question relating to basis functions.
In DCT, basis functions are defined by:

$$\alpha_p\alpha_q\cos\frac{\pi
\left(2m+1\right)p}{2M}\cos\frac{\pi
\left(2n+1\right)q}{2N},\quad \begin{align}0&\leq p\leq M-1\\0&\leq q\leq N-1\end{align}$$

If $M = N = 64$, we'll have $64$ basis functions. These functions have two variables including $p$ and $q$. Now, to verify, I take the first basis function with $m = n = 0$. Clearly, this function changes with different values of $p$ and $q$. Why it's drawn with same color in the following image (the upper left gray square)?

Because the upper left square would correspond to values of $p=0$ and $q=0$, and $cos(\alpha) = 1 $ if $\alpha=0$, so you get a constant term, that is why it is all gray.
Basically as you increment $p$ and $q$ you get basis images that oscillate more either in vertical or horizontal directions.


variance in the time domain versus variance in frequency domain Hi All: I'm trying to better understand the connection between variance of a time series and the integral of the spectral density over all frequencies. Rather than going through all of the relations, if one looks at the link below,

https://en.wikipedia.org/wiki/Spectral_density, the relations are shown there in the power spectral density section.

To summarize what is said there:

The power, $P$, is defined as the limit as $P=\lim_{T \to \infty}
\frac{1}{2T} \int_{-T}^{T} |x(t)|^2$.

$S_{xx}(\omega)$ is defined as the limit as $S_{xx}(\omega)=\lim_{T \to \infty} E|\hat{x}_T(\omega)|^2$. So, I think of it as the expected instantaneous power at frequency $\omega$. I hope that's the way to think of it.

$\hat{x}_{T}(\omega)$ is defined  as $\frac{1}{\sqrt T} \int_0^{T} x(t) e^{-i\omega t} dt$. So the fourier transform of $x_t$ divided by  $\sqrt T$.

Finally,  $\gamma_{0}$ is defined as the power contained in the total frequency band so $\int_{\omega} S_{xx} (\omega)$. My confusion is the following:

In the time domain-econometrics world, $\gamma_{0}$ is defined as  $\mathrm{Var}(X_{t}) = E(X_{t} - \mu)^2$ where $\mu$ is the mean of $X_{t}$.

So, my confusion stems from the fact that when the variance is defined as power, there is no mean subtracted from $X_{t}$. Does that mean that $X_{t}$ has already been de-meaned when it's written as $X_t$. If not, then how can $\gamma_{0}$ be viewed as the variance when statistical variance is the (expected value of the deviation of a random variable from it's mean) squared yet the view of variance in the DSP world is a weighted (by frequency) combination of $E(X_{t}^2)$ over all frequencies.  My only conclusion is that the mean is already being subtracted out or is assumed to be zero ? I basically don't understand how variance is defined in the DSP framework. Is it viewed as expected power rather than expected variability ? Thanks.

Variance is never defined as power. For a wide-sense stationary random process $X(t)$ with zero mean

$$\mu_X=E\{X(t)\}=0\tag{1}$$

the variance of $X(t)$ equals its power.

The autocorrelation of $X(t)$ is defined by

$$R_X(\tau)=E\{X^*(t)X(t+\tau)\}\tag{2}$$

The power of $X(t)$ is

$$P_X=E\{|X(t)|^2\}=R_X(0)\tag{3}$$

The variance of $X(t)$ is

$$\sigma^2_X=E\{|X(t)-\mu_X|^2\}=E\{|X(t)|^2\}-|\mu_X|^2=R_X(0)-|\mu_X|^2\tag{4}$$

Clearly, if $\mu_X=0$, we have $\sigma^2_X=R_X(0)=P_X$.

Note that the autocovariance of $X(t)$ is defined by

$$C_X(\tau)=E\{[X(t)-\mu_X]^*[X(t+\tau)-\mu_X]\}=R_X(\tau)-|\mu_X|^2\tag{5}$$

Consequently,

$$\sigma^2_X=C_X(0)\tag{6}$$

But since the power spectral density $S_X(\omega)$ is the Fourier transform of $R_X(\tau)$, the integral over $S_X(\omega)$ equals the power of $X(t)$ and not its variance, unless $\mu_X=0$:

$$P_X=R_X(0)=\frac{1}{2\pi}\int_{-\infty}^{\infty}S_X(\omega)d\omega\tag{7}$$

Also take a look at this related answer.

For the difference between the definition of autocorrelation in statistics and in signal processing look here.


Signal analysis using a buffer I am performing an analysis of an output signal in the interval $t<0$ from the operational amplifier given in the following image

where $v_{C_1}(0^-)$ and $v_{C_2}(0^-)$ are the initial conditions to be computed of two capacitors. As you can see, the solution is just below the image.

What I can't understand is why all the measured voltages are $1V$ without knowing the values of the resistances and just knowing that $V_-=V_+$. For sure I am missing something that I don't know. Can someone give me a hint? Thank you.

You should know charaterstics of an opamp it has infinite input resistance and zero output resistance hence opamp doesn't draw any current in the input side hence it is an open circuit and current through is zero apply kvl(Kirchoff voltage law) you will know


Inverse z-transform. Where is mistake? I've already wrote about that trouble (link here), but I don't understand where I've made a mistake.

Full description of the task is as follows:

Z-transform of sequence {x(k)} describe by the equation.
$$
X(z) = \frac{2.5 - 3.15z^{-1} + 1.2z^{-2}}{1 - 2.3z^{-1} + 1.2z^{-2}}
$$

Calculate of samples x(-2), x(-1), x(0), x(1) and x(2) double-side of sequence {x(k)}, which correspond to this z-transform.
Note -
Double-side sequence is sequence to consist samples before zero and after zero.

Step 1

I calculate of the integer part of the expression (just divided the numerator by denominator),

$$
X(z) = \frac{2.5 - 3.15z^{-1} + 1.2z^{-2}}{1 - 2.3z^{-1} + 1.2z^{-2}}\qquad(1.1)
$$
transform to,
$$
X(z) = 2.5 + \frac{2.6z^{-1} - 1.8z^{-2}}{1 - 2.3z^{-1} + 1.2z^{-2}}\qquad(1.2)
$$

Step 2

I transformed denominator like $$\frac{(z – p1)(z – p2)}{z^{-2}}\qquad(2.1)$$.
$$
1 - 2.3z^{-1} + 1.2z^{-2} = (1 - 1.5z^{-1})(1 - 0.8z^{-1})\qquad(2.2)
$$

Step 3

I transformed of the expression from:
$$
X(z) = 2.5 + \frac{2.6z^{-1} - 1.8z^{-2}}{1 - 2.3z^{-1} + 1.2z^{-2}}\qquad(3.1)
$$
to
$$
X(z) = 2.5 + X{_1}(z)\qquad(3.2)
$$
and then $$X{_1}(z)\qquad(3.3)$$ transform to
$$
X{_1}(z) = \frac{A}{1 - 1.5z^{-1}} + \frac{B}{1 - 0.8z^{-1}}\qquad(3.4)
$$

Step 4

Calculate of coefficient A and B

Calculate of A

$$
(1 - 1.5z^{-1})X{_1}(z) = A + \frac{B(1 - 1.5z^{-1})}{1 - 0.8z^{-1}},\qquad(4.1)\qquad if\qquad z=1.5
$$
hence
$$
A = \frac{2.6z^{-1} - 1.8z^{-2}}{1 - 0.8z^{-1}},\qquad(4.2)\qquad if\qquad z = 1.5
$$
so
$$
A = 2
$$

Calculate B

$$
(1 - 0.8z^{-1})X{_1}(z) = \frac{A(1 - 0.8z^{-1})}{1 - 1.5z^{-1}},\qquad(4.3)\qquad if\qquad z=0.8
$$
hence
$$
B = \frac{2.6z^{-1} - 1.8z^{-2}}{1 - 1.5z^{-1}},\qquad(4.4)\qquad if\qquad z = 0.8
$$
so
$$
B = -0.5
$$

Step 5

$$
X{_1}(z) = \frac{2}{1 - 1.5z^{-1}} + \frac{-0.5}{1 - 0.8z^{-1}}\qquad(5.1)
$$
and
$$
X(z) = 2.5 + \frac{2}{1 - 1.5z^{-1}} + \frac{-0.5}{1 - 0.8z^{-1}}\qquad(5.2)
$$
Sequence of {x(k)} has the following representation
$$
x(k) = 2(1.5)^{k}, if k<0;
x(k) = 2.5, if k=0;
x(k) = -0.5(0.8)^{k}, if k>0;
$$
And after that we can calculate Sequence of {x(k)} - x(-2), x( 1), x(0), x(1) и x(2)
 x(-2) = -0.88889;
 x(-1) = -1.33333;
 x(0) = 2.50;
 x(1) = -0.40;
 x(2) = -0.32.

Could you check my solve and say where I made mistake.
Thank you very much!

$$X_1(z) = \frac{2.6z^{-1} - 1.8z^{-2}}{1 - 2.3z^{-1} + 1.2z^{-2}} =\frac{A}{1 - 1.5z^{-1}} + \frac{B}{1 - 0.8z^{-1}} \tag{1}$$

  Mistake is that $A \neq 2$ and $B \neq -0.5$.

Well then, what are they ? There aren't any $A,B$ where the above equation is true. Why? Take the limit $z \rightarrow \infty$. The term $\frac{A}{1 - 1.5z^{-1}} + \frac{B}{1 - 0.8z^{-1}}$ factors to maximally $z^{-1}$ in the numerator whereas $\frac{2.6z^{-1} - 1.8z^{-2}}{1 - 2.3z^{-1} + 1.2z^{-2}} $ contains a $z^{-2}$ term!!!!. So a contradiction !! The correct decomposition is as follows:
$$X_1(z) = \frac{2.6z^{-1} - 1.8z^{-2}}{1 - 2.3z^{-1} + 1.2z^{-2}} =\frac{A}{1 - 1.5z^{-1}} + \frac{B+Cz^{-1}}{1 - 0.8z^{-1}}$$

Taking a common denominator in $\frac{A}{1 - 1.5z^{-1}} + \frac{B}{1 - 0.8z^{-1}}$, equation $(1)$ is now equivalent to
$$2.6z^{-1} - 1.8z^{-2} =A+B + (C -1.5B - 0.8)z^{-1} -1.5Cz^{-1}$$
By identification, we get the following system of equations
$$\begin{bmatrix}
    1.0000 &   1.0000   &      0\\
         0 &  -1.5000  &  1.0000\\
         0  &       0  & -1.5000
\end{bmatrix}
\begin{bmatrix}
A\\B\\C
\end{bmatrix}
=
\begin{bmatrix}
         0\\
    3.4000\\
   -1.8000
\end{bmatrix}
$$
which will give you
\begin{align}
A &= 1.4667\\
B &=  -1.4667 \\
C &= 1.2
\end{align}
Can you continue ?


What does it mean to deconvolve the impulse response I am learning DPS and I came across the problem of deconvolution and removing the impulse response from a signal? This still does little sense to me.

My understanding of the impulse response is to find out the response, usually in frequency, of a device, a filter, a room etc. If I don't mess things here, to me, the impulse response in signal processing is equivalent to numeric 1. Finding the impulse response is convolving a signal with something that is equivalent to 1 that is the impulse response. So if I need to find what the value of x is, I simply multiply x by 1:

$$y = x\times1$$

Therefore if $y = 0.8$, I know $x = 0.8$ as well.

In DSP, we have the equation

$$Y = XH$$

So to know the impulse response of $X$, I multiply $X$ by $H$ and obtain $Y$ which is equal to $X$. I've gone that far.

I don't understand what deconvolution means. I though by applying the impulse response we have already an answer about the frequency of anything we need. What does deconvolution mean here and why would we look for it?

Thanks.

A linear time-invariant (LTI) system is completely characterized by its impulse response $h(t)$. The output signal $y(t)$ given an arbitrary input signal $x(t)$ is given by the convolution of this input signal with the system's impulse response:

$$y(t)=(x\star h)(t)\tag{1}$$

where $\star$ denotes convolution. In the frequency domain, Eq. (1) becomes a simple multiplication:

$$Y(f)=X(f)H(f)$$

where $Y(f)$, $X(f)$, and $H(f)$ are the Fourier transforms of the signals $y(t)$, $x(t)$, and $h(t)$, respectively.

Deconvolution now refers to the process of "undoing" the effect of an LTI system. I.e. you're given the signal $y(t)$ and you want to recover the original signal $x(t)$. E.g., you might want to remove the reverberation added to a signal by a given room. In this case you need to undo the effect of the convolution of the room's impulse response with the original signal.

Naive deconvolution can in principle be done by applying an inverse filter to the signal $y(t)$. The inverse filter ideally has a frequency response $G(f)=1/H(f)$. Note that such a filter does not necessarily exist, e.g. if $H(f)$ has zeros at certain frequencies. Another disadvantage of this naive deconvolution is the effect it has on measurement noise. In frequency regions where $H(f)$ has a small magnitude, the magnitude of $G(f)$ becomes large, and the noise is highly amplified. This is why applying an inverse filter is usually not a good solution. If the system's impulse response is known and if the noise can be modeled as (wide-sense) stationary, Wiener deconvolution can be used. In that case, the original signal is estimated using a mean square error criterion, which takes into account the impulse response as well as the noise. If the impulse response is unknown (or can't be estimated), there are methods available referred to as blind deconvolution.


Affine transformation on a pixel grid I'm implementing affine transformations on images. Currently I'm doing rotation, but general transformation will be added later.

Basically each pixel is mapped by the rotation matrix R:

$\begin{bmatrix} cos(a) & -sin(a) \\ sin(a)  & cos(a) \end{bmatrix} $

The problem is that mapping back does not always come back to a grid point, so I basically rounding it using a brute force approach. What is a better approach (similarly to Bresenham's line drawing algorithm)?

What you need is an interpolation method.

The method you described is called nearest-neighbor, because you pick the pixel that is nearest to the place you actually wanted.

Other methods include:

bi-linear interpolation (Select 4 nearest points, interpolate by x and y according to distance) (see here)
barycentric interpolation (Select 3 nearest points, interpolate in the system of coordinates of the triangle)
bi-cubic interpolation (see here)


Trying to translate formula for blending mode I am using OpenCV C++ for making the blending mode like in photoshop, I want to make overlay mode in it. I search its alternate in OpenCV in which I found this blending way, but its not the overlay as I want to use the overlay method in it.

Overlay method formula from this documentation

(Target > 0.5) * (1 - (1-2*(Target-0.5)) * (1-Blend)) +
(Target <= 0.5) * ((2*Target) * Blend)

Can anyone please explain this formula for implementation in OpenCV C++ , how I can easily understand it for implementation or is there any already build in function for it or any other easy way out.

You can find many ready solution in the Internet. For example here following macro is given:

#define ChannelBlend_Overlay(A,B)    ((uint8)((B < 128) ? (2 * A * B / 255):(255 - 2 * (255 - A) * (255 - B) / 255)))

Additionally you might want to check this package: blImage for source code as you can very easilly translate it. I am sure you will find what you are looking for.


Digital receivers for the reception of non pulse shaped signals We are used to deal with pulse shaped signals, for example signals that go through root raised cosine filters at the transmit, and then  matched with root raised cosine filters at the receive

In addition to these functions there are carrier phase and frequency synchronization, symbol timing synchronization at the output of the matched filter , equalization,  etc..

Very briefly please, how would things be different at the receiver, if any, when receiving a signal that has not been pulse shaped by any conventional pulse shaping filter at the transmit(for whatever reason)?

Since you are talking about digital communication you are going to have a symbols.
When u send a symbol $s(t)$ you are going to send it with a symbol time $T_{sym}$
so we can say that any symbol that is sent equals to $s(t)\Pi(\dfrac{t}{2T_{sym}})$ where $\Pi$ is the rectangular function.
So basically any discrete time signal has a rectangular shaping pulse.

E.g.

In BPSK you will be sending one of the following symbols:
$$
    S_1(t)=1
$$
$$
    S_0(t)=-1
$$
Since I am sending the symbols for a known time $T$

I am practically multiplying them with a rectangular window with the length of $T$

so for the data series $a$ the signal $S(t)$ will be :
$$
   S(t) =  \sum_{k=0}^\infty S_{a_k}(t)\Pi(\dfrac{t-kT_{sym}}{ 2T_{sym} })
$$

And the matched filter in the receiver will be the average over $T_{sym}$ (or integration over $T_{sym}$) which is otherwise known as a rect.

Which is the same way you would do with other shaping pulses


Removing the Multiples of Decrease and Increase in a Spike I have a signal shown below. I got this signal from a glove, where I want to detect the movement of a finger. What I am planning to do is to threshold the signal, so I can get a binary signal indicates the presence of the spike. My problem is that I have multiple increase and decrease within a single spike. So, using threshold technique will produce a "fluctuation" of zeros and ones for each spike. This means that the finger was moved multiple times within that time interval, which is uncorrect.I tried using statistical measures to solve this problem, and I also tried to get the envelop of the signal. However, none of them was able to remove those increases and decreases. What do you suggest to solve this problem

You can use a moving average filter followed by your threshold detector. You can change the span of your moving average filter based on the time length in which you want multiple spikes to be combined.


Understanding of Histogram-of-Oriented-Gradients I am learning about HOG and I understand it from here. A well-explained page with an example. I am not understanding this concept that how it works

  A 16×16 block has 4 histograms which can be concatenated to form a 36
  x 1 element vector and it can be normalized just the way a 3×1 vector
  is normalized.

How this 36*1 came and how we calculated it? and is it compulsory that we always need 9 bin vector? Is it a fixed size for HOG?

came?

How this 36*1 came and how we calculated it?

HOG is an algorithm which:

works on a portion of the image, called "detection window";
divides the "detection window" in a certain number of cells;
associates an histogram (of oriented gradients) to each cell.
Each histogram has N orientation bins. WLOG, here, we'll consider N=9);
"Normalization step" (to make histogram values not affected by lighting
variations): it creates blocks (where each block is made up of 4 cells
(2x2)), and it overlaps them;
creates a vector for each block (so, 4 cells x 9 orientation = 36 components for block)
creates a final vector, associated to the overall "detection window";
the final vector is passed to an SVM classifier.

At each iterations, the "detection window" scans the whole original image, and the HOG algorithm is applied at each scan.

  and is it compulsory that we always need 9 bin vector? Is it a fixed size for HOG?

No, the usage of 9 orientation bins for each histogram is not mandatory; it's only a trade-off to which the authors came through tests, as stated by MimSaad.

In this paper, it is explained that increasing the number of orientation bins improves performance significantly up to about 9 bins; beyond this value, perfomance doesn't improve considerably.
Figure 4b shows this result.


Vector Quantization in Hidden Markov Models In the field of Gesture Recognition and Speech Recognition (using HMM classifiers), do multiple HMM classifiers (one for each class) share common Codebook (i.e common set of symbols) in Vector Quantization?

If that is not the case, how does the test/classification stage work?

To which Codebook the test input vector belongs to (we need this to determine the symbol, in turn which is fed to the HMM for classification)?

Some background info:

One of the most earliest works done in Action Recognition used different set of symbols for each of the HMM classifier.
This journal in speech recognition by Rabiner suggest using common Codebook (with different HMM classifiers) for identifying 10 digits uttered by large number of users (in page 12),

do multiple HMM classifiers (one for each class) share common Codebook
  (i.e common set of symbols) in Vector Quantization?

Yes, for example in large vocabulary speech recognition there is a big codebook of about 4000-10000 gaussians which are shared across 20-30k HMMs

  To which Codebook the test input vector belongs to (we need this to
  determine the symbol, in turn which is fed to the HMM for
  classification)?

In modern ASR systems codebook defines a set of GMM probability distributions, not the discrete symbols. So "belongs" doesn't really have sense, it's more a probability value. But during the classification a test vector is usually evaluated for the subset of all codebooks which is active for this stage.

Active codebooks are tracked across test vectors based on active HMMs and preceeding history. Whole HMM set is pruned in order to get active ones. Usually there are about 10k active HMMs on each step and about 1k codebooks.


Pyramid Match Kernel: How to fit grid around data points? I am planning to use Pyramid Match Kernel method for object recognition from depth images: I am going to extract a set of feature vectors $x \in \mathbb{R}^d$ for each object instance and then I want to use the Pyramid Match Kernel as described in http://jmlr.csail.mit.edu/papers/volume8/grauman07a/grauman07a.pdf in Grauman's paper, in order to use the set of feature vectors for each object with a SVM classifier.

I have a problem to understand how exactly we are supposed to build bins around the data points to begin with. The paper gives the following definition:

What does "inter-vector" distance mean here? I think of the following: We consider all $x_i$ values for each $X$ in $S$. Then we find two vectors with smallest distance $d$ ,(unique ones) and then scale all vectors with $1/d$. Then $D$ is the maximum absolute value of a vector element in $S$.

Is this what the authors meant here? What about if we are going to test a new vector $y$, how should we scale it?

No, that's not it. Look at the last sentence of your quote:

  ... which may be enforced by scaling the data to some precision and
  truncating to integer values.

So this is exactly what you do. Take all your feature vectors, multiply them by 1000 (or some other factor), and truncate to integer. Then the distance between any two unique feature vectors will have to be >= 1. If you have a new vector y, then you multiply it by the same factor and truncate.


polyBLAMP anti-aliasing in C++ How should 4 point polyBLAMP work. It does not really seem to dampen aliasing.

I have an floating point digital audio oscillator that makes a +1.0f jump at the end of every period.

My question is, have I misunderstood how it should work?

I have been using these as reference:
http://www.dafx17.eca.ed.ac.uk/papers/DAFx17_paper_100.pdf
https://www.researchgate.net/publication/307990687_Rounding_Corners_with_BLAMP

Here is how I thought it should be done: (the code gets run every sample)

float t = phase / period; //normalized position in period
float dt = deltaPhase / period; // normalized phasestep

float polyBLAMP = 0.0f;
if (t < dt) {
    float d = (t)/dt; // fractional delay
    polyBLAMP = (3.0f*powf(d,5.0f)-10.0f*powf(d,4.0f)+40.0f*powf(d,2.0f)-60.0f*d/2.0f+28.0f)/120.0f;
} else if (t < dt*2.0f) {
    float d = (t-dt)/dt;
    polyBLAMP = (-powf(d,5.0f)+5.0f*powf(d,4.0f)-10.0f*powf(d,3.0f)+10.0f*powf(d,2.0f)-5.0f*d+1.0f)/120.0f;
} else if (t > 1.0f - dt) {
    float d = (dt-1.0f+t)/dt;
    polyBLAMP = (-3.0f*powf(d,5.0f)+5.0f*powf(d,4.0f)+10.0f*powf(d,3.0f)+10.0f*powf(d,2.0f)+5.0f*d+1.0f)/120.0f;
} else if (t > 1.0f - dt*2.0f) {
    float d = (2.0f*dt-1.0f+t)/dt;
    polyBLAMP = powf(d,5.0f)/120.0f;
}
float output = osc(phase)+polyBLAMP;

First of all, if your phasor oscillator has a step discontinuity (as opposed to a ramp one) at $\phi\bmod 1 = 0$, then you should use the four-point, fourth order polyBLEP residual, not the fifth-order polyBLAMP one.

Second, the easiest way to simplify the piecewise polynomials is to see that the residuals are the successive integrals of the piecewise approximation to the band-limited pulse, or in the terminology of Esqueda, Välimäki & Bilbao [EVB2016], the “third-order B-spline basis function.” This basis function is needlessly complicated as given; it can be expressed as the difference of two overlapping piecewise polynomials:

$$
\frac{1}{6}\,\begin{cases}
(2+d)^3 & \text{if}\ {-2} \le d < 0 \\
(2-d)^3 & \text{if}\ \phantom{-}0 \le d < 2 \\
0 & \text{otherwise}
\end{cases} -\
\frac{1}{6}\,\begin{cases}
4(1+d)^3 & \text{if}\ {-1} \le d < 0 \\
4(1-d)^3 & \text{if}\ \phantom{-}0 \le d < 1 \\
0 & \text{otherwise}
\end{cases}
$$

From this, the successive residuals are trivial integrals of polynomials.

Third, the residuals are centered on $d=0$, where $d$ is the offset from the discontinuity to the immediately previous sample point, and the function extends two samples on either side of $0$. In order to apply them to a function of a phasor, the positive portion of the piecewise definition must be shifted to the $\phi=0$ discontinuity while the negative portion must be shifted to the $\phi=1$ discontinuity. Let $0\le\phi<1$ be the phase and $d\phi = f_0/f_s$ the phase increment. The modified polyBLEP in single-sided form is:

$$
\frac{1}{12}\,\begin{cases}
\left(2-\frac{\phi}{d\phi}\right)^4 & \text{if}\ 0\le \phi < 2d\phi \\
0 & \text{otherwise}
\end{cases} -\
\frac{1}{12}\,\begin{cases}
4\left(1-\frac{\phi}{d\phi}\right)^4 & \text{if}\ 0\le \phi < d\phi \\
0 & \text{otherwise}
\end{cases}
$$

(for the double-sided form, add two further terms under the substitution $\phi := 1 - \phi$.) In Java:

public static double blep(double t, final double dt) {
    double y = 0;
    if (0 <= t && t < 2 * dt) {
        final double x = t / dt;
        double u = 2 - x;
        u *= u;
        u *= u;
        y += u;
        if (t < dt) {
            double v = 1 - x;
            v *= v;
            v *= v;
            y -= 4 * v;
        }
    }
    return y / 12;
}

and this must be added on the right of the discontinuity and subtracted on its left. Hence, the sawtooth is:

public static double saw(double t, double dt) {
    double v = 2 * t - 1;
    v += blep(t, dt);
    v -= blep(1 - t, dt);
    return v;
}

(note that this is a positive-going ramp. For a negative going one, built on $1-2\phi$, subtract and then add the BLEPs.) A square wave or pulse with variable duty cycle $2d\phi < d < 1 - 2d\phi$ can be built from a difference of phase shifted sawtooth waves, or directly by “polyBLEPping” on both discontinuities in the cycle:

public static double sqr(double t, double dt) {
    double v = Math.signum(2 * t - 1);
    v += blep(t, dt);
    v -= blep(1 - t, dt);
    v += blep(0.5 - t, dt);
    v -= blep(t - 0.5, dt);
    return v;
}

If your oscillator is a triangle wave, then you do need to use the polyBLAMP residual. This can be found by integrating the polyBLEP residual, and scaling appropriately to the slope of the ramp it will be applied to. In single-sided form:

$$
-\frac{d\phi}{15}\,\begin{cases}
\left(2-\frac{\phi}{d\phi}\right)^5 & \text{if}\ 0 \le \phi < 2\phi \\
0 & \text{otherwise}
\end{cases} +\
\frac{d\phi}{15}\,\begin{cases}
4\left(1-\frac{\phi}{d\phi}\right)^5 & \text{if}\ 0 \le \phi < d\phi \\
0 & \text{otherwise}
\end{cases}
$$

(for the double-sided form, add two further terms under the substitution $\phi := 1 - \phi$.) The polyBLAMP residual in Java:

public static double blamp(double t, double dt) {
    double y = 0;
    if (0 <= t && t < 2 * dt) {
        final double x = t / dt;
        double u = 2 - x, u2 = u * u;
        u *= u2 * u2;
        y -= u;
        if (t < dt) {
            double v = 1 - x, v2 = v * v;
            v *= v2 * v2;
            y += 4 * v;
        }
    }
    return y * dt / 15;
}

Then the triangle oscillator is:

public static double tri(double t, double dt) {
    double v = 2 * Math.abs(2 * t - 1) - 1;
    v += blamp(t, dt);
    v += blamp(1 - t, dt);
    t += 0.5;
    t -= Math.floor(t);
    v -= blamp(t, dt);
    v -= blamp(1 - t, dt);
    return v;
}


Does this make sense for cross-correlation in frequency domain? If I have a signal A at 2 kHz and signal B mixed at 2 kHz + 3kHz, and then i take the FFT of them both, and then cross-correlate them many times, the 3kHz peak should disappear? Since this peak is not present in signal A.

Thanks

If your intention is to cross-correlate the two signals, they will only disappear (zero result) if you are judicious with your time duration of the signals.

The description and question is somewhat confusing however since "cross-correlation in the frequency domain" typically means using the following relationship to perform (circular) cross-correlation:

$${\tt xcorr}(a,b) = {\tt ifft}\left({\tt fft}(a) {\tt fft}(b)^*\right)$$

where the * represents a complex conjugate. And the intention here is to cross-correlate time domain signals through multiplication in the frequency domain.

To cross correlate the actual frequency waveforms meaning you do want to correlate the frequency spectrums themselves (treated as arbitrary waveforms) is simply then cross correlating impulses, so would have the expected result for doing that: the cross correlation function of impulses results in impulses with their related shifts, in this case shifted in frequency.

To achieve what I think you are trying to accomplish given your final sentence and expected result, you would multiply in the frequency domain (not convolve!), in which case the expected signal would "disappear" with proper processing as detailed below since this would indeed represent the equivalent of convolution in the time domain.

Note that the magnitude of correlation versus a frequency offset is a Sinc function with the first null at $1/T$ where $T$ is the time duration of the signal for which you take the FFT.

For example consider a signal that is a single tone with frequency 1 MHz and correlate it to another single tone with 1 KHz offset (1.001 MHz). If your time duration of the two signals is exactly 1 ms then you will indeed get a zero result (this is the foundation of the “Orthogonal” in OFDM). But if for example the duration was 1.5  ms, you will still see a significant result in the cross correlation of the two!

This ultimately defines the resolution bandwidth of the correlation and because of the relatively high sidelobes that don’t rapidly decease, we often use windowing techniques prior to correlating (and taking FFT’s) to decrease the cross correlation of more distant signals when a high dynamic range is desired (but windowing increases the main lobe—- pushes further out the location of that first correlation null beyond the frequency offset given by $1/T$).


Decomposition of a signal to slow and fast components I am currently recording biological responses which are triggered by different events. The figure below shows the original signal (black) and the occurrence of input events (colored dots). There is typically few seconds of a delay in  response to different events, and only some of the events generate a response. These phasic responses occur simultaneously with slowly changing fluctuations (tonic component), which are also biologically important (rather than measurement noise). So it is important to come up with a model that can decompose these signals with few constraints that are biologically plausible: signal = tonic + phasic.

In this system, phasic responses can be well modeled with typical linear system characterization methods based on impulse response function. Also note that, these phasic responses are always positive, that is biologically speaking there is no inhibition. For this reason, if there were no slow fluctuations in the recordings, the recorded signal would never go below zero, except for measurement noise.

What is more problematic is to account for baseline shifts (the tonic component). I would like to come up with a method that can be used to model the tonic component and phasic responses simultaneously based on the following assumptions:

1/ tonic component changes slowly.
2/ phasic component changes fast.
3/ phasic component is always positive, that is: signal - tonic should contain as negative values as possible.

The problem boils down to extracting slow and fast components of a signal. I have already 4 methods in mind, and I would like to see what you think is best suited or come up with an alternative.

The first thing that comes to mind is to low-pass filter the signal to obtain an estimation of the tonic component (see blue line). The main problem here is that the estimated tonic component violates the above-mentionned constraint #3. According to this assumption, no points in the recorded signal should be below the tonic component. Instead the tonic component should pass through all the data points where the phasic response is close to zero (schematically shown with the red curve in the figure below).

Another approach would be to make an assumption about the duration of phasic response. In doing so, one could use recorded samples that are at least as far as the duration of phasic response as keypoints and interpolate using cubic splines for the remaining data points. This is actually how I have drawn the red curve, but I feel like the interpolation of data points between the key points is problematic, because there are not many key points.

A third approach would be to work on first or second derivative space and characterize the system in this domain, as it is supposed to be less sensitive to low-variations.

A fourth approach consists of extracting slow components. One could use rely more on the history, when the signal is changing fast, and update the history with new data when the signal changes slowly. I am not sure though how to mathematically tackle this one.

Your problem is pervasive in signal and image processing: being able to separate a trend (with some assumptions) from a signal (with some other assumptions) and potentially additional noise. This is an instance of (semi-)informed source separation: a sine from a drift,  a texture from a cartoon image, a biological information from instrumental artifacts, etc, that requires a more involved modeling and optimization framework.

Your choices are related to fundamental questions related to the model you can build based on physical and mathematical assumptions, aside from the algorithmic tools you can program.

1) Low-pass filtering:
Standard low-pass filtering is linear and time-invariant. Linearity is very consistent with "signal = tonic + phasic". But phasic could be (on theory) negative, as "$2+(-1) = 1$". As your graphic does not include x/y scales, "above" is clear, but i am not sure of what you can model with "no points in the recorded signal should be below the tonic component", since additional noise can be fluctuating too. But positivity is a property which can be taken into account, with some non-linearity caveats.

2) Duration:
"that are at least as far as the duration of phasic response" : your assumption seems to me to related to a time-delay, yet i cannot interpret it precisely so far.

3) High-order derivatives:
"A third approach would be to work on first or second derivative space and characterize the system in this domain, as it is supposed to be less sensitive to low-variations" Indeed, slow variations are expected to be with limited derivatives. This could result in a "sparsity" prior assumption.

4) Causality: "rely more on the history" there is some causality assumption behind that. This too is a tough stuff.

All your four approaches seem consistent to me, and may help you ponder and state your assumptions. They are otherwise termed as “morphological component analysis (MCA)”, “geometric separation” or “clustered sparsity”, adding positivity. Such topics are addressed for instance in "Blind source separation of positive and partially correlated data" http://dx.doi.org/10.1016/j.sigpro.2005.03.006.

You may find additional history for instance in the introduction of (mong other papers "Chromatogram baseline estimation and denoising using sparsity (BEADS)"
http://dx.doi.org/10.1016/j.chemolab.2014.09.014
which also deals with positivity, low-pass-filtering and derivative sparisty, yet the causal approach is not considered. So you can specify your model in more precise terms (you did do talk about the noise yet), you may find some tracks to solve your questions


How to calculate LF/HF power ratio from heart rate RR-intervals I've been searching everywhere, but I'm having trouble finding a step-by-step guide for someone with little knowledge of signal processing.

I have a file with a list of RR-intervals in milliseconds. The RR-interval is the time between successive heartbeats, so if you took the cumulative sum, you could think of it as a timestamp for each heart beat.

One method of calculating heart rate variability is to find the ratio of power between the low frequency band (0.04 - 0.15 Hz) and the high frequency band (0.15 - 0.40 Hz) in the heart rate signal.

How can I find these values? The RR-intervals aren't like an evenly sampled signal, so I think resampling would be involved. Thank you very much.

Here is a link with a bunch of stuff relating to HVR and how to calculate various features using python discrete-heart-rate-signal-using-python-part-1. Hope it helps.


Do I need an overlapping in STFT if I want to play the signal back? I want to write an algorithm to implement a pitch-shifter, like in this article: http://www.guitarpitchshifter.com/matlab.html. And I need to do it in realtime, not with a signal already defined earlier. I don't understand if I need overlapping in STFT to play the sound back.

Why can't I just record, say, 1024 samples, then process them and then output it and so on (without overlapping)? This seems to be easier to implement but I guess there is a reason they added an overlapping which I can't understand.

Can you help me understand, is overlapping needed in there? And if yes, why?

Note that the algorithm in the article uses a windowed FFT for analysis.  Window functions help avoid certain FFT artifacts, but they strongly attenuate any signal at the ends of the window.  So if you don't overlap, almost all the audio (any pitch changes, etc.) at both edges of each window would be lost, and you would end up with a choppy re-synthesized result.  Also, the window function usually peaks in the middle of the window, and is thus produces more accurate analysis there.  So the closer together the window centers, the more accurate the analysis involved in the pitch shift algorithm.  But the window needs to be longer than this distance for accuracy in frequency estimation.  Thus the overlap of 75% instead of just 50% in the article.

Synthesis for playback often involves the phase vocoder algorithm plus cross-fading of overlapped windows to avoid discontinuity ticks and other glitches at the synthesis window boundaries.


Geometric explanation of a methodology in the article about Image Denoising In article Ghimpeteanu G., et al. - A Decomposition Framework for Image Denoising Algorithms, I found as below:

  Let $\displaystyle I : \Omega \subset \mathbb{R}^2\mapsto \mathbb{R}$ be a gray-level image, and $(x, y)$ be the standard coordinate system of $\mathbb{R}^2$.We denote by $Ix
 \hspace{.2cm}resp.Iy$ the derivative of $I$ with respect to $x
 \hspace{.2cm}resp. y$, and by $\nabla I$ the gradient of $I$ . Our
  image decomposition model for $I$ is a two-stages approach: first, we
  construct an orthonormal moving frame $(Z_1, Z_2, N)$ of $(\mathbb{R}^3,\|
\|_2)$ over $\Omega$ that encodes the local geometry of $I$ . Then, we
  compute the components $(J^1, J^2, J^3)$ of the $\mathbb{R}^3-$valued function
  $(0, 0, I )$ in that moving frame. More precisely, we consider a
  scaled version $\mu I$ of $I$ , for $\mu \in ]0, 1]$, and its graph,
  which is the surface $S$ in $\mathbb{R}^3$ parametrized by $\psi : (x, y)
\mapsto (x, y,μ I (x, y))$

See snapshot-

Please someone explain (possibly visualize geometrically) 4th and 5th line of the excerpt. i.e

  Our image decomposition model for $I$ is a two-stages approach: first, we construct an orthonormal moving frame $(Z_1, Z_2, N)$ of $(\mathbb{R}^3,\| \|_2)$ over $\Omega$ that encodes the local geometry of $I$ . Then, we compute the components $(J^1, J^2, J^3)$ of the $\mathbb{R}^3-$valued function $(0, 0, I )$ in that moving frame.

The explanation below the section you quote seems to be clear to me:

It's forming a vector field made up of the largest gradient direction, smallest gradient direction, and the normal to them both.


How exactly does 8PSK work on each I and Q leg? I am working on making a modulator in VHDL and am implementing 8PSK. Currently, I have BPSK and QPSK looking good, and I have verified the spectrums. However, now that I am moving to 8PSK, I am stuck on how the data mux to the IQ legs should work?

Currently, my implementation is two-fold for QPSK. My first method splits the bits into the I and Q legs for even/odd bits (I am not pulse shaping, purely PCM), then I basically take this I or Q bit and xor it with the MSB of the signed DCO value, in-phase for even bits and quadrature for odd bits (I am using a digitally controlled oscillator as my carrier). The second method I have is I just generate the symbols for QPSK and feed my DCO with a phase shift at 45, 135, 225, and 315 (this saves logic as I only have to use 1 DCO). I know I can use the phase request for 8PSK as well, however, this method produces a lot of spurs.

If I want to perform the following block diagram for 8PSK, how does this work? Am I accumulating 3 bits, or one symbol, for each leg, and performing a phase shift? I'm really struggling to figure out how to do this. I know I can just accumulate one symbol and send a phase shift to the DCO, but this doesn't use I and Q and feels like a crappy way of doing this, especially considering it doesn't work well when I tried it for QPSK. Here is the fundamental question, when I go from serial to parallel through the data MUX, what exactly is going onto the I and the Q leg? If it is still 1 bit, how is this any different than QPSK?

You are mapping a group of 3 bits (an 8PSK symbol) to one of 8 coordinates in the I-Q plane.  I'm not an FPGA guy, but you can do this with a "switch" statement or a lookup table.  I don't know your particular mapping, but a straightforward mapping from bits to I & Q values would be:

$$(0,0,0) \rightarrow (1, 0)$$
$$(0,0,1) \rightarrow (1/\sqrt{2},1/\sqrt{2})$$
$$(0,1,0) \rightarrow (0, 1)$$
$$(0,1,1) \rightarrow (-1/\sqrt{2},1/\sqrt{2})$$
$$(1,0,0) \rightarrow (-1, 0)$$
$$(1,0,1) \rightarrow (-1/\sqrt{2},-1/\sqrt{2})$$
$$(1,1,0) \rightarrow (0, -1)$$
$$(1,1,1) \rightarrow (1/\sqrt{2},-1/\sqrt{2})$$

Note that, in the above mapping, if you take entries 2, 4, 6, and 8, and ignore the least significant bit, you get a mapping for QPSK.


Normalized Wavelet power spectrum I am analyzing the spectral components of a time series using the continuous wavelet transform following Torrence and Compo (1998). I would like to partition the signal variability or spectral power across different scales (e.g. 50% of the variability is in the scales X to Y).

As pointed out here the global wavelet spectrum is a biased estimator of the power spectrum, so is there a way to plot the power spectrum in a normalized way, where I can compare the power at different frequencies?

Would it be correct to just reconstruct each scale range I am interested in and calculate the variance over time?

I've used the Liu-Liang-Weisberg paper before, and it worked pretty well. There's some example code on this page.

In short, they propose that squaring wavelet coefficients and dividing them by their respective scales normalizes the spectrum. The paper provides proof and motivation.


Why is $\pi$ in the domain of filters? I'm just beginning to learn (FIR) filter design and I was wondering, why is $\pi$ involved in the domain of such filters?

I assume it's related to some trigonometric functions, but why do I see e.g. cutoff frequency being related to $\pi$?

You can see for example here that the variable $\omega \in [-\pi, \pi]$ in function $z=e^{i \omega} $ of the transfer function of a FIR filter:

https://en.wikipedia.org/wiki/Finite_impulse_response#Transfer_function

In discrete time, distinct frequencies only have sense in a range of $2\pi$ (where the units are radians per sample).

A sinusoid (or complex exponential) of frequency $2\pi$ is indistinguishable from one with frequency 0. If in doubt, evaluate $x[n] = \cos(2\pi \times n)$ and $x[n] = \cos(0 \times n)$. Same thing with any two frequencies separated $2\pi$ or multiples of $2\pi$.

For example, you can easily show that $\cos( \theta_0 \times n ) = \cos((\theta_0 + 2k\pi) \times n)$ for any integer $k$.

So a range of size $2\pi$ is enough not to repeat frequencies.

In general, the range chosen is $-\pi \cdots \pi$, where $\pi$ is the maximum frequency ( $\cos(\pi \times n ) = (-1)^n$ ).

In the $z$ domain, the frequency response is the transfer function ($H(z)$) evaluated on the unit circle ($z = e^{j\theta}$), so here again you have the range $2\pi$ before repeating your path.


Proof of transfer function factorization $\frac{b_0}{a_0} \frac{\prod_{k=1}^M (1-c_kz^{-1})}{\prod_{k=1}^N(1-d_kz^{-1})}$ This is from Oppenheim's Discrete-Time Signal Processing, but the book doesn't seem to describe how the factorization is done.

The transfer function:

$$H(z)=\frac{Y(z)}{X(z)} = \frac{\sum_{k=0}^M b_kz^{-k}}{\sum_{k=0}^Na_kz^{-k}}$$

can be factored to the form:

$$\frac{b_0}{a_0} \frac{\prod_{k=1}^M (1-c_kz^{-1})}{\prod_{k=1}^N(1-d_kz^{-1})}$$

where each of the factors $(1-c_kz^{-1})$ contribute a zero at $z=c_k$ and pole at $z=0$. Similarly $(1-d_kz^{-1})$ contributes a zero at $z=0$ and a pole at $z=d_k$.

How is the above factorization done?
Is this the only representation?

It's a mere application of the fundamental theorem of algebra of polynomials in 1D, which states that any polynomial $p(x)$ of $x$ of degree $n$, can be written (expressed) in either increasing (or decreasing) powers of $x$ as in $$p(x)=a_0x^n+a_1x^{n-1}+...+a_{n-1}x+a_n = \sum_{k=0}^{n}{a_kx^{n-k}}$$  or as a product of its $n$ roots in the factors of $(x-r_k)$ as in $$p(x) = K(x-r_1)(x-r_2)...(x-r_n)=K\prod_{k=1}^{n}(x-r_n)$$ where each such factor represents a root of the polynomial $p(x)$.

The $n^{th}$ degree poynomial is assumed to have $n$ roots, either distinct or repeating or in complex and/or complex conjugate forms, depending on the coefficients $a_k$.

Applying this to the $Z$-Transform of the form discussed:
$$X(z)=\sum_{n=0}^{M}{a_nz^{-n}}=a_0+a_1z^{-1}+...+a_Mz^{-M}$$ this can be seen either as a polynomial in the variable $z^{-1}$, or after factoring out $a_0z^{-M}$ ($a_0$ being nonzero) it will be a polynomials in $z$ as given $$X(z)=a_0z^{-M}\left(z^{M}+\frac{a_1}{a_0}z^{M-1}+...+\frac{a_M}{a_0}\right)$$ And the paranthesis can be written in product form as:
$$X(z)= a_0z^{-M}\prod_{k=1}^{M}{(z-r_k)}$$
You can distribute the factor $z^{-M}$ into each product form as in: $$z^{-1}(z-p_m) = (1-p_mz^{-1}) $$ which yields the final required form as:
$$X(z) = a_0 \prod_{k=1}^{M}(1 - r_k z^{-1})$$ where each $r_k$ is a root.

When $X(z)$ is expressed as a rational traform as in $$X(z)= \frac{P(z)}{Q(z)} = \frac{\sum_{k=0}^{M}{b_kz^{-k}}}{\sum_{k=0}^{N}{a_kz^{-k}}} $$ apply this factorization to both of numerator and denominator to get the final form:$$X(z)= \frac{\sum_{k=0}^{M}{b_kz^{-k}}}{\sum_{k=0}^{N}{a_kz^{-k}}}=\frac{b_0}{a_0} \frac{\prod_{k=1}^{M}{(1-c_kz^{-1})}}{\prod_{k=1}^{N}{(1-d_kz^{-1})}}$$


Calculating Delta MFCC, correct formula? I have seen different formulations of how to calculate delta MFCCs.  In general, I want to be consistent in such a way that if I have other features to look at other than MFCCs, the way I would calculate their deltas would be the same.

According to this link, HTK supported different ways of looking at how to calculate MFCCs.  One was using a regression formula.  Another is just doing a simple first order difference.  Another listed on that page is doing a modified kind of simple difference.  Another link also showcases the regression formula.

Which one of these is the best way to go?
In what situations might one outperform the other?  I can see that doing first order differences can be easily done in code and doesn't require end frame replication, but many sources I've been reading for MFCCs at least do not use the first order difference.
For other types of delta features that arent MFCCs, is that also the case?

The way you compute MFCCs and derivatives has little-to-no influence on the final result, so just pick one and go ahead. Although when applied to speech recognition and similar problems the same algorithm should be used for training and decoding. Thus if you use a pre-trained model you have to know how features for its training were computed. If not sure, you can always take a look at real-world examples - CMU Sphinx and Kaldi.


Situations where complexity is too big to exceed linearity and gaussianity I'm studing about Kalman Filter and Particle Filter in multiple target tracking in computer vision (tracking pedestrians). Reading sientific papers I'm colliding with a lot of sentences like:

  The assumptions of linearity and Gaussianity of the Kalman Filter
  cannote manage complex scenario.

Now.. I'd like to find some real example of scenario where those assumptions are exceeded. Because I don't really understand where I should use Kalman Filter or not. Can anyone give me some example?

Real videos (like youtube or similar) are very good example.

Well, let's look at the two issues: 1) linearity and 2) Gaussianity.

Linearity

If you're imaging moving 3D objects (people) with a single camera, then you're working with a 2D projection of those 3D objects.  That dimensionality reduction can cause non-linearities to appear.

Take a 2D to 1D example:  an object moving in a circle in 2D.  The object is moving with constant velocity and constant acceleration.  All nice and linear.  Project that object onto a 1D surface (a line) and the velocity and acceleration are now (co)sinusoidal.

One of the underlying assumptions in Kalman filtering (well, the sort usually presented) is that the system you're modeling is linear.  When dimensionality reduction (3D $\rightarrow$ 2D) is involved, this may not be true.

Gaussianity

When a system is linear, any noises tend to be additive.  When Gaussian noise is added to Gaussian noise (or, for enough samples, noise of just about any distribution is added to noise of just about any distribution), the result is also Gaussian.

If the system is non-linear, then the operations involved cause Gaussian noise to be transformed into potentially non-Gaussian noise.

For example, if you divide one Gaussian variable with another, the distribution of the result is Cauchy.

What's so bad about that?  The mean (average) of a Cauchy random variable does not exist.  The variance of a Cauchy random variable does not exist.

The Kalman filter is all about tracking means and variances of variables in the system.  If the system is non-linear, it's possible that the relevant quantities' means and variances do not exist, so application of a Kalman filter is meaningless.


Difference between two sine generators I have two sine oscillators and I get these two different spectrum analysis (let's disregard the slight frequency difference). Is there any deductions that can be made about the "quality", or sine approximation error or any other type of information about the sine generators derived by these graphs? Why the second is more steep? Is it normal? Which of the two is closer to the ideal sine FFT spectrum? What to read so to understand differences of this type?

The sample rate is 32k and there are 8192 frequency bins.

The code that produces this difference is:

a. (upper graph)

double freq;
static double sine_rate;

sine_rate += freq;
*var = (long) ( sin( sine_rate ) * 2142483647.0 );

b.

double freq;
static double sine_rate;

sine_rate += freq;
*var = (long) ( sin( sine_rate * PI ) * 2142483647.0 );

double "freq" is provided when the function is called and can be any double between 0.0 and 1.0. The difference is the multiplication by 3.14.

The two oscillator types are equivalent, because the factor PI can be moved to freq and need not be present in the final equation.

The difference in the plots is not due to one oscillator type being better than the other; they are equivalent. The difference seen is spectral analysis artifacts that depend on the frequency of the sine and how that relates to the Discrete Fourier transform (DFT) length. DFT decomposes the time-domain signal into harmonic frequencies (multiples of a fundamental frequency that has a fundamental period equal to the analysis window length). If the sine frequency is not exactly at such a harmonic frequency, the sine "energy" still needs to appear somewhere in the frequency domain and will be spread to multiple frequency bins. The plots in this question on such spectral leakage are very similar to yours. Spectral leakage can be ameliorated by multiplying the data by a suitable window function before DFT.

With the oscillator of type a, if freq is of form some_integer*2*PI/N, where N is the DFT length, then the sine frequency is one of the harmonic frequencies seen by DFT and there will be no spectral leakage.