Analysis of the ratio of two moving average filters I have created a basic segmentation algorithm for 1-D signals (e.g. audio) based on a ratio of energy averages. The classification rule involves comparing the signal energy measured over a short time window divided by that over a longer window, with a fixed threshold.
More precisely, if $x[n] > 0$ is the instantaneous signal power at index $n$, then the quantity described above is:
$$ y[n] = \frac{\frac{1}{K} \displaystyle\sum_{i = 0}^{K-1} x[n-i]}{\frac{1}{L}\displaystyle\sum_{j = 0}^{L-1} x[n-j]} = \frac{L}{K} \left(1 + \frac{\displaystyle\sum_{j=K}^{L-1} x[n-j]}{\displaystyle\sum_{i=0}^{K-1} x[n-i]}\right)^{-1},$$ for fixed integers $K$, $L$ with $L \gg K > 1.$
I am interested in the behaviour of this nonlinear system, since it seems quite useful in practice. Can anyone suggest a way to analyse it theoretically?

This depends a bit on what you want to get our of the analysis
Step 1: Instantaneous  power
You square the signal. That changes the spectrum considerably. In particularly you need to watch out for aliasing. Let's say your signal is a 15kHz sine wave sampled at 48 kHz. If you square this you get one component at DC (0 Hz) and another at 30 kHz. However, the 30 kHz component will alias down to 18 kHz, so your squared spectrum is outright wrong in this case. To be safe, you should either up-sample the signal by a factor of 2 before squaring OR lowpass filter it at a quarter of the sample rate. Many "natural" audio signals don't have a lot of energy at high frequencies, so could also choose to just live with the aliasing error, but that's an application specific decision.
The physical interpretation of instantaneous power is complicated. Let's say your input signal is the voltage at the terminals of a loudspeaker: Squaring it does NOT represent the instantaneous power. The actual instantaneous power is the product of the voltage and the current waveform. Since a loudspeaker is a reactive load, this will be quite different from the squared voltage. For example, it's perfectly normal for the instantaneous power to be negative which the square of a signal can never be. Negative instantaneous power simple means that the load has energy storage and that at times it can feed stored energy back into the source.
Going back to the sine wave example. A squared sine creates a component at DC and one at double the frequency. One possible interpretation: the DC component reflects the actual power consumed, the double frequency component represents power that gets exchanged between the source and the load but not consumed: it just gets shoved back and forth.
Provided that the signal is the "actual" instantaneous power, than the mean of it will represent the average power and the integral will represent the real energy consumed (exactly).
Step 2: Moving average
Moving average is just a low pass filter where the length of the averaging window determines the "cutoff" frequency. Unfortunately it's a rather poor one. It attenuates the "higher" frequency in the spectrum but very unevenly. Personally I would assume that in most cases a simple IIR filter would work much better.
A low pass filter will eliminate the "back and forth components" of the power signal but keep the "actual consumption". The choice of cutoff frequency (or time constant) should be tied to the physical process you want to observe. I.e. could be to time constant of the power supply storage or the thermal time constant of the system.
Step 3: Dividing averages for different time constants
This creates the ratio of "short term power " vs "long term power" where you can define "short" and "long" by the cutoff frequencies of your low pass filter. Whether that's useful or not depends on your application. It's a relative metric: it shows how "spikey" your energy consumption is.
Since real instantaneous power can be negative, it's conceivable that the lowpass filtered version can be negative or zero as well, which can make the division problematic.


IIR implementation - RuntimeWarning: overflow encountered in double_scalars I have implemented a butter filter in python using this function:
def butter_lowpass(cutoff, fs, order=5):
    nyq = 0.5 * fs
    normal_cutoff = cutoff / nyq
    b, a = butter(order, normal_cutoff, btype='low', analog=False)
    return b, a

b and a:
array([1.53324552e-09, 6.13298208e-09, 9.19947312e-09, 6.13298208e-09, 1.53324552e-09])
array([ 1.        , -3.9671626 ,  5.90202586, -3.90255878,  0.96769554])

Now I have to do an implementation. When testing different algorithms, I have the following problem:
RuntimeWarning: overflow encountered in double_scalars

How to avoid this problem? Perform scaling? The problem, in my opinion, is the b coefficients.
When I call the function lfilter everything is OK
The code and data: https://github.com/Rariusz/IIR-implementation.git
The actual IIR implementation:
N = order1
L = len(data)
x = np.zeros(L,dtype=np.float64)
u = np.zeros(L,dtype=np.float64)
y = np.zeros(L,dtype=np.float64)
a = a1
b = b1

for n in range(L):

  for i in range(N,1,-1):
    x[n-i-1] = x[n-i]
  x[n] = data[n]

  y[n] = 0
  for i in range(N):
    y[n] = y[n] + b[i]*x[n-i]

  for i in range(1,N):
    y[n] = y[n] + a[i]*y[n-i]

Solution:
y = np.zeros(len(data))
x = data
for n in range(4,len(data)):
  y[n] = -(a[1]*y[n-1] + a[2]*y[n-2] + a[3]*y[n-3] + a[4]*y[n-4]) + \
         b[0]*x[n] + b[1]*x[n-1] + b[2]*x[n-2] + b[3]*x[n-3] + b[4]*x[n-4]

First, I assume that it is a that is equal to array([ 1.        , -3.9671626 ,  5.90202586, -3.90255878,  0.96769554]).
Second, that array when you go through all the z-domain stuff, you find that you are representing the difference equation
$$\sum_{k=0}^N a_k y_{n-k} = \sum_{k=0}^N b_k x_{n-k}$$
If you solve that for $y_n$ you get
$$y_n = -\sum_{k=1}^N a_k y_{n-k} + \sum_{k=0}^N b_k x_{n-k}$$
Note the minus sign in front of the summation of the prior terms of $y$ -- your code is missing that.  So you've turned your denominator polynomial from $a(z) = 1 - 3.967z^{-1} + 5.902z^{-2} \cdots$ into $a(z) = 1 + 3.967z^{-1} - 5.902z^{-2} \cdots$.  That polynomial is almost certainly unstable, and you're overflowing even what double-precision floating point can represent.
Note 1:
In general, it's better to break filters down into the smallest sections you can and cascade them.  This is because (a) it is the filter poles and zeros that are important to filter behavior, but polynomial roots get absurdly sensitive to polynomial coefficient variation as the order goes up, and (b) the filter's susceptibility to rounding in the data -- even with perfect coefficients -- goes up pretty much the same as a polynomial's sensitivity to coefficient variation.
So -- small polynomials are good, which means low-order filters are good.
Your filter design function should be able to cough up a form where it gives you poles & zeros instead of polynomials; you can then break those into 2nd-order functions.
Note 2:
You're using Scipy, which -- if I'm not mistaken -- has a built-in filtering function, and it should be compatible with the poles & zeros form of the filter design output.  If you're rolling your own filter code to learn it -- good for you, keep it up, we need people to know what's going on under the hood.  If you just need to filter some @#$% data and get on with life, use the built-in Scipy function.


Eigenfunction of LTI causal system Z-transform I am studying for my DSP final and I came across this question from the Oppenhiem and Schafer book 3rd edition. The question says

  3.18 A casual LTI system has the system function

  $$H(z)=\dfrac{1+2z^{-1}+z^{-2}}{\left(z+\frac{1}{2}z^{-1}
> \right)\left(1-z^{-1} \right)} $$

  (a) Find the impulse response of the system, $h[n]$

  (b) Find the output of this system, $y[n]$, for the input

  $$x[n]=e^{j(\pi/2)n} $$

For the first part, I got
$$h[n] =  -2 \delta[n]+\frac{1}{3}  \left(\frac{-1}{2}\right)^2 u[n]+\frac{8}{3} u[n]$$

which matches the answer at the back of the book, but the second part is difficult. We are asked to find the function $y[n]$ if $x[n] = 2^n$

I know that $2^n$ is an eigenfunction of $y[n]$, but I am not sure how to substitute it in $H(e^{j\omega})$. The answer at the back of the book is $y[n] = \frac{18}{5} 2^n$

Can anyone help me? Thanks in advance

Hint: $y[n] = H(e^{j (\pi / 2)}) x[n]$.

Since $x[n]$ is a complex exponential, it is an eigenfunction of a LTI system with eigenvalue given by the frequency response at the frequency of $x[n]$ (which in this case is $H(e^{j (\pi/2)})$. In less fancy words, if you put a complex exponential into a LTI system, you get the same complex exponential out, except scaled by the transfer function evaluated at the frequency of the complex exponential.


Pulse shaping filter for FSK After a BPSK modulation we have the upsampling and the pulse shaping filter. But for a FSK, do we have upsampling and pulse shaping filter, or just the pulse shaping filter?

Most likely yes. FSK is already generated using many samples per symbol. Look at Matlab fskmod function.


How to Classify Image Filter as Low Pass Filter (LPF) or High Pass Filter (HPF)? I am currently reading on the topics of image filtering and convolution kernels. While the process of convolution is clear to me, I do not understand, how one arrives  at the statement, that a certain filter is a low or a high pass?

How can one derive that? For example, if one averages the value of a pixel respecting its surrounding pixels, does that for all pixels, then each pixel would have a value between the minimum and the maximum value of the previous image. It treads low values and high values equally, however literate often states, that it is a low pass filter.

For that specific case I found this question here. Reading it, however didn't help me, as I seem to lack something to see how these arguments are correct and intuitive. Thus I hope, that a more mathematical approach (can be a general one), would benefit me.

Welcome to DSP Community.

General

We assume 2 types of filters: LPF or HPF.

Classifying Filter Type

Usually, if it is a well planned LPF and well Planned HPF a simple test will do.
Calculate the sum of all coefficients.
The sum of the coefficients is the first element of the DFT of the signal.
It means it is the DC gain and well behaved LPF has gain of 1 and well behaved HPF have DC gain of 0.

So sum the coefficients and:

If the sum is $ 1 $ then this is an LPF.
If the sum is $ 0 $ then this is an HPF.

Pay attention that those method can't differentiate between HPF and Band Pass or between LPF and Band Stop.

The more general method is to look on its DFT Spectrum to see its behavior.

Converting Filter Type

There is very intuitive formula - ALL_PASS_FILTER = LPF_FILTER + HPF_FILTER.
It's not accurate but this is the idea.

Be more accurate in the Frequency Domain All Pass Filter is a filter which has value of 1 for any frequency.
Assuming we have LPF Filter then HPF_FILTER = ALL_PASS_FILTER - LPF_FILTER.
In the spatial domain it means $ h \left[ n \right] = \delta \left[ n \right] - l \left[ n \right] $.

Where $ \delta \left[ n \right] $ is the filter with 1 in its origin and $ 0 $ anywhere else.

In your case $ \delta \left[ n \right] = [0, 0, 0, 0, 1, 0, 0, 0, 0] $.
So the HPF Filter becomes $ h \left[ n \right] = \delta \left[ n \right] - l \left[ n \right] = [-1/9, -1/9, -1/9, -1/9, 8/9, -1/9, -1/9, -1/9, -1/9] $.


How to extract a signal received from a given angle by two microphones Please bear with me, as I'm only a programmer, and have no knowledge of signal processing as a discipline.

I have an application which generates a short pulse of 12 kHz frequency sound from a single central source, then captures any reflected sound with a pair of stereo microphones a few cm apart.

I want to be able to, effectively, scan across a range of angles, and extract the reflected sound at each angle, in order to build a picture of what is reflecting the sound.  Essentially, sonar.

I've tried extracting the signal at different angles by sampling the incoming waveforms at different points corresponding to the time difference each angle would represent; but I'm not sure how to get something sensible from that.
(I've tried calculating what I'm calling the correlation by dividing the absolute difference between the left and right signals and dividing by the maximum possible to normalise it; but all I seem to get, if I plot it, is noise).

I've spent a couple of days trying to find it online; but I'm hampered by my utter lack of knowledge of what even the most basic terms mean, and it's also a long time since I studied maths.

Cross correlation is most suitable for low SNR long term observation of continuous signals.

I recommend you try to identify the leading edges of your short duration received  Signals. Leading edges are usually direct paths. The classic text book solution is to use a matched filter, but often a high Signal level is enough. With the speed of sound and time difference of leading edges, simple geometry should give you an estimate of the direction of arrival.

In actual SONAR a single multi mode  ring hydrophone produces an arrival direction. I don't know if there is a comparable sensor in air acoustics.

You should also try to see if there are built in offsets in your collections. Not all audio cards have actual multichannel sampling.

Good luck


How to improve accuracy while converting floating point coefficients to fixed point in the case of an all pole IIR filter I am having an all pole IIR filter of order 5 with floating point filter coefficients
[1,
-4.708642687,
9.963818327 ,
-11.99038368,
8.647611478,
-3.553230696,
0.659454607].

When an impulse input was applied to the IIR filter, the first 8 output samples obtained were
1,
4.7086,
12.20,
22.5509,
32.38,
36.9400,
32.26087.

To implement the same system in hardware, I have scaled the filter coefficients by a factor of $2^{10}$ and down-scaled back at the output. Similarly all the operations were done in fixed-point representation. However, the first 8 output samples obtained at the output were
1,
5,
14,
29,
49,
70,
86,
91,
82 etc.

How to reduce this large difference in the output when trying to convert from floating point to fixed point

Where are the poles of your filter? If they are close to the unit circle, then quantizing them can cause some of them to cross beyond the unit circle, causing the quantized filter to be unstable. You should check the poles of your filter before quantizing and after quantizing. See how close they match.
Secondly, it's risky, albeit not impossible, to implement an order-5 IIR filter directly in fixed-point arithmetic. Usually, it's a good idead to split it in cascaded order-2 filters, each of these order-2 IIR filters is usually called a "biquad". In your case you would have 2 cascaded biquads with one order-1 IIR filter.
There are many references that you can check out such as this one.
https://www.dsprelated.com/showarticle/1137.php


Using NMF for guitar pitch detection using fixed pitch template I've seen this guide showing a basic NMF of a spectrogram.

My aim is detect guitar pitch from the decomposition. I want to do this by explicitly create a fixed template of pitches and construct the W matrix from these.

I've seen some papers giving an overview of this, but I am not sure how it's practically done.

Do I record every note on the guitar, then do the FFT for each, take its spectrum, and paste it inside a matrix?

I guess my question is basically: how do I create the fixed W matrix, given a set of recordings with each pitch.

This paper says:

  The whole note sample k is first processed in a short-time sound
  representation supposed to be non-negative and approximatively
  additive (e.g. a shorttime magnitude spectrum). The representations
  are stacked in a matrix $V^{(k)}$ where each column $v^{(k)}_{j}$ is the sound
  representation of the j-th time-frame.

So, does that mean we compute the STFT of the recorded sample? Then we use the resulting spectrogram as $V^{(k)}$?

Also, does the duration of the sample matter? If there is absolute silence in the beginning of the sample, will that affect $v^{k}_{j}$?

The goal of non-negative matrix factorization (NMF) is to factorize (as the name says) a matrix, denoted by $\mathbf{X}\in\mathbb{R}^{(N\times K)}$ in the reference you mention, into the product of two matrices:
$$
\mathbf{W}\mathbf{H}\approx\mathbf{X},
$$
with $\mathbf{W}\in\mathbb{R}^{(N\times M)}$ the template, dictionary, or basis matrix and $\mathbf{H}\in\mathbb{R}^{(M\times K)}$ the activation matrix. The approximation sign ($\approx$) in the equation indicates that in many cases it is not possible to actually find a true factorization but rather an approximation.

Inspecting the above equation reveals that each column of $\mathbf{X}$ is approximated by a linear combination of columns of $\mathbf{W}$, with the coefficients of linear combination in $\mathbf{H}$. (Experimenting with some simple examples, e.g., with $N=3,\,M=3, K=5$ and $\mathbf{X}_{ij}\in\left\{0,\,1\right\}$, will surely help gaining insight into this fundamental behaviour).

In many cases it is desired, given a matrix $\mathbf{X}$, to determine both the dictionary matrix $\mathbf{W}$ and the activation matrix $\mathbf{H}$. Optionally, of course, only one can be estimated if the other is given.

Your goal, if I understand correctly, is to keep $\mathbf{W}$ fixed, estimate $\mathbf{H}$ and use the resulting activations to determine when which note has been played. That should work as a first attempt, I think.

After this coarse overview now for your questions:

You are right: to create the $\mathbf{W}$ matrix put an estimate of each note's power spectral density estimate (e.g., obtained with Welch's method) into a column of $\mathbf{W}$ (I'm not sure what you mean by "do the FFT for each, take its spectrum".).
Yes, the matrix $\mathbf{X}$ in the notation used in this answer ($V^{(k)}$ in the paper that you mention) contains the STDFT magnitude spectrogram of the full input recording.
Silence in the beginning of the recording will obviously have an influence on the first columns of $\mathbf{X}$ (or, respectively, $v_j^{(k)}$ for $j<\text{some number}$). However, what do you guess will happen with $\mathbf{H}$ if the first columns of $\mathbf{X}$ contain zeros or only very small values? Correct: the first columns of $\mathbf{H}$ will also receive zeros or very small values.

I hope this answer helps in clarifying the basic functionality of the NMF algorithm. A rather intuitive description of the NMF idea is also given in this famous paper.


Bilinear Transformation Comparison If I have transfer function coefficients,
I can analyze the transfer function in the s-plane and/or the z-plane.

If I wanted to demonstrate that the z-plane and s-plane responses are equivalent:

Would it be possible to use the frequency warping relation
(provided in the otherwise unused variable w.z in in the code)
to depict the z-plane plots in the s-plane or vice versa
in a manner in which the two would completely match?

Matlab code provided:

[Actual filter coefficients not provided.]

%% Input

A.s.cf = [ 0 0 0 0 1 ];
B.s.cf = [ 1 2 3 4 5 ];

T.z   = 2.5e-05;
f_z   = 1/T.z;
f_nyq = f_z / 2;

%% Processing: Matlab SP toolbox
%{
H.s.tf = tf ( A.s.cf, B.s.cf       );
H.z.tf = c2d( H.s.tf, T.z,   'zoh' );
%}

%% Processing: From scratch:
f.s = logspace(-2, +6, 5e3).';
w.s = 2*pi*f.s;
s.s = 1j * w.s;

z.z = exp( 1j * w.s * T.z);

w.z = 2/T.z * atan( w.s * T.z/2 );

switch 0
  case 0; s2z = 1/T.z * log(z.z);
  case 1; s2z = 2/T.z * (z.z-1)./(z.z+1);
end

A.s.jw = s.s.^4 * A.s.cf(1) ...
       + s.s.^3 * A.s.cf(2) ...
       + s.s.^2 * A.s.cf(3) ...
       + s.s.^1 * A.s.cf(4) ...
       + s.s.^0 * A.s.cf(5) ; %

B.s.jw = s.s.^4 * B.s.cf(1) ...
       + s.s.^3 * B.s.cf(2) ...
       + s.s.^2 * B.s.cf(3) ...
       + s.s.^1 * B.s.cf(4) ...
       + s.s.^0 * B.s.cf(5) ; %

A.z.jw = s2z.^4 * A.s.cf(1) ...
       + s2z.^3 * A.s.cf(2) ...
       + s2z.^2 * A.s.cf(3) ...
       + s2z.^1 * A.s.cf(4) ...
       + s2z.^0 * A.s.cf(5) ; %

B.z.jw = s2z.^4 * B.s.cf(1) ...
       + s2z.^3 * B.s.cf(2) ...
       + s2z.^2 * B.s.cf(3) ...
       + s2z.^1 * B.s.cf(4) ...
       + s2z.^0 * B.s.cf(5) ; %

H.s.jw = A.s.jw ./ B.s.jw;
H.z.jw = A.z.jw ./ B.z.jw;

If I plot the two,
the response will be different in either plane,
(as expected).

%% [Plot]: Setup
p = 0;

%% [Plot]: S plane
p = p+1;
figure(p)

% S plane: Magnitude
subplot(2,1,1)
semilogx(w.s, 20 * log10( abs(H.s.jw) ), '.')

title('Bode [From Scratch]')
ylabel('Mag [dB]')
xlabel('w.s [rad/s]')

axis( [ 1e+3 1e+6 -150 +0] )
grid minor

% S Plane: Phase
subplot(2,1,2)
semilogx(w.s, angle( H.s.jw ) * 180/pi, '-')
hold on
semilogx( [ 2*pi*1e-2 2*pi*1e+6 ], [ -180 -180 ], 'k')
semilogx( [ 2*pi*1e-2 2*pi*1e+6 ], [ +180 +180 ], 'k')
hold off

ylabel('Ang [deg]')
xlabel('w.s [rad/s]')

axis( [ 1e+3 1e+6 -200 +200] )
grid minor

%% [Plot]: Z Plane
p = p+1;
figure(p)

% Z Plane: Magnitude
subplot(2,1,1)
semilogx(w.z, 20 * log10( abs(H.z.jw) ), '.')
hold on
semilogx([2*pi*f_nyq 2*pi*f_nyq], [-50 +50], 'k-')
hold off

title('Bode [From Scratch]')
ylabel('Mag [dB]')
xlabel('w.z [rad/s]')

axis( [1e+3 1e+6 -50 0] )
grid minor

% Z Plane: Phase
subplot(2,1,2)
semilogx(w.z, angle( H.z.jw ) * 180/pi, '-')
hold on
semilogx( [ 2*pi*1e-2  2*pi*1e+6 ], [ -180 -180 ], 'k')
semilogx( [ 2*pi*1e-2  2*pi*1e+6 ], [ +180 +180 ], 'k')
semilogx( [ 2*pi*f_nyq 2*pi*f_nyq], [ -180 +180 ], 'k-')
hold off

ylabel('Ang [deg]')
xlabel('w.z [rad/s]')

axis( [1e+3 1e+6 -200 +200] )
grid minor

%% [Plot]: S Plane: Bode

%{
p = p+1;
figure(p)

bode(H.s.tf)
grid minor
%}

%% [Plot]: Z Plane: Bode

%{
p = p+1;
figure(p)

bode(H.z.tf)
grid minor
%}

i really don't wanna slog through your MATLAB code.  is this what you're trying to prove? :

let $H_\text{a}(s)$ be some analog (or "continuous-time") LTI transfer function on the s-plane and $H_\text{d}(z)$ be some digital (or "discrete-time") LTI transfer function in the z-plane.

$$ H_\text{d}(z) = H_\text{a}(s) \Bigg|_{s=\frac{2}{T}\frac{z-1}{z+1}} $$

do you want to prove that

$$ H_\text{d}(e^{j\omega}) = H_\text{a}(j\Omega) \Bigg|_{\Omega=\frac{2}{T}\tan\left(\frac{\omega}{2}\right)} $$

is that what you want proven?


Block diagram reduction I have the following block diagram

I know that it's solution is:

$$\dfrac{10 K_p}{Js^2+(B+2K_pK_v)s+10K_p} $$

But when I try to derive the solution with reduction rules I get different result.

Assumptions:

$5$ and $K_vs$ are in parallel: $(5+K_vs)$

$K_p \cdot \dfrac{2}{(Js+B)s}$

result: $$\dfrac{K_p \cdot \dfrac{2}{(Js+B)s}}{1 + K_p \cdot \dfrac{2}{(Js+B)s} + (5+K_vs)}$$

What am I doing wrong?

I believe that it will be worth for many people to do it from scratch. Let's call the signal after input and feedback summation $E(s)$:

Then question arises what $E(s)$ is equal to? Well that's simple:

$$E(s)=5R(s)-(5+K_vs)C(s) $$

So we have our backward path, now the forward one. What $C(s)$ is equal to?

$$C(s)=E(s)K_p \dfrac{2}{(Js+B)s} $$

Let's substitute previously calculated $E(s)$ to the above equation:

$$C(s)=\dfrac{2}{(Js+B)s}K_p[5R(s)-(5+K_vs)C(s)] $$

Because we have our output $C(s)$ on both sides of the equation, we must eliminate that. Expand the right side of the equation:

$$C(s)=\dfrac{2}{(Js+B)s}5K_p\cdot R(s)-\dfrac{2}{(Js+B)s}K_p(5+K_vs)C(s) $$

Now move the part with $C(s)$ to the left side and move it before parenthesis, this yields:

$$C(s) \left[1+ \dfrac{2}{(Js+B)s}K_p(5+K_vs)\right]=\dfrac{2}{(Js+B)s}5K_p\cdot R(s) $$

Now you can divide right side of the equation by the factor from the left side:

$$C(s) =  \dfrac{\dfrac{2}{(Js+B)s}5K_p}{1+ \dfrac{2}{(Js+B)s}K_p(5+K_vs)}\cdot R(s) $$

Great! Now we have direct relation between the input and the output. Let's make this fraction more pleasant:

$$C(s) =  \dfrac{\dfrac{10K_p}{(Js+B)s}}{1+ \dfrac{10K_p+2K_pK_vs}{(Js+B)s}}\cdot R(s) $$

Now expand the one in denominator to: $\dfrac{(Js+B)s}{(Js+B)s} $, and add it:

$$C(s) =  \dfrac{\dfrac{10K_p}{(Js+B)s}}{\dfrac{(Js+B)s + 10K_p+2K_pK_vs}{(Js+B)s}}\cdot R(s) $$

We can now remove common $(Js+B)s$:

$$C(s) =  \dfrac{10K_p}{(Js+B)s + 10K_p+2K_pK_vs}\cdot R(s) $$

Finally let's do one last multiplication:

$$C(s) =  \dfrac{10K_p}{Js^2+Bs + 10K_p+2K_pK_vs}\cdot R(s) $$

And voila!

$$C(s) =  \dfrac{10K_p}{Js^2+ s(B +2K_pK_v)+ 10K_p}\cdot R(s) $$


Is sinusoids the nature's default signals? I know the basic working principle of sinusoidal oscillators. But I have a doubt on this matter when I came across Gram-Schmidt orthogonalization procedure. I have read that the source for sinusoid signal output from an oscillator is the noise. Due to controlled positive feedback a single frequency component is selected and amplified, right? If the noise is a white noise, we say it has all frequency components. We think in terms of all possible sinusoids with different frequencies. When coming to frequency response too, we think in terms of sinusoids. Why can't we think in terms of other orthogonal basis set? Although it makes life difficult, I would like to know the justification for thinking in terms of sinusoids. Is sinusoids the nature's default signals?

Time invariance plays a huge role in nature. Most systems (including your ear/brain) don't have an absolute time reference but treat all points in time equally. That results in a preference for the description of these systems with essentially time invariant basis functions, which is what (complex) sinusoids are.

For linear time invariant systems, the complex exponentials are even eigenfunctions of the system and allow you to describe the system as a "diagonal" operator, i.e. in a very simple way. That even includes a system that just performs a time delay.

Nonlinear time invariant systems can also be described in terms of sinusoids, but it gets way more complicated there. A nonlinear time invariant system excited by a periodic input will produce a periodic output in general, and you can analyze both in terms of sinusoids.

But you can also perform a more mathematical rigorous analysis by looking at the nonlinear system response to small disturbances of any input signal. And taking sinusoids to describe these disturbances gives you a local description in a time invariant way.

So all in all, if you have a system that is time invariant, no matter of what nature, time shift invariant basis functions are helpful, and complex sinusoids are precisely that.


Shannon-Nyquist theorem reconstruct 1Hz sine wave from 2 samples Lets say I want to set the minimal sampling rate to reconstruct a 1Hz sine wave, according to the Nyquist-Shannon theorem that states that the maximum recoverable frequency is Fs/2 i.e. we must sample the signal 2 times the maximum frequency.

It seems obvious that the limit is a sampling frequency of 2Hz. Lets even say I interpret it as just more than 2Fs, lets say 3Fs. This would imply that 3 samples is enough.

So, in Matlab i generate:

f=1;%my 1Hz freq
Fs=2*f+1;%=3 =>a bit more than the Nyquist freq
t = linspace(0,1,Fs);
y=sin(2*pi*f*t);
subplot(1,2,1);
plot(t,y);
title('sin(2*pi*1*t)');
subplot(1,2,2);
plot(abs(fft(y)));
title('abs(fft(y))');

The ''sine'' does not even go back up to 0, not even mentioning that even if it did it would be a saw-tooth wave rather than a sine but i guess that is not the problem.

What am I missing, why do I need at least 4 samples rather than 2?

I think this is important to understand the theorem 'in practice'.

Although even with 4 samples the spike in the FT is wrong, it is at 3 rather than 2 (2 is the 1Hz since the first would be the DC freq)

In fact I need 5 samples to get finally this spike at 2 in the FFT! Why ?

It's your observation interval which creates the main problem.

Your reasoning based on the Nyquist sampling theorem is ok; of course with a pure sine wave at the exact Nyquist frequency you will have troubles and therefore it's wise to relax the sampling frequency (slightly) above that of Nyquist rate, such 2.2 Hz instead of a strict 2 Hz... So this is one problem you will practically have.

But your main problem that appears on the FFT plot is about the spectral resolution due to short observation interval. Mainly with two samples (or one second observation of a 1 Hz sine wave) your FFT resolution will be limited to even less than one Hz. Please search the site for PSD, spectral resolution, FFT bin frequency to get a better understanding of spectral analysis of practical windowed data.

In order to see sharp frequency peaks (ideally impulses) in the FFT output of your sine wave, you should increase the spectral resolution, which requires you to increase the observation interval.

I have modified and also expanded your code to see the result of ideal sinc based interpolator on the (near) critically sampled data. Note that I've mereley included a digital simulation of ideal sinc based interpolator (not the simulation of analog interpolator) to see that it would indeed reconstruct the pure sinudosidal from its given samples taken at close to Nyquist rate. Note that for the ideal sinc interpolator to work, the original signal should strictly be bandlimited, or at least sufficiently so, which will have a lot of consequences on the success and efficiency of interpolation.

f  = 1;                % 1 Hz. sine wave...
Fs = 4.2*f;            % sampling frequency Fs = 2.2*f  ; a bit more than the Nyquist rate.
Td = 25;               % duration of observation ultimately determines the spectral resolution.
t  = 0:1/Fs:Td;        % observe 25 seconds of this sine wave at Ts = 1/Fs
Td = t(end);           % get the resulting final duration
L  = length(t);        % number of samples in the sequence
M = 2^nextpow2(10*L);  % DFT / FFT length (for smoother spectral display, not better resolution! )

x = sin(2*pi*f*t);     % sinusoidal signal in [0,Td]
%x = x.*hamming(L)';   % hamming window applied for improved spectral display

% Part-II : Approximate a sinc() interpolator :
% ---------------------------------------------
K = 25;                  % expansion factor
xe = zeros(1,K*L);       % expanded signal
xe(1:K:end) = x;

D = 1024*8;
b = K*fir1(D,1/K);     % ideal lowpass filter for interpolation

y = conv(xe,b);
yi = y(D/2+1:D/2+K*L);

subplot(3,1,1);
plot(t,x);
title(['1 Hz sine wave sampled at Fs = ',num2str(Fs),' Hz, Duration : ', num2str(Td), ' s'])
%xlabel(' time [s]');

subplot(3,1,2);
plot(linspace(-Fs/2,Fs/2-Fs/M,M),fftshift(abs(fft(x,M))));
title(['magnitude of ', num2str(M), '-point DFT / FFT of y[n]']);
%xlabel('Frequency [Hz]');

subplot(3,1,3)
plot(linspace(0,Td,length(yi)),yi);
xlabel('approx simulation of ideal sinc interpolation');

Below is a plot for the result of interpolation from a set of near critical sampling.

And below is the same simulation with a more relaxed sampling, as you can see the interpolator performs much better for this set of improved samples (consequence of better bandlimitedness)


How do I eliminate line noise from USB microphone I'm trying to record voice from a microphone into my laptop.  Unfortunately, the laptop does not have a microphone jack, so I have to use the USB input.  This doesn't present a problem when I record voice from my gaming headset, but when I record from the condenser microphone over USB I get a very pronounced line hum.  This occurs even when the laptop is unplugged from power and on its battery.

The hum consists mainly of 55 Hz noise, but with some 110 and 512 Hz components. I can eliminate this using a notch filter in "post production", as it were, but I'd rather not have to do this.  This noise is not the fault of the microphone, by the way (at least not directly), since when I run the microphone signal into my digital voice recorder (using a 3.5mm plug) there is no background hum at all.

Is there some way to process or filter the output from the microphone so that this noise never makes it to the laptop?

You already know the noise frequencies, if you know the phases, you can cancel the noise by adding mods with opposite phase, like adaptive headphones do.


Phase wrap when computing impulse response from frequency response Short:

I have the frequency response of a system as a set of complex numbers.
I want the impulse response for that system.

Problem: The phase of my frequency response wraps around 2pi multiple times. When I do a IDFT of the frequency response to get an impulse response, how can I be sure the group delays match ones of the original system?

Long:

My naive approach would be to use the sampled frequency response and feed that into an IDFT/IFFT to obtain an impulse response. But the phase information in the complex numbers is quantized to 0..2pi. For example, let system A have an arbitrary frequency response. System B has the same response, but with an additional phase shift of 360° for all frequencies above a certain threshold f0. In the space of their frequency responses in the form of complex numbers, both systems would look identical. If I did the IDFT of system Bs frequency response, I would get an impulse response that doesn't reflect that 360° phase shift because it wasn't reflected in the complex data. Is this correct?

Staying with my actual problem:

For the system in question, I can clearly see from my phase response, that the phase wraps around 2 Pi multiple times. I can reconstruct the actual phase response by unwrapping it with Matlabs unwrap() command as the phase response is relatively smooth. But as soon as I feed this into an IDFT/IFFT I have to quantize the phase into the angle of a complex number again and the information gained from unwrap() is lost again.

For my application, the actual group delays are very important. So I'm looking for a way to reconstruct the correct impulse response, taking into account that the phase information for a certain frequency may actually have an additional N full rotations that are not reflected in the complex numbers.

Or am I simply imagining the problem and systems like system B are actually impossible?

No problem at all. Inverse Fourier Transform is totally blind to phase wrapping. You can add any multiple integers of two pi to any phase and you will still get the exact same impulse response


How to implement Matlab welch in python with the same result? I want to implement Matlab pwelch function in python.
The Matlab code is shown blew:

[PSD,f] = pwelch(data,length(data),[],[],2000);

which data is a row vector of double type with size (1 * 1547600)

My python code is:

f, Pxx_den = signal.welch(data, fs=2000, window=signal.hanning(len(data)), return_onesided=True, detrend=False)

The output of python code doesn't equal to the output of the Matlab code and I don't know how to set python welch arguments to reach Matlab output results.

Thanks for your attention.

As written in Matlab pwelch document, by using an integer for window, a Hamming window of length window(integer) will be used. Also the nfft in python should set to the next power of two greater than length of window.

freqs, PSD = signal.welch(data, fs=2000, window=signal.windows.hamming(len(window)), nfft=num_fft, scaling='density', return_onesided=True, detrend=False)


difference between wavelet filters Some may think this question is little primarily,I just want to know the differences between different wavelet filters like decomposition high pass,low pass,reconstruction high pass,low pass filters..Which one is best suited for an ECG type signal

Comparative Analysis of Wavelet Basis Functions for ECG Signal Compression through Compressive Sensing, 2012 provides a comparative analysis for different wavelets. On a more generic basis, Compressed Sensing for Bioelectric Signals: A Review, 2015, addresses other kinds on splines, Gabor wavelets, etc. The references therein could provide you with a good starting point.

Fixed percentage of wavelet coefficients to be zeroed for ECG ccompression, 2003, could be interesting as well, but I cannot access the paper.


Filter away sinusoidal noise properly I have a stereo music signal corrupted by strong sinusoidal noise that varies over time. Here is the spectrogram of Left channel I plotted with Matlab.
As you can see there are 3 or 4 strong harmonics with frequency that varies over time.

As a first try I computed the difference of left and right channels and it seems most of the harmonic noise disappears, so my guess is the noise in left and right are almost identical. Here is the difference spectrogram to prove this.

Which technique could I use now to filter away all the interfering harmonics while preserving the original music spectrum as well as possible?

Notch filters?
Tracking the frequencies over time and adding a sinusoid with reversed phase?

How do I estimate and use the phase information? I think phase is important to have good denoising, but would not know how to proceed.

Thanks a lot! :)

Definitely, try to reconstruct the sine and then subtract it.  A notch filter will create artifacts as the desired signal changes, with notes stopping and starting, and percussion. Even if those notes are at pitches of frequencies outside the range of the notch filter.

To extract the unwanted sine, use a bandpass filter. While this has the same problem as the notch filter, with notes stopping and starting, you do know that the sine is steady (or slowly varying) over time, while the desired signal is always changing.  Or may include some steadily droning elements - bagpipes, anyone?   Hope that's outside the passband of the filter.  You also have found how the sine and the signal related between left and right, so make use of that to help extract the sine.  The electronics engineer part of my brain wants to send the extracted sine to a PLL, to filter out any accidentally included signal, and provide smooth amplitude and phase data. Some sort of mathematical flywheel to steady the cycles.

The exact methods for extracting and purifying the sine depend on if you need to do this in real time, or near real time with a fraction of a second delay acceptable, or if you have the whole span of signal from beginning to end all at once such as with an .mp3 file.

I know this is the best way to do it, because Cassini's NAC camera had a similar problem. You have the audio version of what I had to fix years ago.


Can FFT information be used to obtain the DTFT? If $x[n] = \left\{\frac 14, \frac 14, \frac 14, \frac 14\right\}$
and the resulting DTFT is
$$
\frac 12\exp\left(\frac 32\omega\right)\left(\cos\left(\frac 32\omega\right)+\cos\left(\frac 12\omega\right)\right)\tag{1}
$$

I'm looking for a general, computational method for obtaining the general DTFT formula, that works in almost all cases, that doesn't involve doing anything by hand. Is there enough information in a FFT magnitude and phase result, assuming there were enough bins used and appropriate zero-padding was done, to work out the DTFT?

To reiterate, pretend like all you have is the graph below. I'm asking if you can somehow derive the formula equation $(1)$ without doing anything else but looking at that graph. I'm asking if that is possible using just FFT information.

Your question is a bit difficult to comprehend. It seems that you are requesting a closed form solution to the DTFT given the DFT of a sequence. This is of course possible. If there is no time-domain nor frequency-domain aliasing, then we can simply apply the formula for the IDFT on the samples (or equivalently perform the IFFT). Following this, we can apply the formula for DTFT on the time-domain samples. This may not be an answer in the simplified form you want, but it is a closed form expression for the DTFT and can be evaluated at any point. Furthermore, a very good computer algebra system (e.g., Maple, Mathematica, etc.) may be able to reduce it to the form that you have requested (meeting your requirement for no computation by hand).

The DFT is a frequency-domain sampled version of the DTFT. The FFT computes the DFT. Because of this, you could simply apply the Shannon-Nyquist Sampling Theorem and interpolate the points of the DFT using sinc functions,

$$\DeclareMathOperator{\sinc}{sinc}
\begin{align} X_{DTFT}(f) &= \sum_{k=-\infty}^{\infty} X_{DFT}[k] \sinc(f - k) \\
&= \sum_{k=0}^{N-1} X_{DFT}[k] \frac{\sin(\pi(f-k))}{N \sin\left(\tfrac{\pi}{N}(f-k)\right)}
\end{align}$$

where $\sinc(x) = \frac{\sin(\pi x)}{\pi x}$. This expression is equal to your simplified expression (i.e., it evaluates to the same value at every point). This expression assumes a unit length sampling interval (i.e., $T = 1$). This at least suggests a method for interpolating the DTFT from the DFT samples.

Looking for an exact simplified expression is difficult because different schools of thought will simplify it in different ways. Who is to say that your expression involving trigonometric functions is any better than one involving complex exponentials? For example, isn't the expression
$$ X_{DTFT}(f) = \frac{1}{4} + \frac{1}{4} \exp(-j 2 \pi f) + \frac{1}{4} \exp(-j 4 \pi f) + \frac{1}{4} \exp(-j 8 \pi f)$$
exactly what you are after? This is exactly what we would get if we used the approach in the first paragraph. Take the IFFT of the DFT, and then apply the formula for the DTFT. A computer program could do all of this and could evaluate this expression at any point $f$, requiring no computations by hand.


cancellation controller transfer function design This chapter offers a free overview on cancellation controllers, https://link.springer.com/chapter/10.1007%2F978-3-662-02319-8_6
However, i did not get equation (6.4),  how the transfer function of $G_R$ is determined?


Basic question: Why is the output of a system the convolution between the impulse response and the input? I forgot a very simple fact and I am now struggling to find reference that proves this basic property?

How would you prove that for a single in single out system, the system output is the impulse response convoluted with the input?

Because it is the response of the system when an unit impulse (delta function) is applied to the input of the system. So basically, If you  multiply this output impulse response of the system to each of the input samples and then add them all, you get the overall output of the system. And this is only true, If the system agrees with the linearity properties such as homogeneity, additivity and shift invariance. For more details and to understand the mathematics behind this, you can go through this link: http://www.dspguide.com/ch6.htm

I tried to put this as comment instead of answer, since its not detailed enough, however, due to my lack of reputation, I could not comment. Anyways, the link I have provided is very useful. Hope it helps!


Conceptual question on entropy and its relation to information Learning
Informative
Statistics:
A
Nonparametric
Approach  paper presents an approach to parameter estimation by entropy minimization. There are other related works "Minimum-entropy estimation in semi-parametric models"  download link ( http://dl.acm.org/citation.cfm?id=1195853). The rationale provided is that minimization of error entropy is equivalent to maximization of likelihood. I am new to this area and find it hard to understand the intuition behind why entropy of the error minimization will yield the parameters. What happens when entropy is minimized?

What happens when Shannon Entropy is maximized? Entropy (Shannon's) is the uncertainty  = average information or uncertainty (unsure).
And what happens when entropy is minimized and
What is the meaning of minimizing entropy of error?

The basic idea behind maximum entropy models is that you want to make the least assumption about the data. This is considered equivalent to retaining as much unpredictability as possible, as quantified by entropy. For more information I refer you to this Wikipedia article.

Information and unpredictability are closely related. If you think of the signal as a stochastic process, its information content is determined by unpredictable it is. At one extreme, if the signal tends towards being deterministic, you can tell what it will be any an arbitrary time so there is no value in sampling the signal; it has no information content. Entropy formalizes this notion.

To find the density that maximizes the entropy we have to use the calculus of variations, and that involves determining several coefficients called Lagrange multipliers. This is where the moments come in: they're the constraints that determine the multipliers. We set the $n$'th absolute moment for n=1..N to constants and we get the corresponding maximum entropy distribution. For example, the maximum entropy distribution whose first two moments are $\mu$ and $\sigma^2$ is the Gaussian distribution $N(\mu, \sigma^2)$. For details of the derivation please refer to a textbook.

I haven't read the paper but I'm guessing the rationale is that if you're trying to fit a model, you want to capture as much information of the signal as possible. Equivalently, you want to leave as little residual information as possible on the table; i.e., in the error. For details see 3.4 Minimum error entropy algorithm in Information Theoretic Learning: Renyi's Entropy and Kernel Perspectives.


Measurement of moving objects using computer vision I am working on a project to measure lengths of the black objects and distance between them using computer vision (Shown below in the fig.). These objects are moving on a conveyor shown in the blue color. The speed of the conveyor belt is known(0-1m/s) at every instance.Due to limited field of view of the camera, its not possible to have complete object length in a single frame , how can I measure the object lengths in this case?
As I have very little knowledge of computer vision, I am finding it difficult to search the algorithms relevant for solving this problem. Any help with keywords, algorithms, resources is appreciated.

My first instinct is to use line detection algorithms for detecting object edges and calculate distance between them, based on speed and time.Also, how difficult is it to stitch multiple images based on speed to get a complete object?

how difficult is it to stitch multiple images based on speed to get a complete object?

Trivial. If you know the speed and the time at which an image was taken, it's just a linear translation...

  My first instinct is to use line detection algorithms for detecting object edges and calculate distance between them, based on speed and time.

Excellent approach!

Just: reduce your image to a single row, maybe the average of the 10 center rows from your picture.

That reduces your problem from 2D to 1D.

Then, the question you have is "when does that start and stop being black?", which is trivial to answer if you just convert that color row to its brightness, and use a threshold.


Get problem to understand "Topological Structural Analysis of Digitized Binary Images by Border Following" Suzuki Abe I'm confuse with their explanation about definition 1 notation of 4- (8-) and 8- (4-) in Suzuki Abe paper "Topological Structural Analysis of Digitized Binary Images by Border Following" pg 34, and what they mean about NBD and LNBD in pg 35? Is NBD=neighborhood?

In the paper 4-(8-) refers to the 4 and 8 neighbor connectivity between the two pixels and in the algorithm NBD is the newest border found during the raster scan of the binary image and LNBD refers to the last newest border, NBD and LNBD is stored to obtain the surroundness relation between the borders.


How do I convert a sampled complex signal to a real signal? I’m writing a Soapy SDR plugin for the FL2K-SDR.  Soapy SDR expects plugins to produce and consume IQ samples. However, the FL2K-SDR driver consumes only real samples. I know you can convert a real signal to a complex signal using a Hilbert transform. How do I convert IQ samples to real samples?  Do I use an inverse Hilbert transform? Do I just drop every other sample to get only the I samples? One of my goals is to minimize the processing latency of this plugin (I.e. use lookup tables if possible).

Part of my quandary is that I don't know how the user is going to configure their flow-graph. They may:

Simply feed me a complex cosine and expect some sort of real cosine output.
Sample the PC microphone, convert that to a complex signal, and send that to me (no modulation involved).
Attach an analog low-pass filter (~ 70 MHz), do their modulation in their flow-graph, and use the FL2K USB VGA dongle for an HF Ham radio transmitter
Generate an IF in their flow graph, and then attach an analog band-pass filter to transmit an image frequency higher than 70 MHz.

Unfortunately, I don't know which of these they might do.

So the point is that "classically", communications theory tends to be done in complex baseband, i.e. signals centered around 0 Hz, but not necessarily symmetric in spectrum.

When you want to represent such signals, you need complex values. When you want to transmit such signals over the air, you need to:

Convert them from digital to analog,
mix the I(nphase) and Q(uadrature) components with the carrier (say, a cosine) and a 90° shifted version of the carrier, respectively,
and add up the result.

Since I and Q themselves are just real values, the result of 3. is just a real signal.

Now, if you do the above in the order that I've written it in, then you need two DACs running in step at the baseband sampling rate, and you can produce any RF frequency that the carrier synthesizer can generate.

But, you don't have to do it in the order 1., 2. and last 3.: You can also first digitally upsample to a higher sampling rate, then multiply the complex baseband signal with a complex sinusoid of frequency $f_\text{IF}$ to move it up to an intermediate frequency, and just discard the imaginary part of the signal. Then, you DAC just that real signal.

You could then use an LO of frequency $f_\text{carrier}-f_\text{IF}$ to move it to the final carrier frequency. But:

The FL2K is just a DAC. So it doesn't have a mixer. So you can't do that.

However, every DAC produces images at multiples of the sampling rate (that's why you always need an analog reconstruction filter, which classically is a low-pass cutting off at half the sampling rate). So, if, say your IF is $f_\text{IF}=10\,\text{MHz}$ and your sampling rate is 50 MHz, then you can produce a signal at 210 MHz simply by not low-passing the output of the DAC to below 25 MHz, but bandpassing for say 200-220 MHz.


Audio codec analysis in MATLAB: crosstalk I am currently making an analysis of a codec system that can encode up to 2 wave files into 1 file and decode it. Here is a simple illustration:

L+R ---> ENC --->L'+R'
with L'=L+∂1 & R'=R+∂2

The codec is not lossless and we are analyzing if there is any leakage from one source to another, some kind of crosstalk, if you will.

So what I have to do is find out if ∂1 has traces of R & ∂2 has traces of L.

I feel a general approach would be to subtract L/Ldec to get ∂ and then compare that ∂ to R.

I have done some reading on correlation but it's all a little vague at the moment.

So, onto my question(s):

Is this possible with limited knowledge of DSP (reading up but it's a
lot & not that easy)?
If not, on what should I do some reading (From what I've gathered so far: correlation & Pearson's coefficient)?
Is there some sort of standardized test for this - I presume there is
but I can't find any!

Even though code snippets would be handy since I only just started to use MATLAB & still am somewhat inexperienced, the main point of this topic is understanding so I can implement it. If someone else implements it I won't learn!

Most lossy stereo encoders (MP3, AAC, etc) operate in a joint-stereo mode, i.e. there will be a lot of cross talk between the two channels specifically at higher frequencies and in general they will allocate more bits to the sum then do to difference in the channels which also introduces cross talk.
The amount of cross talk will depend heavily on the signal and the bit rate. There isn't an east metric for it.
Subtracting the decoded version from the original is also not a great idea since you can't easily distinguish the encoding/decoding artifacts from the cross talk.


DFT of Transient vs. Analytical Solution I'm trying to compute the Fourier Transform of a lightning impulse waveform.  The python code I'm using is below.  The problem that I am running into is that the magnitude of the DFT and the magnitude of the analytically derived Fourier Transform are off by a factor equal to 1 / time period, where the time period is the length of time that I am using to model the signal for the DFT.

from math import *
import numpy as np
import matplotlib.pyplot as plt

def Impulse(t, V_0=1.04, alpha=0.0146e6, beta=2.467e6, t_chop=float("inf")):
    if t<=t_chop:
        return V_0 * (exp(-alpha*t)-exp(-beta*t))
    else:
        return 0

npImpulse = np.vectorize(Impulse, excluded=['V_0', 'alpha', 'beta', 't_chop'])

# Sampling rate and time vector
total_time = 500e-6 # seconds
n_samples = 2**24
sampling_rate = n_samples/total_time # Hz (samples/s)

t_values = np.linspace(0, total_time, n_samples)

# Time constants
tau1 = 68.2e-6
tau2 = 0.405e-6

y_values = npImpulse(t_values, V_0=1.0, alpha=1/tau1, beta=1/tau2)

# plot with various axes scales
plt.figure(1)
plt.plot(t_values, y_values)
plt.grid(True)
plt.show()

# FFT calculation
sp = np.fft.rfft(y_values) # real FFT, + freq only
freq = np.fft.rfftfreq(y_values.size, t_values[1]-t_values[0])

# Normalization
# This is where things are not right.  Should be dividing by n_samples, not samples/sec
sp = sp / n_samples

# Calculation of Fourier Transform analytically
sp_analytic_real = tau1/(1+(2*pi*freq*tau1)**2)-tau2/(1+(2*pi*freq*tau2)**2)
sp_analytic_imag = -(2*pi*freq*(tau1**2)/(1+(2*pi*freq*tau1)**2)-2*pi*freq*(tau2**2)/(1+(2*pi*freq*tau2)**2))
sp_analytic_mag = np.sqrt(sp_analytic_real**2 + sp_analytic_imag**2)

plt.figure(2)
plt.plot(freq, np.abs(sp), label='FFT')
plt.plot(freq, sp_analytic_mag, label='Analytic')
plt.xscale("log")
plt.yscale("log")
plt.legend()
plt.grid(True)
plt.show()

If I set my calculation time period to 1s, the two curves match.  If I divide the magnitude of the DFT return by rfft by the sampling rate (or alternately multiply by the 1 / time period), they match.  I can't for the life of me figure out why that should be.  Can someone clue me in to the obvious that I am missing?

That's the case. Assuming you have a band-limited continuous-time signal $x_c(t)$, then the relation between its (analytically derived) CTFT, $X_c(j\Omega)$, and the (numerically computed) DTFT, $X_d(e^{j\omega})$, (or actually its DFT $X_d[k]$ ) of the samples, $x_d[n] = x_c(n T_s)$, of it taken at a period $T_s$ will be :

$$ X_d(e^{j\omega}) = \frac{1}{T_s} X_c( \frac{\omega}{T_s}) ~~~, |\omega| < \pi $$

In practice you are computing the N-point DFT $X[k]$ of $x_d[n]$ which is the uniform samples of $X_d(e^{j \omega})|_{\omega = \frac{2\pi}{N} k}$, for $k=0,1,...,N-1$. And they have the same magnitude relation too.

$$ X[k] =X_d(e^{j\frac{2\pi}{N} k}) = \frac{1}{T_s} X_c( \frac{\omega}{T_s}) ~~~, |\omega| < \pi $$

In other words, you will have a weight of
$$ K = \frac{1}{T_s}$$

between the computed DFT / DTFT and the theoretical CTFT of the signal.


Significance of $ \lambda $ in Basis Pursuit In basis Pursuit, L1 minimization is done to perform compressed sensing. In the literature there is a $ \lambda $ parameter used as a regularizer.

What is its significance?

There are 2 forms of the Basis Pursuit problem:

$$\begin{align*}
\text{The $ \lambda $ Form:} & \quad && \arg \min_{x} \frac{1}{2} {\left\| A x - b \right\|}_{2}^{2} + \lambda {\left\| x \right\|}_{1} \\
\text{The $ \epsilon$ Form:} && \quad & \arg \min_{x} {\left\| x \right\|}_{1} \\
& \text{subject to} &&  \frac{1}{2} {\left\| A x - b \right\|}_{2}^{2} \leq \epsilon
\end{align*}$$

Sometimes it is more intuitive to have a look on the $ \epsilon $ form as it tell us exactly the trade off.
The lower $ \epsilon $ it means we require teh solution to hold the Linear Equation and we pay for that with the $ {L}_{1} $ Norm of the solution.

The nice things is that for any value $ \epsilon $ we can find $ \lambda $ for the other form in such way the generate the same solution:

$$ \forall \epsilon, \; \exists \lambda : x \left( \epsilon \right) = x \left
( \lambda \right) $$

Where $ x \left( \epsilon \right) $ and $ x \left ( \lambda \right) $ are the solution of the $ \epsilon $ form and the $ \lambda $ form respectively.

You may see an example in my answer at Cross Validated Question 291962.


Are high pass/low pass filters FIR filters? I am trying to understand the basic categories of DSP filters, and it seems that the two main groups are FIR and IIR, however I'm not sure how to categorize high pass, low pass, band pass, and band stop (the 4 basic filters I learned about in university) and I can't find material on this either even though it seems so basic. Any help with making categories out of or distinctions with these groups would be greatly appreciated!

The expressions Finite Impulse Response (FIR) and Infinite Impulse Response (IIR) refer to the time-domain response of the filter. When an impulse is fed into the FIR filter in the time domain, it's response goes to zero within a finite amount of time. The response of the IIR filter does not.
Low Pass Filter (LPF), Band Pass Filter (BPF), High Pass Filter (HPF), Band Stop Filter (BSF) and Notch Filter (NF) are expressions that refer to the frequency-domain response of the filter - which frequencies it passes, stops or attenuates.
Both domains can be used to classify a filter's behavior. You can implement a LPF as a FIR filter. In the frequency domain, it will pass low frequencies and stop high frequencies. Simultaneously, in the time domain, the impulse response of the filter will return to zero after a finite amount of time. You can also implement a LPF as an IIR filter. It will have the same behavior in the frequency domain, but its impulse response will never reach zero in the time domain.
The filter's response in the time and frequency domains complement each other.
I hope this helps.


Showing that the sum of zero-mean noise is zero. Then computing the convolution of zero-mean noise with a given function This is likely to be a quick fix for people with experience in stochastic processes.
Let $ \eta[k] $ be a sequence of Uniform noise, $ \eta \sim U([-M,M]) $. I want to test if the following is correct theoretically
\begin{equation*}
\exists B>0\ s.t. \left|\sum_{k\in\mathbb{Z}} \eta[k] \right| < B.
\end{equation*}
In the general case, I would also want to try to test if
\begin{equation*}
\exists B>0\ s.t. \left|\sum_{k\in\mathbb{Z}} \eta[k] \psi [n-k] \right| <B,
\end{equation*}
What I tried:
Clearly, for $\psi[k]=1$ the two equations above are equal.
For the first one, I tried using ergodicity, which leads to $$ \lim_{N\rightarrow\infty} \frac{1}{2N}\sum_{k=-(N-1)}^{N} \eta[k] =0,$$ but this does not imply $$ \sum_{k\in\mathbb{Z}} \eta[k] =0.$$
For the second one, it's easy to show that $$\left|\sum_{k\in\mathbb{Z}} \eta[k] \psi [n-k] \right|<M\sum_{k\in\mathbb{Z}} \left| \psi [n-k] \right|, $$ but I don't want to assume that $$\sum_{k\in\mathbb{Z}} \left| \psi [n-k] \right|<\infty.$$
Even if the reader might not have an answer, it would be very helpful to know if this is a trivial result in random processes I am not aware of, or if it is a more complex problem. Thank you!

Strictly: No such $B$ exists. You could simply have a sufficient streak of "bad luck" and draw positive $\eta>\epsilon>0$ continuously, for example. Obviously, $B<\lim_{N\to\infty}\sum_{n=-N}^{N-1}\epsilon_n<\lim_{N\to\infty}\left\lvert\sum_{n=-N}^{N-1}\eta_n\right\rvert\,\forall B\in \mathbb R$.
Granted, the event that every $\eta_n>0$ has probability 0 ($=\lim_{N\to\infty}\prod_{n=-N}^{N-1} P(\eta_n>0)=\lim_{N\to\infty}\left(\frac12\right)^N$). But you asked for the existence of such $B$, not for a proof of you being able to be sure that you'll be able to pick a finite $B$ which your sum-absolute never exceeds.
Sadly, there exist uncountably infinitely many cases for which the $B$ can't exist, and all have probability 0. Whether or not that makes a non-zero sum probability isn't trivial to say, far as I can tell.
What you could show is that there's a 100% probability that your sum is below some $B$; i.e. show that:
\begin{align}
1&=P\left(\lim_{N\to\infty}\left\lvert\sum_{n=-N}^{N-1}\eta_n\right\rvert<B\right)\\
&=P\left(\lim_{N\to\infty}\left\lvert\underbrace{\sum_{n=-N}^{N-1}\frac{\eta_n}{2N}}_{=:X_N\sim\mathcal N(0,\sigma_\eta^2)}\right\rvert<\frac{B}{2N} \right)\\
&=P\left(\lim_{N\to\infty}\left\lvert X_N \right\rvert<\frac{B}{2N} \right)\\
&=P\left(\lim_{N\to\infty}X_N <\frac{B}{2N} \right)-P\left(\lim_{N\to\infty}X_N <-\frac{B}{2N} \right)\\
&=F_{\mathcal N(0,\sigma_\eta^2)}\left(\lim_{N\to\infty}\frac{B}{2N}\right) - F_{\mathcal N(0,\sigma_\eta^2)}\left(\lim_{N\to\infty}-\frac{B}{2N}\right)\\
&=\frac12-\frac12\\
&=0\text,
\end{align}
which really isn't a probability of 1. Quite the opposite is true: You can be certain that if you let a random walk go on for infinity, that you're infinitely far from where you started.
By the way, the above proof doesn't use the fact that $\eta$ is uniform – it just "needs" that all $\eta_n$ are i.i.d. and that they have bounded variance.
You could construct the very same proof even without knowledge of the central limit theorem, just using Chebyshev's inequality.


Transfer function of open loop system Suppose i have the system equation

$Y(s) = G(s)X(s)+ 3T(s)$

Then what is the transfer function of the system?

I know that the transfer function is $Y(s)/X(s)$, but i can't get that expression.

If you consider $X(s)$ as the only input of the system then the transfer function only exists for $T(s)=0$, and it is given by $G(s)$. Strictly speaking, with $T(s)\neq 0$, the system cannot be characterized by a transfer function because it is no longer linear. Note that a linear system must produce an output $Y(s)=0$ for $X(s)=0$. However, this is not the case for $T(s)\neq 0$.

You could also treat that system as a (linear) multiple input single output (MISO) system with two inputs $X(s)$ and $T(s)$. In that case there are two transfer functions:
$$H_X(s)=\frac{Y(s)}{X(s)}\Big|_{T(s)=0}=G(s)$$

and $$H_T(s)=\frac{Y(s)}{T(s)}\Big|_{X(s)=0}=3$$


Why must the order of a Band-Pass and Notch filter always be even? My professor mentioned that the order of a band-pass and a notch filter must always be even, when showing an example of designing a digital filter using the bilinear transformation.
Then he also mentions that the MATLAB code ellipord returns half the actual order, so I figured the two statements must be connected somehow. I can't ask him since it's currently 3 AM and I couldn't find anything alone as to why these two statements are true.
Any thoughts? Thank you!

In principle there is no reason why the filter order of a general bandpass or bandstop filter must be even. Such a restriction is a consequence of a specific design procedure. In classic IIR filter design (Butterworth, Chebyshev, Cauer) you start with an analog prototype lowpass filter. Bandpass or bandstop filters are then obtained by a frequency transformation. And it is this frequency transformation that doubles the order of the prototype lowpass filter, hence the even filter order for that specific design method. Note that the bilinear transform has nothing to do with that restriction on the filter order.
Of course, there are specific designs for which an odd filter order doesn't make much sense. E.g., for a notch filter (with a notch frequency greater than zero) you want exactly one zero at the notch frequency plus its mirror image at the negative notch frequency. So for each notch frequency you get two zeros, and, consequently, two poles.


why should I scale the fft using 1/N? I am writing a report, and my advisor asked me to explain why I scale the fft by a factor 1/N (where N is the length of the array).

I used to use the scaling convention of multipling the fft by the time increment (dt), this convention was good for me, because it ensures the check of the Parseval theorem. Unfortunately I had discussion with one of my advisors, because, since this convention does not give you the right amplitude, he thinks it is wrong. As I have read online, there are not right convention or wrong convention. If I use the factor 1/N the amplitude is checked, if I use the factor dt then the parseval identity is preserved.
I have 2 question now:

Why, doing the fft, I cannot have both: amplitude checked and energy checked?
I have already demonstrated (in my report) that if I scale the fft by a factor 1/N, I obtain the right amplitude, since the first value of my fft is equal to the time average of my function. Now I would like to show with formulas why this convention gives me the right amplitude? I have already searched online and here on the forum, but I did not find any good answer that explain every passage.

OK, let us go for a 2-point DFT. It should be noted that, depending on the software used, scaling can be different, and ought to be checked.
The standard unscaled version does multiply the input vector by:

$$\left[\begin{matrix} 1 & 1\\1 & -1\end{matrix}\right]$$

We have two more options: preserve the average, thus we need:

$$\frac{1}{{2}}\left[\begin{matrix} 1 & 1\\1 & -1\end{matrix}\right]$$

as $\frac{1}{{2}}\left[\begin{matrix} 1 & 1\\1 & -1\end{matrix}\right]\left[\begin{matrix} x_1\\x_2\end{matrix}\right]$ has the same average as $\left[\begin{matrix} x_1\\x_2\end{matrix}\right]$, or

$$\frac{1}{\sqrt{2}}\left[\begin{matrix} 1 & 1\\1 & -1\end{matrix}\right]$$

to preserve energy (thus orthogonality) in both domains. More generally, you have three standard options:

avoid scaling
scale by $\frac{1}{\sqrt{N}}$
or by $\frac{1}{{N}}$,

depending on the invariance purpose. But you can scale columns with other quantities, as needed for your computational purpose.
The following code aims at showing that the $1/\sqrt{N}$ scaling could work (if I am not too tired) to preserve Parseval-Plancherel or Rayleigh's energy theorem.

nSample = 2^1 ; % One can change the 2^* to other powers
data = randn(nSample,1) ;
dataFFT = fft(data)/sqrt(nSample);
ratioParseval = (norm(data)-norm(dataFFT))/norm(data)


How to validate an estimated model in case of output-only data (in frequency domain)? Before moving to the actual question, I would like to emphasize on the following points (maybe they are obvious to some of you, but I still would like to list them, since they make the difference):

I don't have any reliable information on the input, thus I have only
the measurements of the output $y(t)$.
I have only 1 set of recorded data, thus no repetitive measurements are in this case available.

Objective: Identifying a reliable transfer function model for this data in the frequency domain (e.g. rational fraction polynomial with numerator $B$ and denominator $A$)

Question: After obtaining a transfer function model, how does one validate this model and thus can be assured that it does not contain any systematic/modelling errors? Please, keep in mind the above 2 points!

To my limited experience, calculating the uncertainty of the transfer function residuals will be not possible, since the missing information on the input (and thus also not having the variance of the input). A second approach that I thought of is to do a 'simple' whiteness test of the residuals (sum of squares due to errors (SSE) or the R-square both mentioned in the Matlab manual), but is this approach even significant? I mean, sure... One obtains a result, but it's complex valued due to the data ($Y(k)$) which is also complex valued. Consequently, there is always an imaginary part which makes it fail the criterion (i.e. model is ok if SSE is close to 0). Any thoughts/solutions on this? Or am I mixing things up?

Any help, advise/guidance (tips) or references that I can consult in order to solve this issue is very much appreciated! Because I'm really lost and totally confused.

I don't know of a standard way to do this, but I think taking a leaf out of the training of artificial intelligence algorithms might help:

divide your output data into two disjoint sets,
perform whatever identification procedure you want to obtain the ARMA parameters, and
then compare the two transfer functions using whatever distance metric makes sense for your application.

In the figure below, I've decided that my distance metric is that the poles and zeros of each ARMA model are "not too far" from each other (or the ARMA poles and zeros found using all the data).

The $\color{black}{\bf black} $ markers are the pole ($\times$) and zero ($\circ$) positions for the ARMA model calculated from first half of the data ($G_1$).
The $\color{red}{\bf red} $ markers are the pole ($\times$) and zero ($\circ$) positions for the ARMA model calculated from second half of the data ($G_2$).
The $\color{grey}{\bf grey} $ markers are the pole ($\times$) and zero ($\circ$) positions for the ARMA model calculated from all of the data ($G_{\rm tot}$).

Another way to do it might be to look at the transfer function
$$
E = G_1 - G_2
$$
That is plotted below using the freqz command in R.

R Code Below

#30332
library(signal)
library(tseries)

# Create a time series that is not just white noise
T <- 1000

do_plot <- 1

#for (k in seq(1,100))
{
  x <- rnorm(T, 0, 1)
  bf <- butter(4,0.41)
  y <- filter(bf$b, bf$a, x)

  # Split data in two and compare ARIMA models
  L2 <- length(y)/2
  y1 <- y[100:L2] # Start away from 1 to avoid startup transients
  y2 <- y[(L2+1):length(y)]

  orders <- c(3,0,3)

  # First half
  tryCatch(
  {arima1 <- arima(y1, order = orders)},
  warning = function(w) {  },
  error = function(e) {  },
  finally = {});

  all_1 <- coefficients(arima1)
  num_1 <- all_1[4:6]
  den_1 <- all_1[1:3]

  # Second half
  tryCatch(
  {arima2 <- arima(y2, order = orders)},
  warning = function(w) {  },
  error = function(e) {  },
  finally = {});

  all_2 <- coefficients(arima2)
  num_2 <- all_2[4:6]
  den_2 <- all_2[1:3]

  # All data
  arima_tot <- arima(y, order = orders)
  all_tot <- coefficients(arima_tot)
  num_tot <- all_tot[4:6]
  den_tot <- all_tot[1:3]

  # Plot  boundaries
  xs <- c(-1,1)
  ys <- xs

  if (do_plot == 1)
  {
    par(mfrow=c(1,1), pty='s')
    plot(roots(den_1), xlim=xs, ylim=xs, lwd = 6, pch=4, xlab="Re(z)", ylab="Im(z)")
    do_plot <- 0
  }
  else
  {
    points(roots(den_1), lwd = 6, pch=4)
  }
  points(roots(den_2), lwd = 4, col="red", pch=4)
  points(roots(den_tot), lwd = 2, col="grey", pch=4)

  points(roots(num_1), lwd = 6, pch=1)
  points(roots(num_2), lwd = 4, col="red", pch=1)
  points(roots(num_tot), lwd = 2, col="grey", pch=1)
}

par(pty="m")
subtract_num <- conv(num_1, den_2) - conv(num_2, den_1)
subtract_den <- conv(den_1, den_2)

freqz(subtract_num, subtract_den)


Questions on Weighted Least Square Estimation This is a page from the book linear algebra,geodesy and gps by Gilbert Strang....
the page explains about the justification of the inverse of the of the co variance matrix of measurement vector $b$ in the over determined system $Ax = b$ as the best weight matrix for best estimate of x.

what does the line errors contained in the matrix mean?(underlined in blue color..)
how does that substitution take place marked with a blue curve...
is E{$rr^T$} = E{$bb^T$} ??

any explanation is welcome...

The linear model is $b=Ax+r$, so you can see the errors $r$ as the noise added to the measurements $b$.
The covariance $\Sigma_b=E\{(b-\bar{b})(b-\bar{b})^T\}=E\{(b-Ax)(b-Ax)^T\}=E\{rr^T\}$
$E\{rr^T\} \neq E\{bb^T\}$ unless $Ax=0$


How do you calculate the gain of a matched filter? I know that you can calculate the gain of a low-pass filter by summing the filter coefficients.  I don't know how to calculate the gain of an arbitrary matched filter.

For instance, I have created a matched filter for an 802.11a long training symbol (for the purposes of this question that detail doesn't matter- just thought some of you might be curious).  My MATLAB code for creating the filter is:

longSymFreq = [0 1 -1 -1 1 1 -1 1 -1 1 -1 -1 -1 -1 -1 1 1 -1 -1 1 -1 1 -1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 ...
           1 1 -1 -1 1 1 -1 1 -1 1 1 1 1 1 1 -1 -1 1 1 -1 1 -1 1 1 1 1];
longSymTime = ifft(longSymFreq);
intentionalTimeOffset = 4;
% Barrel shift by intentionalTimeOffset to avoid the corruption of the early data samples.  Due to the
% cyclic prefix doing this is not a problem.
longSymTime = [longSymTime(end-(intentionalTimeOffset-1):end) longSymTime(1:end-intentionalTimeOffset)];
longSymFilter = conj(fliplr(longSymTime));

The sum of this filter is 6.938893903907228e-018 +3.122502256758253e-017i. The average of the absolute value of the taps, though, is  0.107538417337840, so clearly the sum is not giving me any information on what kind of gain I will get when the filter gets the matched signal.

How do I calculate that gain?  The reason I care is because my final filter will be fixed-point, and I can't have an overflow when I get the signal that I'm looking for.
Is the gain the energy of the filter?

There is no single "gain" value for any filter. The value that you're referencing as the sum of the taps is its DC gain:

$$
H(0) = \sum_{n=0}^{N-1} h[n] e^0 = \sum_{n=0}^{N-1} h[n]
$$

where $H(\omega)$ is the frequency response of the filter at angular frequency $\omega \in [-\pi, \pi)$. What you are concerned about is how large the peak can possibly be when your matched filter lines up perfectly with your signal of interest. That's easy to calculate for a constant-envelope signal like yours:

$$
P = G \sum_{n=0}^{N-1} |h[n]|^2
$$

If the matched filter lines up perfectly with the signal you're looking for (both in time and in phase, which is the worst case for this analysis), you're effectively calculating the squared magnitude of each sample and summing them all up. Now, in any practical system, there will be some gain or loss factor $G$ that is not known a priori; this reflects the fact that you won't know the exact power level of the input signal as it hits your ADC. If you have an automatic-gain-control (AGC) loop, you could use its nominal output level to decide what $G$ is and therefore calculate the maximum peak $P$. Otherwise, you can design for worst-case values of $G$ to ensure that you don't overflow with the resulting value of $P$.


Add pathloss to signal I simulate a transmission channel in Matlab. I want to add pathloss ($PL=100\textrm{ dB}$) to the signal TX.
here is the code edited

PL=100;
snrdb=10;
snr=10^(snrdb/10);

L=length(TX1);
REC1=TX1+sqrt(1/(2*snr))*(randn(1,L)+1i*randn(1,L));

can anyone help regarding the pathloss effect!

To simulate 100dB pathloss from the transmitter to the receiver, you can do

tx = .....; % generate the TX signal
PL_db = 100; PL_lin = 10^(PL_db/10);
rx = tx / sqrt(PL_lin);  % the amplitude is scaled with the square root of the power scaling.

% add some noise at the receiver, according to the SNR


Is deep learning the only way to detect humans in a picture? I'm looking for a way to detect humans in a picture. For instance, regarding the picture below, I'd like to coarsely determine how many people are in the scene. I must be able to detect both standing and sitting people. I do not mind not detecting people located behind a physical object (such as the glass in the bus picture).

AFAIK, such a problem can rather easily be solved by training deep neural networks. However, my coworkers would like me to also implement a detection technique based on general image processing techniques. I've spent several days looking for techniques designed by researchers but I couldn't find anything else than saliency-based techniques (which may be fine, but I'd like to test several techniques based on old-fashioned image processing).

I have strong constraints on the nature of images, for instance luminosity: whether they've been taken during the day or at night, the pictures will have very different characteristics (e.g. image contrast can strongly vary).

Long story short, I'd like to know if there's indeed some existing segmentation technique for which it'd be interesting giving a shot, given the fact that the information carried by the images vary a lot?

One of the popular feature descriptors used for human detection is HOG - Histogram of Oriented Gradients. Usually you would train a classifier for recognizing human vs non-human, and then you would implement a sliding window technique. HOG has implementations in most of the scientific libraries in the Python community, for example in skimage there is an example here, which I am copying here because a link could always become broken:

import matplotlib.pyplot as plt

from skimage.feature import hog
from skimage import data, exposure

image = data.astronaut()

fd, hog_image = hog(image, orientations=8, pixels_per_cell=(16, 16),
                    cells_per_block=(1, 1), visualize=True, multichannel=True)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4), sharex=True, sharey=True)

ax1.axis('off')
ax1.imshow(image, cmap=plt.cm.gray)
ax1.set_title('Input image')

# Rescale histogram for better display
hog_image_rescaled = exposure.rescale_intensity(hog_image, in_range=(0, 10))

ax2.axis('off')
ax2.imshow(hog_image_rescaled, cmap=plt.cm.gray)
ax2.set_title('Histogram of Oriented Gradients')
plt.show()


Loudness of PCM stream I'm receiving PCM samples through stream which is one channel, 8000 samples per second. I figured out way of finding out loudness by calculating 'force' of signal by adding abs value of every next 1000 samples and then dividing by samples count which is 1000. It works somewhat good but sometimes my calculations are going crazy, my sum of absolute values of samples seems to constantly build up and reach ridiculous values... so I think that I'm doing something wrong there. Funny thing is that when calculations are going crazy everything will go back to normal when I hit 'sound input device' with my hand ;)

Anyway. My question is: how to properly calculate loudness of set of PCM samples?

To measure the energy (which is closely related to, but not the same as, "loudness"), calculate the RMS (Root-Mean-Square).

$ E = \sqrt{\frac{\displaystyle\sum\limits_{n=0}^{N-1} s[n]^2}{N}}$

$N$ is the number of samples and $s[n]$ is the sample at time $n$.

You can do this in any block size of samples.  You can do it on all the samples at once to get an average energy, or you can break your sample set into chunks to see how the energy changes with time.

There are also more sophisticated methods that take into account the body and brain's quirks.  You can find more information about that here.


How can I get rid of this unexpected minus sign on my inverse Fourier transform of two impulse functions? I'm trying to find the inverse Fourier transform of two impulse functions, which correspond to the Fourier transform of the function $h(t)=A\sin(2πf_0t)$.
The Fourier transform of the above sine function is as follows:
$$H(f)=j \frac{A}{2} \delta\left(f+f_0\right)-j \frac{A}{2} \delta\left(f-f_0\right)$$
When I try to find the inverse Fourier transform of $H(f)$, I keep getting a minus sign which is not in the original time-domain function. Here is my procedure:
\begin{align}
h(t)&=j \frac{A}{2} \int_{-\infty}^\infty\delta(f+f_0 ) e^{j2πft} df-j \frac{A}{2}\int_{-\infty}^\infty\delta(f-f_0 ) e^{j2πft}df\\
&=j\frac{A}{2}e^{j2πf_0t}-j\frac{A}{2}e^{-j2πf_0t}\\
&=jA\left(\frac{e^{j2πf_0t}-e^{-j2πf_0t}}{2}\right)\\
&=-A\left(\frac{e^{j2πf_0t}-e^{-j2πf_0t}}{2j}\right)\\
&=-A\sin(2\pi f_0t)
\end{align}
THAT is the minus sign that, so far, I've been unable to get rid of!
What is it that I am missing here? I hope you can help me!
Thank you very much!

HINT: you should review the sifting property of the Dirac delta impulse:
$$\int_{-\infty}^{\infty}f(t)\delta(t-t_0)dt=f(t_0)\tag{1}$$
if $f(t)$ is continuous at $t=t_0$.


Is it possible to extract peak locations in the time domain using help from fourier/wavelet analysis? The signal I'm studying has fundamental frequencies of 20 and 60 cycles per minute (shown in the Periodogram graph). It is straight forward to extract the peaks in the time domain belonging to the 20 cycles/minute frequency (both are circled), however the peaks associated with the 60 cycles per minute frequency can be noisy.

6 seconds of the time domain signal

Both signals in the time domain will ALWAYS overlap so low pass filtering isn't completely necessary but I'm just curious to know if there's a better way of extracting the 60 cycles per minute peaks (aside from peak-to-trough analysis) using frequency domain information?

Thanks!

You can figure where the peaks would be if the other tone wasn't there by an accurate estimate of the phase.  This location will be a rougher estimate in the presence of the other tone.  When you add a slanted signal (say a line) to a mode (a peak), it will shift the location of the peak somewhat.  Since you are interested in the higher frequency tone, the lower frequency one is going to look like a broad hill in comparison which is going to be less interfering than if you were doing it the other way around.
(I'm using the word "tone" broadly, meaning a pure sinusoidal signal.)


Where can I find a table of $z$-domain coefficients for Butterworth filters? The primary source lists Butterworth polynomials in s-domain and provides a link to bilinear transform for digital implementation. But, who needs analog specification in our digital world? Why should everybody get involved into the quite laborious bilinear transform, which, additionally requires frequency pre-warping? Meantime, looking at example bilinear transform in Butterworth design, I see that sampling frequency, T, cancels out, leaving you with a nice polynomial, after the monstrous transform. Why nobody can just conduct those monstrous computations only once and publish the z-polynomials, likewise s-polynomials are provided? Where can I find a table of z-polynomials not to be bogged down with the terrible pre-warping technicalities?

You can't find a table of $z$-domain coefficients that is analogous to the tables of Laplace-domain Butterworth transfer functions because it's impossible. For a second-order Butterworth filter, for example, there is no single corresponding $z$-domain transfer function, in contrast to the continuous-time case.

The reason: discrete-time transfer functions don't cover an infinitely-long frequency axis like a continuous-time system does. When creating a digital approximation to an analog prototype, techniques like the bilinear transform are used to map the frequency response from the $s$-domain transfer function evaluated along the (infinitely long) $j\omega$ axis to the $z$-domain's (finite-length) unit circle. This obviously can't be done exactly, so the various methods for mapping an analog filter to a digital one are approximations.

One factor that plays into the quality of the approximation is the sample rate of the target digital filter. In the paper you linked, you noted that the sampling interval $T$ cancelled out in the calculations of the digital filter coefficients, which is true. However, the sample interval is still required in order to perform the (necessary) frequency prewarping. Recall the prewarping formula:

$$
\omega_p = \frac{2}{T}\tan{\left(\frac{\omega T}{2}\right)}
$$

While the rest of the $T$'s may cancel out after performing the substitution $s = \frac{2}{T}\frac{z-1}{z+1}$, the $T$ in the argument to the tangent function cannot go away. Therefore, the resulting coefficients will always be a function of the sample rate $T$, so your hope of creating a simple table of digital Butterworth filter coefficients invariant of sample rate is stymied.

A couple other notes regarding this:

The prewarping step is essential to ensure that particular frequencies of interest (i.e. the cutoff frequency for a Butterworth filter) are preserved through the approximation process. Without prewarping, then the stretching of the frequency response that occurs when mapping from an infinitely-long frequency axis to the unit circle would cause the cutoff to appear at a different frequency than desired in the digital filter.
The cookbook linked in niaren's answer also demonstrates the dependence upon the sample rate: you have to trudge through the equations a bit, but you'll see factors in the listed coefficients like $\cos(\omega_0)$ and so on (this is really just a way of slightly hiding the prewarping process). The cookbook defines $\omega_0$ as:
$$
\omega_0 = \frac{2\pi f_0}{F_s}
$$
so when it all comes down to it, you'll have to define your sample rate and prewarp any frequencies of interest in order to design a Butterworth filter approximation. It's not particularly difficult, as the cookbook shows; I'm not sure what the motivation for avoiding those steps is.


DFT for different data lengths I have a sinusoidal signal that is sampled at $f_s=10$ Hz.
I am asked to find the absolute value of the DFT for different data lengths (1 sec., 2 sec. etc.)

I do not understand where will I use the information about the data lengths.

Does the question want me to change the time vector for each attempt?
as such:

** original time vector **

f=10;         % sampling freq
t=0:1/f:5;

rendering

x = cos(2*pi*3*t)

** time vector for each attempt **

t1=0:1/f:f*1; % data length = 1 sec
t2=0:1/f:f*2; % data length = 2 sec

thus rendering

x1 = cos(2*pi*3*t1);
x2 = cos(2*pi*3*t2);

and the rest would be how we usually get the FFT.

Did I understand the request correctly or not?

----Edit ------

The question is specifically for matlab simulation:

The information about the frequency of the sine seems to be missing. This could be an hint for you to compute the DFT exactly (not with Matlab simulation), for a given frequency $f$ and a number $s$ of seconds. Since the sampling frequency $f_0$ gives a integer multiple of samples inside one second, you should be able to find manageable formulas.

As commented by @Juancho, there is a lack of homogeneity in your time base: a frequency times a duration in seconds is not homogeneous to a time (it is unitless). A flexible solution is:

fs=10; % sampling freq
lFreq = [3 3.5 4]; % list of frequencies
lDuration = [1 2 4 10]'; % list of durations

for iDuration = 1:length(lDuration)
    timeBase = [0:1/fs:lDuration(iDuration)]';
    dataCosine = cos(timeBase*lFreq*2*pi); % Three cosines
    data = sum(dataCosine,2); % sum the three cosine components
    figure(iDuration)
    plot(timeBase,data,'x-');
    ylabel('Amplitude (a.u.)')
    xlabel('Time (seconds)')
    title(['Duration: ',num2str(lDuration(iDuration)),' s'])
    grid on
end


Does the Fast Wavelet Transform produce the same coefficient as the Discrete Wavelet Transform? Does the Fast Wavelet Transform(FWT) produce the same coefficients as the Discrete Wavelet Transform(DWT) if configured for the same depths? Or is the the FWT just an approximation of the DWT?

If  the discrete wavelet transform can be implemented with a FIR filter bank, with appropriate extensions, yes, up to numerical precision, coefficients will be the same.
If the discrete wavelet transform possesses a non finite support, then a FIR filter bank implementation would require filter truncation, and the results may differ. On those case, like for spline wavelets, the processing is often performed in the Fourier domain.
There exist IIR discrete wavelets, which may bridge this gap, yet I am not familiar enough to tell.


adding noise to a signal I have an audio file. I took an extract of the file called variable "CALL", and I took another extract from the file which is just filled with background noise, called "BACKGROUND".
Using matlab, I found the signal to noise ratio of CALL using:
snr(CALL,BACKGROUND); % both variables are of the same length

the signal to noise outputted was     5.3949

I want to add noise to CALL, and lower this signal to noise to 3 for example. How would I go about programming this?

I've only seen examples where people add noise to a perfect signal (i.e. they just create some kind of sinewave, then add noise to that), I have a signal which already has noise, and I want to add more noise, and know its signal to noise value.

Also, is the above method for how I found my first signal to noise value correct? (I'm assuming the background noise is somewhat constant throughout the audio file)

maybe the info at the following web page will be of some use to you.  http://www.dsprelated.com/showcode/263.php  [-Rick-]


Does changing the frequency of a wave cause a change of the phase of the same wave? I know that the phase is the delay between two different waves . But if we change the frequency of two waves by a different value does it cause a varying in the phase between them let say that two waves start from the same point.and if I apply that in the single wave and change the frequency does that related to the phase !
My main question related to the Fm and Pm modulation

Let's assume your two sinusoids are at the same phase, $\phi_0$, at $t=0$ and you change the frequency of one of them at exactly that time.

The phases as functions of time, for $t \ge 0$, for those sinusoids at the two different frequencies are

$$\phi_1(t) = 2\pi f_1 t + \phi_0$$
$$\phi_2(t) = 2\pi f_2 t + \phi_0$$

And the difference between those phases is

$$\psi(t) = \phi_1(t)-\phi_2(t) = 2\pi(f_1-f_2)t$$

So the phase difference between them is a function of time and the difference between the two frequencies.


MFCC Deltas, Delta-Deltas I have just a quick question about how to calculate Deltas from MFCC. I have found the formula for calculating deltas here:

$$d_t=\frac{\sum_{n=1}^{N}n\left(c_{t+n}-c_{t-n}\right)}{2\sum_{n=1}^{N}n^2}$$

Another source states it is a regression equation. From what I can see, it takes coefficient from a previous and next frame. For any $N$ e.g. $N=2$, $t-n$ in $c_{t-n}$ can refer to the past where there are no frames. Which values should I substitute in this case? Should it be simply 0?

Similarly, when $t+n$ in $c_{t+n}$ refers to a frame in future that not exists, should it also be 0?

Thanks

It's better to copy first frame and last frame values to extend vector sequence beyond boundaries than to assign 0. This could be implemented just by adjusting indexes:

if (index1 < 0)
    index1 = 0
if (index2 > N - 1)
   index2 = N - 1

delta = v[index1] - v[index2]


Discrete Fourier transform of an exponential decay I have a vector with an exponential decay, using Numpy:

t=np.arange(128)
k=0.1
decay=np.exp(-k*t)

I would like to compute the discrete Fourier transform (DFT) of decay so I get the same result as applying np.fft.rfft(decay, n=128*2). I tried the formula described here plus some similar ones (without square terms in the denominator) but I never got the same result. The final result I want seems to be a one-sided Lorentzian? The reason I want to do this is because computing a Lorentzian is faster than computing an exponential and then applying fft.

To summarize, is there an easy way to compute the DFT for a given exponential decay function of known k?

This boils down to computing the sum of a geometric series. Given

$$x[n]=e^{-\alpha n},\qquad n=0,1,\ldots,N-1\tag{1}$$

the DFT is

$$\begin{align}X[k]&=\sum_{n=0}^{N-1}x[n]e^{-j2\pi nk/N}\\&=\sum_{n=0}^{N-1}e^{-\alpha n}e^{-j2\pi nk/N}\\&=\frac{1-e^{-\alpha N}}{1-e^{-\alpha}e^{-j2\pi k/N}},\qquad\alpha\neq 0,\quad k=0,1,\ldots N-1\tag{2}\end{align}$$

For $\alpha=0$, you have $X[0]=N$ and $X[k]=0$ for $k=1,2,\ldots,N-1$.

If you zero pad with $N$ zeros, i.e., if the DFT length is $2N$, the result can be computed in a very similar manner:

$$\begin{align}X[k]&=\sum_{n=0}^{N-1}x[n]e^{-j\pi nk/N}\\&=\sum_{n=0}^{N-1}e^{-\alpha n}e^{-j\pi nk/N}\\&=\frac{1-e^{-jk\pi}e^{-\alpha N}}{1-e^{-\alpha}e^{-j\pi k/N}}\\&=\frac{1-(-1)^ke^{-\alpha N}}{1-e^{-\alpha}e^{-j\pi k/N}},\qquad\alpha\neq 0,\quad k=0,1,\ldots 2N-1\tag{3}\end{align}$$

For $\alpha=0$ you get

$$X[k]=\begin{cases}N,\quad k=0\\0,\quad k\textrm{ even},k\neq 0\\\displaystyle\frac{2}{1-e^{-j\pi k/N}},
\quad k\textrm{ odd}\end{cases}\tag{4}$$


I am looking for an analytic description of a continuous-time Butteworth High-pass filter in the time domain (=impulse response) having derived the Butterworth Lowpass Time domain response, I am now struggling to find a similar function for a Butterworth Highpass filter. I understand you need to replace s by 1/s. But this leads to a transfer function with as many zeros as there are poles. So you can no longer use the Heaviside 'cover-up' method to decompose the individual exponential functions, to arrive at the time domain response.
Is there perhaps someone who stumbled on the same problem and arrived at a solution ?
The ultimate goal is to arrive at a Butterworth bandpass filter whereby the lowpass and high pass sides can have a different order (slope). So I still need to convolve the LP- with the HP-filter. . .

$\frac{s^n + b_{n-1}s^{n-1}+\cdots+b_0}{s^n + a_{n-1}s^{n-1}+\cdots+a_0} = 1 + \frac{(b_{n-1}-a_{n-1})s^{n-1}+\cdots+b_0 - a_0}{s^n + a_{n-1}s^{n-1}+\cdots+a_0}$, and the inverse Laplace transform of 1 is $\delta(t)$.
Is that enough to get on with?


Why is successive decomposition of a signal performed only on low pass component? The question says it all. In typical (wavelet-like) decomposition of a signal, why is only the low pass component chosen for successive decomposition ?

Wavelets decomposition separates out the details/fluctuations/high-pass information from the image or signal. At each step details are separated out from the remaining of the signal. The processing is therefore only applied on the coarse part not the part that already has the details.


Where to get the Kodak PhotoCD raw image dataset? I've been doing some background reading on demosaicing raw images.
I've seen the Kodak PhotoCD raw image dataset mentioned in a few papers on demosaicing:

Image Demosaicing: A Systematic Survey
High Quality Linear Interpolation for Demosaicing of Bayer-Patterned Images
Demosaicing using optimal recovery

After some googling, I haven't been able to find the Kodak dataset. Where can I download the Kodak PhotoCD dataset?

The dataset appears to be in the public domain, and has been posted at this site here.

The pictures (24 in total) were originally released in 1991 on a CD.  The original media with the actual image data could be requested via interlibrary loan; here is a WorldCat link to the repositories holding this item. The CD was originally distributed by Apple Computer.  The raw version of the files do not appear to be available online.


Image Reconstruction:Phase vs. Magnitude Figure 1.(c) shows the Test image reconstructed from MAGNITUDE spectrum only. We can say that the intensity values of LOW frequency pixels are comparatively more than HIGH frequency pixels.

Figure 1.(d) shows the Test image reconstructed from PHASE spectrum only. We can say that intensity values of HIGH frequency (edges,lines) pixels are comparatively more than LOW frequency pixels.

Why this magical contradiction of intensity change (or exchange) is present between Test image reconstructed from MAGNITUDE spectrum only and Test image reconstructed from PHASE spectrum only, which when combined together form the original Test image?

clc;
clear all;
close all;
i1=imread('C:\Users\Admin\Desktop\rough\Capture1.png');
i1=rgb2gray(i1);

f1=fftn(i1);
mag1=abs(f1);
s=log(1+fftshift(f1));
phase1=angle(f1);

r1=ifftshift(ifftn(mag1));
r2=ifftn(exp(1i*phase1));
figure,imshow(i1);
figure,imshow(s,[]);
figure,imshow(uint8(r1));
figure,imshow(r2,[]);
r2=histeq(r2);
r3=histeq(uint8(r2));
figure,imshow(r2);
figure,imshow(r3);

Figure 1.(c) shows the Test image reconstructed from MAGNITUDE spectrum only. We can say that the intensity values of LOW frequency pixels are comparatively more than HIGH frequency pixels.

Actually, this is not correct. The phase values determine the shift in the sinusoid components of the image. With zero phase, all the sinusoids are centred at the same location and you get a symmetric image whose structure has no real correlation with the original image at all. Being centred at the same location means that the sinusoids are a maximum at that location, and is why there is a big white patch in the middle of Figure 1.c.

The phase-only reconstruction preserve features because of the principle of phase congruency. At the location of edges and lines, most of the sinusoid components have the same phase. See http://homepages.inf.ed.ac.uk/rbf/CVonline/LOCAL_COPIES/OWENS/LECT7/node2.html  This properly alone can be used to detect lines and edges, http://www.csse.uwa.edu.au/~pk/research/pkpapers/phasecorners.pdf , without regard to magnitude. So you can see that the phase information is most important.

Changing the magnitude of the various component sinusoids changes the shape of the feature. When you do a phase-only reconstruction, you set all the magnitudes to one, which changes the shape of the features, but not their location. In many images the low frequency components have a magnitude higher than the high frequency components, so phase-only reconstruction does look like a high-pass filter.

In short, phase contains the information about the locations of features.

You cannot add the phase-only and magnitude-only images to get the original. You can multiply them in the Fourier domain and transform back to get the original.