Beamforming using MATLAB I'm attempting to beamform some data which has signals coming from all directions in a 180 degree range.

I found an example of using the phased array toolbox, in MATLAB, to beamform a signal in https://uk.mathworks.com/help/phased/examples/conventional-and-adaptive-beamformers.html however, when I added a second signal from a different direction, the results aren't what I would expect.

If the two signals are coming from directions which are symmetrical about 90 degrees, i.e. 30 and 150 or 70 and 110, the beam forming algorithm cannot separate them. Is there a way of suppressing the signal coming from the other direction?

The code below is an example of what I mean where two signals are incoming from 30 and 150 degrees and then the beamforming algorithm is trying to isolate the signal from 30 degrees.

clear;
clc;
close;
t = 0:0.001:0.3;                % Time, sampling frequency is 1kHz
s = zeros(size(t));
s = s(:);                       % Signal in column vector
s(201:205) = s(201:205) + 1;    % Define the pulse
s_2 = zeros(size(t));
s_2 = s_2(:);                       % Signal in column vector
s_2(41:47) = s_2(41:47) + 1;    % Define the pulse
figure(1)
plot(t,s, t, s_2);
title('Pulse');xlabel('Time (s)');ylabel('Amplitude (V)');
carrierFreq = 30e6;
wavelength = physconst('LightSpeed')/carrierFreq;
ula = phased.ULA('NumElements',10,'ElementSpacing',wavelength/2);
ula.Element.FrequencyRange = [90e5 110e6];
inputAngle = [30; 0];
x = collectPlaneWave(ula,s,inputAngle,carrierFreq);
inputAngle2 = [150; 0];
x2 = collectPlaneWave(ula,s_2,inputAngle2,carrierFreq);
rs = RandStream.create('mt19937ar','Seed',2008);
noisePwr = .5;   % noise power
noise = sqrt(noisePwr/2)*(randn(rs,size(x))+1i*randn(rs,size(x)));
rxSignal = x + x2 + noise;
figure()
subplot(211);
plot(t,abs(rxSignal(:,1)));axis tight;
title('Pulse at Antenna 1');xlabel('Time (s)');ylabel('Magnitude (V)');
subplot(212);
plot(t,abs(rxSignal(:,2)));axis tight;
title('Pulse at Antenna 2');xlabel('Time (s)');ylabel('Magnitude (V)');
psbeamformer = phased.PhaseShiftBeamformer('SensorArray',ula,...
                'OperatingFrequency',carrierFreq,'Direction',inputAngle,...
                'WeightsOutputPort', true);
[yCbf,w] = step(psbeamformer, rxSignal);
figure()
plot(t,abs(yCbf)); axis tight;
title('Output of Phase Shift Beamformer');
xlabel('Time (s)');ylabel('Magnitude (V)');

Is it possible to get the signal from just one direction if the signal composes of reflected energy from all directions?- And if so, how?

Thanks for reading!

In MATLAB the default direction of orientation of Planar array is y axis, which can be seen as

ula.getElementPosition

ans =

         0         0         0         0         0         0         0         0         0         0
  -22.4844  -17.4879  -12.4914   -7.4948   -2.4983    2.4983    7.4948   12.4914   17.4879   22.4844
         0         0         0         0         0         0         0         0         0         0

which makes your planar array graphical representation of antenna elements as beam pattern to be

The Pattern is identical in front and back of the planar array, so this bring us to the region of operations for your case either front of array or back of array.

In case you want to shift the beam to 30 degrees then you can not distinguish between +30 and +150 degrees because as per matlab definition of azimuth angle which can be found here we are looking at identical beams in different regions as can be seen from the figure below

I have tried to change the azimuth angles as per MATLAB convention +30 and -30 degrees and the signal is now separable. So you can work from +90 to -90 degree for Region 1

The modified code

clear;
clc;
close all;
t = 0:0.001:0.3;                % Time, sampling frequency is 1kHz
s = zeros(size(t));
s = s(:);                       % Signal in column vector
s(201:205) = s(201:205) + 1;    % Define the pulse
s_2 = zeros(size(t));
s_2 = s_2(:);                       % Signal in column vector
s_2(41:47) = s_2(41:47) + 1;    % Define the pulse
figure(1)
plot(t,s, t, s_2);
title('Pulse');xlabel('Time (s)');ylabel('Amplitude (V)');
carrierFreq = 30e6;
wavelength = physconst('LightSpeed')/carrierFreq;
ula = phased.ULA('NumElements',10,'ElementSpacing',wavelength/2);
ula.Element.FrequencyRange = [90e5 110e6];
inputAngle = [30; 0];
x = collectPlaneWave(ula,s,inputAngle,carrierFreq);
inputAngleTwo = [-30; 0];
x2 = collectPlaneWave(ula,s_2,inputAngleTwo,carrierFreq);
rs = RandStream.create('mt19937ar','Seed',2008);
noisePwr = .5;   % noise power
noise = sqrt(noisePwr/2)*(randn(rs,size(x))+1i*randn(rs,size(x)));
rxSignal = x + x2 + noise;
figure()
subplot(211);
plot(t,abs(rxSignal(:,1)));axis tight;
title('Pulse at Antenna 1');xlabel('Time (s)');ylabel('Magnitude (V)');
subplot(212);
plot(t,abs(rxSignal(:,2)));axis tight;
title('Pulse at Antenna 2');xlabel('Time (s)');ylabel('Magnitude (V)');
psbeamformerOne = phased.PhaseShiftBeamformer('SensorArray',ula,...
                'OperatingFrequency',carrierFreq,'Direction',inputAngle,...
                'WeightsOutputPort', true);
[yCbf,w] = step(psbeamformerOne, rxSignal);
figure()
subplot(2,1,1)
plot(t,abs(yCbf)); axis tight;
title(['Output of Phase Shift Beamformer when beam is at: +30 degrees' ] );
xlabel('Time (s)');ylabel('Magnitude (V)');

clear yCbf w;
psbeamformerTwo = phased.PhaseShiftBeamformer('SensorArray',ula,...
                'OperatingFrequency',carrierFreq,'Direction',inputAngleTwo,...
                'WeightsOutputPort', true);
[yCbf,w] = step(psbeamformerTwo, rxSignal);
subplot(2,1,2)
plot(t,abs(yCbf)); axis tight;
title(['Output of Phase Shift Beamformer when beam is at: -30 degrees' ] );
xlabel('Time (s)');ylabel('Magnitude (V)');

Result:

Now the two signals are separable


Given Gain function, how to design a causal, stable and minimum phase IIR filter? I am given $|H(\omega)|$, I wonder if minimum phase stable causal filter is unique and how to calculate it.

If $H(\omega)=e^{\alpha(\omega)+j\phi(\omega)}$ is a minimum phase frequency response, then the attenuation $\alpha(\omega)$ and the phase $\phi(\omega)$ are related by the following Hilbert transform relationship:

$$\phi(\omega)=-\frac{1}{\pi}\int_{-\infty}^{\infty}\frac{\alpha(\Omega)}{\omega-\Omega}d\Omega$$

So $\phi(\omega)$ is uniquely determined by $\alpha(\omega)$. This is the corresponding wikipedia entry.


Respiration (breathing) signal movement artifact removal Situation: I have breathing data of participants who have done different things during the study which caused their breathing change; watching mindless video to which resulted in unconscious regular baseline breathing; attending to a pacer with causes slow-paced breathing in range of 5 to 9 breath per minute; performing a stressor task which cause breathing to become fast and irregular.
Goal: I need to extract breath per minute (BPM) and breathing irregularity (using FFT).
Challanges. Participants' movement caused the signals look noisy . Note: this is a baseline data which supposed to be regular.  I need to remove the artifacts related to movement but I don't know how. Why is this important? My BPM values can be wrong and the FFT gives a large power to frequencies very close to zero.
File examples: http://stanford.edu/~parism/HPC_physio_split_undownsampled_256_example/
These are files from a person in the treatment group. The sequence of the columns in the file is as follows:
'Time','BVP', 'EDA', 'Temp','ABD', 'CST','Label'
The treatment group received the breathing pacer intervention during pre-, stressor, and post- block2. That means the pacer created vibration patterns at frequency of 6breath per minute (for this participant) to breathe with. During the pre- and post- stressors, the participant is not experiencing any cognitive tasks. Therefore they are very likely to synchronize their breathing with the pacer. During the the stressor2 tasks, they are less likely to do so. The pre- post- and stressor 1 blocks are similar to those in block 2. The only difference is that the participant is not receiving any intervention. Therefore they are less likely to do slow-paced breathing.
During the baseline the participant is not receiving any vibrations. They are just breathing at their resting heart rate. During the meditation phase, the participant is listening to an audio to pace to practice slow-paced breathing. We derived the pace of the pacer (in this example 6) from the last 30 seconds of this file. During the Practice, the participant is breathing with the pacer for about 90 seconds to ensure that they are comfortable with the pacer's pace.
Thank you in advance for help.

I figured out how to tackle the problem.

I used the EMD algorithm and then smoothed the signal. I removed the first two IMF[i] and the last IMF[i] and reconstructed the signal again. This removes posture movement.
I then  identified 3sd away outliers and chopped off the signal to pieces. For each piece I calculated the FFT and found the dominant frequency.


How to calculate the noise power for a non-stationary noise? With stationary noise we have constant mean and variance (let's assume it is Gaussian noise). My first question is, how is the noise power calculated and how it is related to the variance?

Now, I have a non-stationary noise which its variance is not constant and changes periodically over time. My second question is, how can I calculate the noise power for this noise and how it is related to the time-variant variance?

Any help and remark on the topic would be appreciated.

Thank you all.

If your signal is just $n(t)$ distributed as $N(\mu, \sigma^2)$, then the expected signal power is then:
$$
p_{\rm ensemble} = E\left [ \left| n(t) \right|^2 \right] = \mu^2 + \sigma^2
$$

For the time-varying case, I'd suggest that the answer is:
$$
p_{\rm time\ varying} = \mu^2(t) + \sigma^2(t)
$$


How does Delta-Sigma Modulation convert a 1-bit signal to higher resolution signal? I've looked at various resources about Delta-Sigma modulators, and I find them fascinating.  However, I'm confused on a particular point and and am looking for an explanation.

In my understanding, the Delta-Sigma modulators operates by sampling at a very high sampling rate (much higher than the ulitmate rate) an input analog signal.  During this process the digitized signal converted back to into an analog waveform and fed back into the input signal (possibly a number of times) which has the effect of high pass filtering the quantization noise -- so-called 'Noise Shaping'.

So far, so good.  To me, this makes sense.  However, the next stage, I find confusing.

The resulting digitized signal has been quantized using only single bit, and somehow, a digital filter is then applied to this signal which does two things:

(1) Filters out the shaped quantization noise and
(2) Re-quantizes the signal to a much higher resolution (e.g. 1-bit => 24-bit)

The resulting signal is then decimated to the desired sampling rate.

My question is about the re-quantization.  How is this performed?  Any pointers appreciated.  If you wish to poke any other holes in my understanding of DSM, that's also fine.  Cheers!

The output of the low-pass filter consists, in principle, of real numbers. This is because, unlike the input values, the filter coefficients are not quantized to one bit and several of them are added to yield one output sample, as well.

Of course, there is no such thing as real numbers in a finite computer; but you may compute the filter output to any desired precision, such as 24-bit.


change the central frequency of a wavelet I am using the discrete meyer wavelet in MATLAB for time series decomposition. I would like to change the central frequency of this wavelet but I do not know if this is possible or if I would need to design a new wavelet.
Thank you very much in advance.

MATLAB documentation indicates that discretization of the continuous and infinite-support Meyer wavelet is done simply by regular sampling. I doubt that there are any unwanted consequences if you sample it over a longer interval than the nominal -8 to 8 effective support, to have a different central frequency.

The standard values for the first two parameters to meyer(LB,UB,N) are LB = -8, UB = 8, defining the interval to be covered by sampling. Multiply those parameters by the same factor by which you want to multiply the central frequency. If factor < 1, you should double LB and UB to ensure that you sample the whole effective support of the Meyer wavelet and scaling functions. Then you will find your wanted frequency on the next scale level.


wide sense stationary of dynamic process I am trying to understand the definition of wide sense stationary on my own and probably have some silly questions. Wikipedia says, wide sense stationary is a process with constant mean and  autocorrelation function over time. Now, lets assume a process,
$$X(n) = 0.9X(n-1) + V(n) + U(n),$$
where $V$ is white noise and $U$ is an external input.

Lets assume $U(n) = 0\quad \forall n$, i.e. $X_1(n) = 0.9X_1(n-1)+V(n)$. In this case can $X_1$ be called WSS? What does mean being constant actually implies? Can I assume $0.9X(n-1)$ is mean since it is deterministic and say mean is decreasing over time, and this is not WSS?
Lets assume $U(n) = C \text{ const.}$, i.e. $X_2(n) = 0.9X_2(n-1) + V(n) + C$
If $X_1$ was WSS, I believe this should be WSS as well. Please confirm.
If $U(n)$ is an external signal changing over time, i.e. $X_3=X$ then would still be WSS, if $X_p= 0.9X_p(n-1)+V(n)$ were originally? In this case $X$ is depended on a signal which is not part of the original process.

Ok, let's rule out 3. first. Whether that is WSS depends on $U$ and you don't know that. For example, if $U$ has unbounded variance, then so does $X_3$, and bounded variance is a necessary condition for WSS.
For the others:
The aspect of the definition you're looking for is that you can define an autocorrelation function that only depends on the time difference:
$$\phi_{X_kX_k}(n,m)= E\left(X_k(n)X_k^*(m)\right) \overset != E\left(X_k(n)X_k^*(n-h)\right)=:\phi_{X_kX_k}(h)$$
for all times $n$ (and $m$).
Now, let's apply this definition of the autocorrelation function to your $X_1$:
\begin{align}
\phi_{X_1X_1}(n,m) &= E\left(X_1(n)X_1^*(m)\right) \\
&=E\left[\left(0.9X_1(n-1)+V(n)\right)\left(0.9X_1(m-1)+V(m)\right)^*\right]\\
&=E\left[0.81 X_1(n-1)X^*_1(m-1)\right]+E\left[X_1(n-1)V^*(m)\right]+E\left[X^*_1(m-1)V(n)\right]+E\left[V(n)V^*(m)\right]\\
&\text{due to noise being }\color{blue}{\text{independent}}\text{ from signal and } \color{gray}{\text{white}}\\
&=0.81E\left[X_1(n-1)X^*_1(m-1)\right]+\color{blue}{E\left[X_1(n-1)\right]E\left[V^*(m)\right]}+\color{blue}{E\left[X^*_1(m-1)\right]E\left[V(n)\right]}+\color{grey}{0}\\
&\text{white noise must have a }\color{green}{\text{zero mean}}\\
&=0.81E\left[X_1(n-1)X^*_1(m-1)\right]+E\left[X_1(n-1)\right]\color{green}{0}+E\left[X^*_1(m-1)\right]\color{green}{0}\\
&=0.81E\left[X_1(n-1)X^*_1(m-1)\right]\tag Ä\label{fail}
\end{align}
But as seen in \eqref{fail}, $\phi_{X_1X_1}$ depends on the actual points in time $n$ and $m$, and hence, $X_1$ can't be WSS.
Let's address

If $X_1$ was WSS, I believe this should be WSS as well. Please confirm.

Neither $X_1$ nor $X_2$ are WSS, same derivation as above, just that the expectations of noise aren't zero, but still time independent whereas the $X_2$-product isn't.


Camera Pose Estimation from Vanishing Points I want to estimate the pose of camera (Location and Orientation). Can I use Vanishing Points to get the location and orientation of the camera?

There is object in the image which I also know his size, can it help?

You can use the vanishing points to calibrate the camera.
I recommend reading  Chapter 8.6 of Multiple View Geometry in Computer Vision [ http://www.cambridge.org/9780521540513 ]

Even better, if you have those points on a known sized plane, then you can compute an homography H and the it is possible to compute both intrinsic and extrinsic camera matrices [K],[R|T] respectively.

The size of the object helps determining the scene scale, or more precisely the scale of the column vector t in the extrinsic matrix.
(Zhang's calibration method)[ https://www.microsoft.com/en-us/research/wp-content/uploads/2016/02/tr98-71.pdf ]


Fourier series of a square pulse (simple simplification) Hi guys I'm currently studying for finals and I came across a square pulse problem.

I ended up getting this:

I'm wondering if I can simplify this to a sin function. If so how? And please correct me if I did the problem wrong in the first place.

You could rewrite $x(t)$ as a rectangular function before doing the Fourier transform. You could deduce the width and center of the rect signal directly from the given expression but it would be a easier to just draw the signal and deduce the values from the graph.

Your rect signal will be in the following form:

$x(t) = rect(\frac{t-t_o}{W})$

where:

$t_0$ - center of the rectangular signal

W - width(duration) of the pulse

From this point you would use a Fourier transform table and some Fourier transform theorems (involving scaling and delay) to get the $sinc$ version of the $X(\omega)$.

Or, if you prefer, you could just modify a little bit the expression you already have to get the expression where you could use Euler's formula as Dilip Sarwate commented.

$X(\omega)=1 \cdot \frac{1}{j\omega}(e^{-j\omega} - e^{-j2\omega})= e^{-j1.5\omega}e^{+j1.5\omega} \frac{1}{j\omega}(e^{-j\omega} - e^{-j2\omega}) = e^{-j1.5\omega}\frac{1}{j\omega}(e^{+j0.5\omega} - e^{-j0.5\omega})$

From this point it should be easy to get the result.


Why does FFT generate a jaggedy signal I have an application that fills a buffer from my computer's microphone and does an FFT. My objective is to figure out what frequencies are in the sample.

I've tested the application by whistling as well as by playing this video from my cellphone into my computer's mic (I've tried two different mics to be sure). What I expect (taking the absolute frequency) is a smooth Gaussian-like curve with a mean at the frequency of the incoming tone, but what I'm getting is this:

What is the cause for the jaggedness of the signal's frequency spectrum? And what technique can be used to reliably identify the main frequencies present in the signal while filtering out the unwanted bumps?

What you are seeing is a peak from sinusoid convolved with a spectrum of rectangular window sinc function. That is due to property:

$$x(t)\star w(t) = X(f)\cdot W(f) $$

Where your signal $x(t)$ is a sinusoid with a spectrum being a Kronecker (Dirac) Delta function, $w(t)$ is a rectangular window function. In DFT by default we are using rectangular window (taking a time snippet of a signal). Please check out this paper for more info: A Comprehensive Window Tutorial


Is this a power signal? $x[n]=\sin(\pi n)$

Since this is a discrete signal, and n can take only values like ...,1,2,3 ... the value of the power is getting zero, and since it gets zero does it mean that it's not a power signal?

There a few variations on the definitions of "energy signal" and "power signal". For example, Sklar in his textook "Digital Communications" defines "energy signal" as a signal with energy $E$ such that $0 < E < \infty$, and a "power signal" as a signal with power $P$ such that $0 < P <\infty$.

However, others (see for example this answer) define "energy signal" as a signal with $E < \infty$.

In your case, $x[n]=0$ for all $n$; then this signal has zero energy and zero power. Whether it is classified as an energy or a power signal depends on which specific definition you adopt. For instance, if you go with Sklar's, then this signal is neither an energy signal nor a power signal.


FFT results for short burst of signal I am sampling short bursts of signal (20ms length) in 30ms buffers in time domain.
The 30ms buffer consists from 4096 I/Q samples from which roughly 2550 cover the actual signal. Now I do in-place FFT on that.
If carrier signal is in the sampled data then it is found by FFT but the power of frequency bin varies depending where inside of the samples the actual signal is located.
If I have more noise after signal then FFT power is higher, if I have more before signal then it is weaker.
How I can overcome this problem ? I understand FFT will interpolate signal over the noise samples in this case.
I would like to have the same FFT power for signal frequency despite of location inside the processed samples.
Thanks
Michal

Eventhough you are asking for an answer about your problem with FFT usage, It seems to me, based on your definition of what you are looking for (detecting the presence of a carrier -a known signal?- inside a buffer of 4096 samples) that you would better use a "matched-filter" to detect the presence of that carrier, a known signal, inside the noisy 30ms measurement buffer.

Of course FFT would still help especially for high SNR cases, but for noisy data, matched filter is the optimal choice (under suitable conditions). And I have never heard of FFT by itself reducing noise. May be I am wrong.

To implement a matched-filter, you would simply convolve the time reversed carrier signal with the noisy observation and look for a thresholded peak.


Signal Estimation after detection Part 2 I have a binary information source signal $s(t)$ that is corrupted by additive White Gaussian noise $w(t)$ at a particular SNR. The received signal is:
$$x(t) = s(t) + w(t)$$
Then I have created a matched filter $h(t)$ as the time reversal of the source binary signal $s(t)$.
Question: After signal detection, how to estimate the clean signal $s(t)$? I have used the sign() operation. Is that the correct way or should I use sophisticated methods like MLE, LMS etc?
In many implementations, a hard threshold zero is used to decode the signal in order to get back the transmitted source symbols $\hat{s}(t)$. Is that the correct way? Is my implementation correct where I have used sign function to estimate the symbols.
This is my implementation. Please correct me where wrong.
clear all
N = 50;

input = rand(1,N)>0.5;
s=(2*input-1); %input
x = awgn(s,15,'measured'); %received noisy signal

matched_filter_h = flipud(s);
s_hat = sign(filter(matched_filter_h,1,x));

The code output shows that $s(t)$ and $\hat{s}(t)$ are identical. So, it seems that the estimation is possible.

After signal detection, how to estimate the clean signal $s(t)$?

Matched filtering is used to detect the presence of a known signal in noise. There is no estimation part when you are talking about a matched filter. The estimate part comes after you have done the matched filter and need to estimate the symbols.
It looks like you are talking about a communication system context, but the matched filtering in those systems is not done by matched filtering the symbols as in your example code. For each symbol, the true symbol is unknown and the pulse shape is known, and that is why the matched filter uses the pulse shape to matched filter. Perhaps all that you want to know about matched filtering can be found here, https://dsp.stackexchange.com/a/9389/31316, and you should give it a read.
Edit
The OP questions about the zero threshold and whether there are other methods. That threshold is set through using the maximum likelihood solution, or minimum distance decoding. For the case of BPSK, you want to decide if $x$ contains $s_1=+1$ or $s_2=-1$, and the ML rule is to choose $\hat{s}=s_1$ if:
$$ |x-s_1|<|x-s_2| $$
This can also be interpreted as partitioning the IQ diagram into different regions (Voronoi regions), a region for $s_1$ and a different region for $s_2$. These regions end up looking like this (note how for BPSK the threshold line is at zero):

For full derivation, see Tse's Wireless Comms book section A.2.1, https://web.stanford.edu/~dntse/Chapters_PDF/Fundamentals_Wireless_Communication_AppendixA.pdf.


Capacity and channel estimation algorithms: explanation of results I study a spectral efficiency of a system with a precoding scheme with different channel schemes.
Algorithm A for channel estimation gives a smaller estimation error instead of Algorithm B.
I have computed the precoding schemes using the resulted estimated from Alg. A and Alg.B.
As a result, the capacity with precoding for Algorithm A is better ( bigger) than for Algorithm B.

How can I explain the results?
Can I say that the reason is in estimation error? ( estimate of Alg A
is better)

The channel capacity expression you cited is the one with CSI at receiver that means the receiver know perfectly the realization of fading channel, denoted $\mathbf{H}$, but not the realization of additive white Gaussian noise.
In that case, the channel capacity is the max (or sup) data rate at which there exists a receiver that achieves decoding error probability smaller than any infinitesimal $\epsilon > 0$, and can be computed using the expression in your links (with adding expected value operator depending assumptions) see chapter 8.
Therefore, I would not be convinced if people say their receiver needs to estimate the channel while it is assumed to be known (that estimation is not illegal though).
That said, it is legit if the channel estimation is not at the receiver. One such case is that the channel estimation is done at transmitter. But in any case, receiver must known perfectly the channel so that you can use the channel capacity expression.
Thus, if your transmitter estimates the channel and calculate precoding vectors to transmit data and your receiver knows perfectly the precoding and the channel, you can use the capacity expression.


Can we broadcast two messages using NOMA? Is it possible for a transmitter to broadcast two messages to a set of receivers at the same time using superposition coding (Non-Orthogonal Multiple Access - NOMA)?

A transmitter has two messages $m_1$ and $m_2$ that wants to broadcast to a set of $n$ receivers. I think this not possible even for $n=2$. If the transmitter broadcasts $P_1m_1 + P_2m_2$ with transmit power factors $P_1$ and $P_2$, then only the receiver with higher channel gain can decode both $m_1$ and $m_2$. Is this true?

There is a THRESHOLD, both can receive. If that would not have been possible why would NOMA be so popular, What is the use of NOMA then? Well, theoretically is possible however practically such channel state information is not available, that is what we read, however, I haven't seen whether that practically possible or not, Energy-Efficient NOMA Systems tells exactly what you have just said.


Frequency Axis problem in a DTFT I have a doubt related to calculating the Discrete Time Fourier Transform (DTFT) by hand. Specifically in how calculate the frequency axis of the spectrum. My signal has N values and was sampled at FS Hz, the spectrum would have N entries too (where N/2 values are a mirror of the other half). The maximum representable frequency is FS/2 (by the Nyquist theorem), that means that I have to multiply each entry by FS/N, so for the last entry (N/2) I have this:

(FS/N) * (N/2) = FS/2

But when I sample at a higher FS Hz the spectrum is shifted. Think for instance in this function:

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

Where "f0 = 1/T", T is the period of the cosine and "t" means each entry of the time axis. Then the spectrum is two pulses in "-f0" and "f0". But doing this in python:

f = range(-N/2,N/2)
f = [float(FS)/(float(N)) * i for i in f]

And sampling at a higher FS, then the pulses are shifted. But the correct behavior is that the pulses remain in the same location ("-f0" and "f0"), because the cosine function period didn't change. Am I doing something wrong?

Thanks in advance ;)

PD: I know that increasing the sampling rate would increase the density of the spectrum and the time signal too, so N is going bigger automatically, because I would have more samples per second.

Generally I like to compose my sinusoids using this format:

x = cos(2*pi*f/fs*(0:num_samps-1))

Depending on the FFT routine you will use it will provide either the onesided or twosided result. Calculating the domain of the frequency axis is as follows. It should be noted that num_samps and NFFT may not be the same.

% Nfft - FFT size
% Fs - Sampling frequency in Hz
% oneside - flag indicating FFT results will be onesided or two sided

% Compute the frequency vector
if oneside == true
  f = fs * (0:Nfft/2).'./Nfft;
else
  f = fs * (-Nfft/2:Nfft/2-1).'./Nfft;
end

Have a look at this paper for a nice explanation on spectrum estimation.


Compute group delay of an audio file from STFT I'd like to reproduce results from this paper in which the authors used Group Delay-grams as an input to the Machine Learning model. The Machine Learning side is not the one that I have problems with. I have a problem with computing a Group Delay function of an audio file.

To begin with, I know that a group delay per se is a property of the filter. Authors of the paper, however, wrote something like this

  Group delay [10] is defined as the negative derivative of the
  phase spectrum of STFT:
  $$\tau(\omega,t)=-\frac{d(\theta(\omega,t))}{d\omega}\quad(2)$$
  As the implementation of Equation (2) requires the unwrap- ping  of
  the  phase  spectrum,  the  group  delay  function  can  be
  alternatively calculated using only the amplitude values:
  $$\tau(\omega,t)=\frac{X_{R}(\omega,t)Y_{R}(\omega,t)+Y_{I}(\omega,t)X_{I}(\omega,t)}{|X(\omega,t)|^2}\quad(3)$$
  where $R$ and $I$ denote the real and imaginary parts. $X(\omega,t)$ and $Y(\omega,t)$ denote the STFT of $x(n)$ and $nx(n)$, respectively.

This, I believe, is supposed to generate an image like in Fig. 3 in the paper, which looks like this (after either truncating or padding to the length of 256 along the time axis)

The problem is that when I try to compute this GD-gram from the equation above, the matrices multiplication cannot be done since both $X$ and $Y$ are non-symmetric. This is the original paper, that is referenced in both this paper and [10], where the authors are deriving the equation (3).

As I understood, the matrices multiplication is a standard one, not element-wise (element-wise multiplication doesn't produce the result that is in the paper). In my implementation below I'm using a window of length 800 samples since the authors of the paper use 50 ms window, which corresponds to 800 samples for 16 kHz sampling rate audio file. I tried to implement that in python as follows:

from scipy.io import wavfile
from scipy.signal import stft, get_window
import numpy as np
import matplotlib.pyplot as plt

rate, data = wavfile.read("test.wav")

n_data = np.multiply(data, np.arange(len(data)))
f, t, X = stft(data, rate, window="hamming", nperseg=800, return_onesided=False)
f_n, t_n, Y = stft(n_data, rate, window="hamming", nperseg=800, return_onesided=False)
group_delay = (np.dot(X.real, Y.real) + np.dot(Y.imag, X.imag)) / np.power(np.abs(X), 2)

but the code fails with

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-8-0058b9785caf> in <module>
----> 1 group_delay = (np.dot(X.real, Y.real) + np.dot(Y.imag, X.imag)) / np.power(np.abs(X), 2)

ValueError: shapes (800,109) and (800,109) not aligned: 109 (dim 1) != 800 (dim 0)

which is rather obvious. I tried to implement that in Matlab as well, hoping that maybe it's a python fault but it fails as well.

[x,fs] = audioread('test.wav');
n_x = x .* [1:size(x,1)];
X = spectrogram(x, 800, 400, 2048, fs);
Y = spectrogram(n_x, 800, 400, 2048, fs);
T = (real(X) * real(Y) + imag(Y) * imag(X)) / abs(X)^2;

How can I compute the Group Delay function using mentioned equation

The problem is that the stft function is splitting the signal up into different windows. That means that the signal from time $n$ to $n+N_{w}-1$ is multiplied by
$$ n, n+1, n+2, \ldots, n+N_w-1$$
instead of
$$
0, 1, 2, \ldots, N_w-1
$$
which is causing the scaling problem.

If I apply the group delay calculation from this derivation, I get:

where the top plot is the "data", the middle plot is the FFT of the data, and the bottom plot is the group delay. The group delay plot seems to be correct ($\tau_g = 49.5$, except for where the data is zero in the frequency domain.

To get this to work with the stft function, you're going to have to modulate the data with a sawtooth function rather than just np.arange(len(data)).

If I slice the sound file bachfugue.wav using the new, new code below I get the following image, which is more like what I'd expect.

New New Code Below

#!/opt/local/bin/python2.7
from scipy.signal import firwin
import numpy as np
import matplotlib.pyplot as plt
from scipy.io import wavfile

def group_delay(sig):
    b = np.fft.fft(sig)
    n_sig = np.multiply(sig, np.arange(len(sig)))
    br = np.fft.fft(n_sig)
    return np.divide(br, b + 0.01).real

rate, data = wavfile.read('bachfugue.wav')

N_w = 1024
L = np.floor(len(data)/N_w).astype(int)

full_delay = np.zeros((L, N_w))
for index in range(L):
    window = data[index:(index+N_w), 1].astype(float)
    full_delay[index][:] = group_delay(window)

plt.imshow(full_delay[0::10, 0::10], cmap=plt.get_cmap('hot'),aspect='auto')
plt.show()

New Code Below

#!/opt/local/bin/python2.7
from scipy.signal import firwin
import numpy as np
import matplotlib.pyplot as plt

N = 100
data = firwin(N, [0.05, 0.95], width=0.05, pass_zero=False)
n_data = np.multiply(data, np.arange(len(data)))
rate = 1

B = np.fft.fft(data)
Br = np.fft.fft(n_data)
group_delay2 = np.divide(Br, B + 0.01).real

plt.figure(1)
plt.subplot(311)
plt.plot(data)
plt.subplot(312)
plt.plot(np.abs(B))
plt.subplot(313)
plt.plot(group_delay2)

plt.show()

--

Original Answer

I've implemented it as I think it should be. I don't have your data, so I've just made a sinusoid.

The resulting plot is below. Code at end.

So I'm wondering what you should expect to see... more thought needed. I'll read that paper in more detail and see what I can make of it.

Now that I've had a chance to read through it, I'm wondering if that is correct.  Based on this derivation... I think the real way to do it might be different.  I'll try to code that up later.

Python Code

#!/opt/local/bin/python2.7
from scipy.signal import stft
import numpy as np
import matplotlib.pyplot as plt

N = 102400
rate = 1.0
data = []
for phi in ([2.0*3.14159265*x for x in range(0, N-1, 1)]):
    data.append(np.sin(phi))

n_data = np.multiply(data, np.arange(len(data)))
f, t, X = stft(data, rate, window="hamming", nperseg=800, return_onesided=False)
f_n, t_n, Y = stft(n_data, rate, window="hamming", nperseg=800, return_onesided=False)
# group_delay = (np.dot(X.real, Y.real) + np.dot(Y.imag, X.imag)) / np.power(np.abs(X), 2)

group_delay = np.divide(np.multiply(X.real, Y.real) + np.multiply(Y.imag, X.imag), np.power(np.abs(X), 2))

plt.imshow(group_delay)
plt.show()


Bandpass filter to get EEG frequency bands? I'm doing a project with an EEG to detect driver drowsiness and learned about EEG frequency bands (alpha, beta, gamma...).

As the frequency bands are simple frequency ranges, I wonder if I can use several bandpass filters to get them (instead of using WPT / FFT)?

Is this approach correct? Is there any reason not to do it (performance)?

As the frequency bands are simple frequency ranges, I wonder if I can use several bandpass filters to get them (instead of using WPT / FFT)?

Sure! That's how it's usually done!

  Is there any reason not to do it (performance)?

diabolical laugther as it happens, I've prepared just the blog post for you… TL;DR: If you don't have to process more than 20 Million EEG samples per second, your PC should do fine performance-wise.

Filtering Time-Series Data on the GNU Radio blog.

I came across an EEG problem that tried to do exactly that, filter things into multiple bands. And lo, I had 5 minutes of fun designing this bandpass filterbank for signals  with $f_\text{sample}=220\,\text{Hz}$:

$$
\begin{align}
\Delta:& [1,3]\text{ Hz}\\
\theta:& [4,7]\text{ Hz}\\
\alpha_1:& [8,9]\text{ Hz}\\
\alpha_2:& [10,12]\text{ Hz}\\
\beta_1:& [13,17]\text{ Hz}\\
\beta_2:& [18,30]\text{ Hz}\\
\gamma_1:& [31,40]\text{ Hz}\\
\gamma_2:& [41,50]\text{ Hz}
\end{align}
$$

in this GNU Radio Companion flow graph

The result was that, even in this totally naive filterbank implementation, things are so much faster than real-time that you can assume that if done a little more cleverly, it'd work on any halfway DSP-affine microcontroller in real-time.


Discrete Fourier transform - finding the fundamental quickly? First, I appologize because I'm a software developper and it's been a very long time I didn't dive into pure mathematics, so my question might seem dumb. I hope not.

The context is pitch recognition in music.

If you take a musical note, and apply a Fourier transform to it, you will have and infinite sum of amplitudes for given frequencies. For example, if I play a note whose fundamental is $F$, on any instrument, after Fourier transform, I will have harmonics at at $F, 2F, 3F,\ldots,nF$. Every frequency will have a given amplitude which defines the timbre of the instrument (piano, voice, trumpet, ... all follow this loaw, but you'll have different amplitudes for every harmonic)

Now what I'd like to do is from a given audio signal, find $F$. Just that.
It's more complicated than it seems because you will always have background noise and so on ... Further more, $F$ isn't necessarilly the frequency with the highest amplitude !

So my idea for finding $F$ is to apply a DFT (well actually an FFT for speed) and find a frenquency $F$, so that $F + 2F +3F + \ldots + nF$ is maximal in the FFT output.

Do you think that's possible at all ?
Do you think that's possible in a very short time (let's say < 5 milliseconds) ?

What you are describing is very similar the the Harmonic Product Spectrum method of pitch estimation, as listed in this Stanford CCRMA paper.

An FFT does not give you an "infinite sum of amplitudes", but a finite number of result bins depending on the length of the FFT.

5 mS is only 1 period of a 200 Hz note, and only a fraction of a period below 200 Hz.  Musical pitch recognition usually requires hearing or analyzing a multiple number of periods of the periodicity of a pitched sound.  And a lot of music uses notes below G2. If you have a sufficient length of data, computing a pitch estimate from that data might take only on the order of microseconds rather than milliseconds on a modern PC or mobile device.


Modified Periodogram function in MATLAB Given below is a function written in MATLAB for computing the modified periodogram, taken from Hayes, Statistical Digital Signal Processing and Modeling

function Px = mper(x, win, n1, n2)

x = x(:);
if nargin == 2
    n1 = 1; n2 = length(x);
end

N = n2 - n1 + 1;

switch win
    case 1
        w = rectwin(N);
    case 2
        w = hamming(N);
    case 3
        w = hanning(N);
    case 4
        w = bartlett(N);
    case 5
        w = blackman(N);
end;

xw = x(n1:n2).*w/norm(w);
Px = N*mypdg(xw);

Can anyone explain what is going on in second last line of the code
xw = x(n1:n2).*w/norm(w)

Is dividing by norm necessary to compute xw?

Yes it's explicitly necessary to normalize the periodogram for non-rectangular windows. However I cannnot see the mypdg script in the book ?

From Monson Hayes's book Statistical Digital Signal Processing and Modeling , chapter 8 (spectrum estimation), Section 8.2.3 (modified periodogram), equations 8.49

$$ \boxed{ \hat{P}_M(e^{j\omega}) = \frac{1}{NU} \left| \sum_{n=-\infty}^{\infty} x(n)~ w(n) ~ e^{-j\omega n} \right|^2 }$$

and 8.50

$$ \boxed{ U = \frac{1}{N} \sum_{n=0}^{N-1} |w(n)|^2 }$$

define the scaling to correct the biasing of the periodogram estimate due to using nonrectangular windows.

Note that in the program norm is used during windowing stage before the FFT, so that after FFT and squaring, you will get norm square which is identical to $U$ above.


Can the FFT reduce its computation load for specific sparse data patterns? I am working on a signal processing application, which involves an IFT.
My background is CS and not signal processing maths.
Consider an N-FFT configured with 256 "bins" (N = 256) for processing a
100Hz bandwidth. The application is in the a-priori position of knowing at
a given time where the input data can be sparse, in terms of only a subset
of the total bandwidth being active.
Ex: at time T, there is a signal in a band B of [26,40]Hz, while at [0,25)
and (40,100] Hz the signal is zero.
My questions are whether :

in such cases the 256-FFT could be configured to operate only
as a 16-FFT (as B = 40 - 26 = 15, 2^2 < 15 <= 2^4) , and potentially avoid
doing reams of computation as a 256-FFT where the add/multiply would be on
zero values (and thus redundant) .

if so, professional impls of the FFT already do this (a quick scan of
the input vector etc to determine whether such sparsity exists etc before
beginning the computation proper) .

Can the FFT reduce its computation load for specific sparse data patterns?

Sort of but only of it's very sparse. You can always NOT use an FFT but implement the inverse FFT directly and omit the 0 terms from the sum.
We can guestimate the break even point. An FFT takes about 10 operations per butterfly and for a length 256 FFT you have 8 stages with 128 butterflies each, so you end up with about 10000 operations per FFT
For a direct inverse DFT you'll have for each real output sample about 3 operations per non-zero frequency coefficient. In this case the break even point comes out to be 13 or thereabouts.
Your case of 16 is probably close to a wash, so I wouldn't bother.


How to perform dot annotation for crowd counting in a image file using Matlab? I have some crowd images having more than 300 people in each image. I want to label each person(i.e. head) with dot annotation such that coordinate of each dot(x,y) is saved in a .mat file. I hope there is something in Matlab for this task.

The Viola Jones algorithm does the face detection for you using the Cascade Object Detector.

%Read the image
x= imread('C:\~\crowd.jpeg');
x=rgb2gray(x);

%Face detection
faceDetector = vision.CascadeObjectDetector;
bboxes = step(faceDetector,x);

The first two columns of bboxes gives the (x,y) coordinates of the faces detected

If you want to label the people in the picture:

x1 = bboxes(:,1) %x co-ordinates
x2 = bboxes(:,2) %y co-ordinates
figure;
imshow(x);
for i=1:size(x1,1)
    plot(x1(i),x2(i),'*');
end

Saving the data in a .mat format:

save('dot.mat','x1','x2');

Accessing the data:

load('dot.mat');

To save data points from a given image, you can use ginput method. First read the image and then use [x,y] = ginput(n) where n is the number of points you need to save.


Python bandpass filter - singular matrix error I've been trying to design a bandpass filter using scipy but I keep getting a LinAlg Singular Matrix error. I read that a singular matrix is one that is not invertable, but I'm not sure how that error is coming up and what I can do to fix it

The code takes in an EEG signal (which, in the code below, I have just replaced with an int array for testing) and filters out frequencies < 8Hz and > 12Hz (alpha band)

Can anyone shed some light on where the singular matrix error is coming from? Or alternatively, if you know of a better way to filter a signal like this I'd love to test out other options too

from scipy import signal
from scipy.signal import filter_design as fd
import matplotlib.pylab as plt

#bandpass
Wp = [8, 12]   # Cutoff frequency
Ws = [7.5, 12.5]   # Stop frequency
Rp = 1             # passband maximum loss (gpass)
As = 100              # stoppand min attenuation (gstop)

b,a = fd.iirdesign(Wp,Ws,Rp,As,ftype='butter')
w,H = signal.freqz(b,a)  # filter response
plt.plot(w,H)

t = np.linspace(1,256,256)
x = np.arange(256)
plt.plot(t,x)

y = signal.filtfilt(b,a,x)
plt.plot(t,y)

I'd suggest checking again your parameter choices for

b,a = fd.iirdesign(Wp,Ws,Rp,As,ftype='butter')

Looking at the resulting filter coefficients, you have pole/zero cancellation at DC, i.e. at $z=1$.


How do I get the phase angle from Cross Correlation? Im trying to cross correlate two signals in matlab and get the phase difference between the signals.

For cross correlation (the idea is to do it without xcorr) I used:

Cxx=fftshift(ifft(fft(x,N).*conj(fft(y,N))))/(norm(x) * norm(y));

I get the result and not sure about the reference point I have to take for phase calculation Now which point should I take as zero on time scale and how do I get the phase difference from this result.

Thanks.

There is no such thing as the phase angle between two signals unless they both consist of a single sinusoid at the same frequency, that is, $x(t) = A\cos(\omega t+\psi)$
and $y(t) = B\cos(\omega t + \phi)$.  If you have $N$ samples of these signals $x(t)$ and
$y(t)$, taken at times $0$, $T$, $2T, \ldots$, $(N-1)T$, so that
$$x[n] = x(nT), ~~ y[n] = y(nT), 0 \leq n < N,$$ and either $N\omega T$ is an
integer multiple of $2\pi$ or $N\omega T \gg 1$
then the phase angle between the two sinusoids is
$$\theta
= \arccos\left(\frac{\displaystyle \sum_{n=0}^{N-1} x[n]y[n]}{\sqrt{\displaystyle
\sum_{n=0}^{N-1}(x[n])^2\displaystyle \sum_{n=0}^{N-1}(y[n])^2}}\right).
\tag{1}$$
For complex-valued sinusoids $x(t) = Ae^{j(\omega t+\psi)}$ and
$y(t) = BAe^{j(\omega t+\phi)}$, $(1)$ should be replaced by
$$\theta
= \arccos\left(\frac{\displaystyle \sum_{n=0}^{N-1} x[n](y[n])^*}{\sqrt{\displaystyle
\sum_{n=0}^{N-1}|x[n]|^2\displaystyle \sum_{n=0}^{N-1}|y[n]|^2}}\right).
\tag{2}$$

Of course, these formula can be used for arbitrary signals, not just for pure sinusoids,
but then, what you get is the _angle between the two vectors
$\mathbf x = (x[0], x[1], \ldots, x[N-1])$ and
$\mathbf y =(y[0], y[1], \ldots, y[N-1])$ in the $N$-dimensional spaces
$\mathbb R^N$ or $\mathbb C^N$, and not a phase angle between two sinusoids
at the same frequency.
Note that the three points $\mathbf 0 = (0,0,\ldots,0)$, $\mathbf x$ and
$\mathbf y$ lie in a (two-dimensional) plane in $N$-dimensional space and
the $\theta$ that you get is the angle between the line segments with endpoints
$\mathbf 0$ and $\mathbf x$ and endpoints $\mathbf 0$ and $\mathbf y$ which
lie in this plane. Another way to think about this is that
$$\langle\mathbf x, \mathbf y\rangle = ||\mathbf x||\cdot||\mathbf y||\cdot\cos(\theta)$$
and thus $(1)$ is obtained from
$$\cos(\theta)
= \frac{\langle\mathbf x, \mathbf y\rangle}{||\mathbf x||\cdot||\mathbf y||}
= \frac{\displaystyle \sum_{n=0}^{N-1} x[n](y[n])^*}{\sqrt{\displaystyle
\sum_{n=0}^{N-1}|x[n]|^2}\sqrt{\displaystyle \sum_{n=0}^{N-1}|y[n]|^2}}.$$

If you must use FFTs because that's the way you have been told to do it, then you have
$$\theta = \arccos\left(\frac{\displaystyle \sum_{n=0}^{N-1} X[n](Y[n])^*}{\sqrt{\displaystyle
\sum_{n=0}^{N-1}|X[n]|^2\displaystyle \sum_{n=0}^{N-1}|Y[n]|^2}}\right)\tag{3}$$
so that you have the ineffable pleasure of not only needing to
compute two FFTs first, but also of using complex multiplications in $(3)$
instead of the real multiplications in $(1)$ (for real-valued
signals).  This is overkill in my estimation,
but as usual, YMMV, and what your boss insists on is always right, regardless
of what people write on Internet forums.


How to determine whether a filter is high/low or band pass from the Z transform? How to solve questions of these kind?

I have tried by replacing $z=re^{jw}$ and taking the limits from $0$ to $\infty$. But I am not sure what $e^{j\infty}$ is.

HINT: There's no need to compute the frequency response by evaluating the transfer functions for $z=e^{j\omega}$. Just compute the poles and zeros and see where they are. Estimate the effect of the poles and zeros on the frequency response by trying to see how they influence the behavior of the transfer function on the unit circle.


What kind of power spectrum implies uncorrelation in time? I am new in signal processing field, and actually I just need some basic tools for what I am working on. I've learnt that the autocorrelation function is the inverse Fourier transform of the signal's power spectrum, and that if the power spectrum is flat then the autocorrelation function is the Dirac delta, so that the signal is uncorrelated in time.
But what if the power spectrum is constant only within the range of the cutoff frequencies?
In particular, I am given the following power spectrum for a signal $x(t)$:

$ S_{xx}(f)= \frac{1}{2(f_u-f_l)}$, for $\vert f\vert \in [f_l,f_u]$,

where $f_u$ and $f_l$ are the upper and the lower cutoff frequency respectively.

Does the uncorrelation still hold? If not, does the signal satisfy any property anyways?

Thank you in advance.

Since this does seem to answer the question, I'm reposting my comment as an answer:

I guess you may have to distinguish continuous (analog) and discrete (digital) signals. For analog signals, the autocorrelation function is a Dirac if and only if the spectrum is flat for all frequencies. Typically though what we care about is whether or not our digital samples are uncorrelated. It is easier to safisfy this. If your low-pass filter before your ADC is, say an ideal one, your ACF will be a sinc. If you get things right, your samples are in the zero crossings of the sinc and the discrete ACF looks like a discrete delta...

As the OP mentioned, the trick is to sample the signals in the zero crossings of the sinc function. For an ideal (rectangular) lowpass filter with (two-sided) bandwidth $B$ this means sampling exactly at a rate $1/B$.


Can I run a GA to optimize wavelet transform? I am running a wavelet transform (cmor) to estimate damping and frequencies that exists in a signal.cmor has 2 parameters that I can change them to get more accurate results. center frequency(Fc) and bandwidth frequency(Fb). If I construct a signal with few freqs and damping then I can measure the error of my estimation(fig 2). but in actual case I have a signal and I don't know its freqs and dampings so I can't measure the error.so a friend in here suggested me to reconstruct the signal and find error by measuring the difference between the original and reconstructed signal e(t)=|x(t)−x^(t)|.
so my question is:
Does anyone know a better function to find the error between reconstructed and original signal,rather than e(t)=|x(t)−x^(t)|.
can I use GA to search for Fb and Fc? or do you know a better search method?
Hope this picture shows what I mean, the actual case is last one. others are for explanations

Thanks in advance

You should be able to use GA or any other derivative-free optimization methods to solve your problem [1]. The correct parameter values can be determined using your application. If I'm trying to optimize for compression then I will look into the entropy of the resultant coefficients. Deciding on the actual values for your parameters will depend on how you plan to use them. The correct parameter values are simply the ones that work better.


how can I plot a power spectrum for 150 range bins with 512 fft points This is my code for plotting... However when I plot my spectrum data using this code, it is showing only 1 row. How can I plot all 150 rows? Can anyone help me to solve this?

 for i=1:length(dataheader.nrgb)
   plot(freq,arraydata(i,:)+height(i))
end

By default, matlab's plot command clears the plot area before drawing the new curve.

To prevent this, and allow multiple curves to appear on the same figure, you can activate the hold mode:

hold off; % 1st plot will clear figure
for i = 1:length(data header.nrgb)
   plot(freq,arraydata(i,:)+height(I));
   hold on; % set hold mode for subsequent figures
end

Alternatively you can stack multiple copies of freq and the plot command will draw one curve per column. Since your data is organized by rows, you'll need to transpose your data matrix:

% just make sure freq & height are column vectors. This may not be necessary if they already are.
if size(freq,1) < size(freq,2)
   freq = transpose(freq);
end
if size(height,1) < size(height,2)
   height = transpose (height);
end

% plot one curve per column on same figure
N = length(data header.nrgb);
plot(repmat(freq,1,N), transpose (arraydata)+repmat(height,1,N))


Can I use Fourier transforms instead of Laplace transforms (analyzing RC circuit)? I don't study electrical engineering or something related but I was assigned a problem on transfer functions, impulse responses, and in general, everything related to this post. (Specifically, I'm analyzing an RC circuit)

I do know the details of the Laplace transform but I prefer the Fourier transform. I wanted to know if there is any problem in substituting the Laplace transform for the Fourier transform in the analysis of LTI systems, impulse responses and transfer functions.

The important thing these transforms share is a Convolution theorem, which is highly relevant to the analysis of the things I've already mentioned above.

What do you think? What am I sacrificing by substituting transforms? Why is the Laplace transform favored over the Fourier transform?

Both transforms have a large overlap in their applications. So you can use both to analyze an RC circuit. However, with the unilateral Laplace transform it's much more straightforward to take initial conditions into account, such as an initially charged capacitor. This has to do with the unilateral Laplace transform of the derivative of a function:

$$\mathcal{L}\{x'(t)\}=s\mathcal{L}\{x(t)\}-x(0^-)$$

with the initial condition $x(0^-)$ occurring explicitly.

On the other hand, if you're interested in the behavior of the circuit when exposed to periodic input signals, then it's most natural to use the Fourier transform. Note that periodic signals (extending from $-\infty$ to $+\infty)$ don't have a Laplace transform.


How do I perform Spectral Analysis (FFT, Windowing, Detrending) on Sonic Anemometer Data? I am having trouble with all of the forums that I have looked at thus far.
I have one day's worth of Sonic Anemometer data. I want to see the spectral analysis of this data. When the anemometer data was being logged it appears that the frequency ranges between 31 and 32 hz. I am looking to analyze this data in thirty minute chunks and I would just like some guidance as to how to go about this.

I have already removed the bad data (stray data points with 98 m/s winds x.x).  I am not sure if I am supposed to detrend the data then apply a Hamming window and take the Fourier transform, or if I'm supposed to just apply the Hamming window and then take the Fourier transform.

Also should I calculate the sampling frequency for each 30 minute interval since it appears to change or is that change insignificant? (or should I analyze the whole data set together?)

Sorry for the long post but I have been hacking away on forums all week and I cannot seem to find a consistent solution.

take fft of data

let's say data is x = [ 3 4 5 6 6 9  7 ]

take fft

y = fft(x);
shift so -fs/2  to fs/2
y = fftshift(x);
plot mag
plot(abs(y))

see link for help
fft Help


Interpreting Coefficients from wavelet The wavelet transform informs  what all derived wavelets (from mother wavelets) with different scales at different point of time , are present  (along with the associated strength) in my original signal.  I am not able to understand how do I now correlate it with FFT that simply gives what all sine and cosine waves are present.

Now how do I make use of this information that I get from WT? I cannot say unlike as in case of FFT or STFT " we have a 100hz frequency in the window of 4 seconds". I will have to say "we have a db4 wavelet present  in some window". How it will make sense. How do I interpret this in a meaningful way?

Wavelets are used for time-scale, or time-frequency, analysis, just like STFT. You can map the coefficients of CWT to frequencies through the center frequency of the wavelet. See for example http://www.mathworks.se/help/wavelet/ref/scal2frq.html


Ground plane estimation from a single image (was: confusion about homography and calibration camera parameters) I'm new to this multiple view geometry issues and I'm a little bit lost, and some concept are not clear at all.

I have a calibrated camera ( = I have distortion coefficents, field of view, and with those I can find what is called intrinsic matrix).

With this camera I'd like to make an application that can estimate the ground plane from a picture. I have those constrains:

In the picture there is no chessboard (or other kind of known pattern used in computer vision for calibrating cameras)
I can manually define some coplanar points (like points that I know are the projection of a square) or some

So I need something that take the 2D coordinate of a trapezoid in the picture, knowing that in the real world it is a unit-square on the ground plane in the 3D real world coordinates and give me the transformation.

From what I have understand this is procedure is "computation of extrinsic parameter given given intrinsic parameters, a few 3D points, and their projections".

Now, my questions:

Is it possible?
Are 4 points enough?
This is not related with homography, right? (because homography is not related with 3d points but only with 2d).. but how can homography help me? (I can find the homography that transform the trapezoid in a square)
Is there an alternative way to do it (like with triangles and normals)?

For the sake of curiosity I'm using openCV for the implementation.

Yes it is possible to compute the extrinsics given the intrisics, some points in 3D and their projections in the image.  If all your 3D points are in the same plane, then the math for computing the extrinsics is explained in the paper by Zhengyou Zhang, which is the basis for the camera calibration code in OpenCV.

If your 3D points are not co-planar, then you need to use a different algorithm that is described in Hartley and Zisserman.


What is energy in image processing? There are several web articles on the topic, but, after having read them, I still don't know what "energy" is, in the context of computer vision (CV) and image processing (IP). It is related to optimisation, given that, apparently, there are the so-called "energy minimisation problems".

Why is it called "energy" (that is, why the word "energy")?  I read it is because it is related to the concept of "energy" in physics, but assume that I have no idea of what energy in physics is. How would you explain what energy is in IP and CV? What is the intuition behind it?

Why do we care about energy? Given that there are the energy minimisation problems, we want to minimise energy, in such cases. I have also read that there are energy maximisation problems, so sometimes we also want to maximise energy.

I have read that it has different meanings depending on the context, but I just want to know its meaning in image processing and computer vision, that is, in the context of images.

The usage of the term "energy" in image processing context has historical reasons. It can be mapped to "brightness" or "intensity".

Imagine you have a (for matters of simplicity) greyscale image, that you want to transmit. Historically, that was done via some kind of analog modulation, where, simply speaking, dark "pixels", low intensity meant "low voltage" and bright "pixels", high intensity meant "high voltage". When looking at a CRT television, you could take the word "energy" quite literally, as an all black image will produce significantly lower energy consumption than an all white image.

The term "engergy minimisation" is usually employed, when trying to optimize some kind of algorithm or codec by comparing the difference of original signal and encoded-decoded signal. The codec would be perfect, if this difference was all zeros, so the target is "to minimise the energy" of the difference signal. Usage of "energy maximisation" is analog.

So, think of "energy" as "intensity" or "brightness" in the context of image processing.


Does the system matrix being singular tell us anything about the system? If we have a linear system, represented in State Space and the A matrix is singular (det(A) == 0), can we expect any special properties from the system?

If the system matrix $A$ is singular then it has at least one eigenvalue equal to $0$. Since the eigenvalues of $A$ represent the poles of the transfer function of the system, there is at least one pole at $s=0$, which makes the system unstable (or marginally stable if there is a single pole at $s=0$, but in any case not asymptotically stable).


Interpolation by factor of 2 I want to up-sample my input signal by a factor of $2$. I saw zero padding followed by low pass filter method being used in few cases. But still I need some help in this.

Say I have $10$ input samples and my interpolating factor is $2$, then if I am doing interpolation using above method, will I get $20$ samples at the output?
Also how to increase the pass-band gain of the low pass filter to $6.02\textrm{ dB}$?

When you upsample by a factor of two by inserting a zero between each sample you create an alias of your signal centered at the new Nyquist frequency (half the new sample rate).  You get rid of the alias by low-pass filtering.

The specifications of the filter that you need depend on your signal.  Say that your signal has one-sided bandwidth $B$ (by one-sided I mean that we are only talking about the positive frequencies, not the negative ones).  $B$ must have been less than the Nyquist frequency of the old sample rate, $f_{sOld}$, to avoid aliasing.

Once you upsample the alias appears.

Hopefully these pictures will help to gain an intuitive grasp of how to characterize the required low-pass filter.  The passband must be flat from 0 Hz to $B$ Hz, and it must cutoff at $\frac{f_{sNew}}{2} - B = f_{sOld} - B$.  The amount of attenuation that is required in the cutoff region depends on your application.  Usually 50 dB is more than sufficient.


Noise variance calculation for 16APSK Demodulator I need to use approximate LLR algorithm for 16APSK demodulation. For that, I need to calculate noise variance which is given as

$Noise\;variance = \frac {Signal\;power\;\times\;Symbol\;time}{Sample\;time\;\times\;10^{\frac{Es/No}{10}}}$

I have the following informations
SNR = 10, Signal Power = unity.
Number of bits/symbol = 4,
Code rate = 8/9

How can I know symbol time and sample time values? If I have to assume, what are the typical assumptions?
And if Es/No = SNR? or it is dependent on  symbol time and sample time?
So far I have assumed Symbol time = Sample time.
Though I can understand symbol time, I am unable to comprehend the meaning of sample time in  this context.

This is to take into account the oversampling operation.

Symbol time = sample time implies no oversampling.

See AWGN model for more details about the conversion among EbN0, EsN0 and SNR.

  For example, if a complex baseband signal is oversampled by a factor of 4, then EsNo exceeds the corresponding SNR by $10\log_{10}(4)$.


Frequency response of a (rectangular) integrator A very simple question which can't really figure out:

I have a simple discrete time rectangular integrator with Z transform H(z) = 1/(1-Z^(-1))

Plotting the frequency response with Matlab using freqz shows what it looks like a kind of low pass filter response.

However, when finding the freq response by sending an delta pulse through the integrator via filter Matlab function, and performing the FFT on its output, I just get a pulse (with amplitude equal to number of FFT points) followed of zeroes.

I guess this kind of makes sense mathematically, as the output of the integrator is a step function, which, after FFT, is the addition of as many samples of the step output as the FFT length.

However,

Why does the freqz method does not match the impulse response + FFT approach?
Am I wrong in trying to look at an integrator as a digital filter?

Thanks a mill

PD:
I forgot to mention that this is the closest to this question I could find on the site... However it is a leaky integrator not a rectangular (i.e. feedback multiplier is 1) integrator

Is a leaky integrator the same thing as a low pass filter?

Those are good questions! Here's my "two cents." Hilmar is correct, the integrator is not stable because there's a pole on the z-plane's unit circle. But I'll bet a pint of beer that the region of convergence is outside the pole (and therefore outside the unit circle) because the integrator's impulse response is causal.

Answer to question# 1: Matlab's freqz() command computes the freq response at zero Hz by evaluating a ratio of polynomials at various frequencies. And that ratio has a zero-valued denominator at zero Hz. So freqz() gives a freq mag response of infinity at zero Hz.  The DFT of an impulse response that's a finite-length sequence of ones, evaluated at zero Hz, is merely the sum of the impulse response samples.  The two approaches result in freq magnitude responses that are zero-valued at all frequencies except at zero Hz.  Both approaches give similar results.

Answer to question# 2: You could, if you wish, call an integrator a "filter." But it's unstable, and unstable filters are of no value in trying to filter real-world signals.


MU-MIMO Sensitivity to Line of Sight (LOS) Why do researchers talk about MU-MIMO being "less sensitive to line of sight" / "LOS Okay"? For example in this video, time=[00:40,00:42]

Similarly, "User separation based on spatial channel properties is
particularly  difficult  in  situations  where  the  users  are  located
close to each other and experience LOS propagation conditions
to the BS antenna array"

My perception is that LOS is always preferred to the opposite, multipath propagation.

Can someone please explain this? Thank you.

There is no mention of light-of-sight (LOS) on the time you mentioned. There is a mention of point-to-point, however. Point-to-point channels could be LOS or multipath.

MU-MIMO is better than point-to-point MIMO because users are usually well-separated in space. MU-MIMO is particularly better when the channels are correlated at the base station (BS). The reason for that is that the spatial separation of the users make the channel correlation matrices from different users to be different. In the ideal case where these correlation matrices eigenspaces are orthogonal, the data from different users can be perfectly separated without any interference. In point-to-point MIMO, on the other hand, the data from different antennas from the multiple-antenna transmitter will have the same correlation matrix at the BS, and thus there will be more interference from other data streams.

Read section 2.4 from the book Massive MIMO Networks: Spectral, Energy, and Hardware Efficiency.

EDIT: In the paper, they are talking about separating users that have LOS channels with the BS in massive MIMO, that are difficult to separate in traditional MIMO.


Trying to compare two speech samples that have been recorded by smartphones in pockets I'm trying to compare two audio samples of speech from two different Android smartphones to determine if the people are co-located. I'm doing this using the frequency bins of an FFT between 300 Hz - 3400 Hz. When the two devices are sitting on a table, I have good results.

However, when the devices are placed within pockets, the audio is recorded differently than when they are placed on a table.

Does anyone know exactly what happens to the audio frequencies/signal when the microphones are covered by a pocket or a bag? Are certain frequencies filtered?

The distance emitter-receptor and the presence of obstacles in the trajectory(clothes for example) will afect the recorded signal for sure. So some frequencies can be attenuated or just disappear.


Convolution vs multiplication vs polynomial multiplication? Apparently these 3 terms are some how similar but what is exact difference between them?please kindly explain with example, because the number of terms in result/output are different in all 3 cases

Especially in  context of signal processing, when do we need to multiplication, when we need to do convolution and when we need to do polynomial multiplication?

The convolution theorem is a fundamental signal processing property that links the frequency and time domains. Specifically convolution in one domain is multiplication in the other. Often this can be used to simplify the number of operations required for solving signal processing problems. The utility of this is amplified given the significance of the Fast Fourier Transform in how efficiently it can be used to translate between domains.

With this in mind, consider that the output of a linear system in time is related to its input using convolution:

$$y(t) = x(t)\star h(t)$$

Where $\star$ represents the convolution operator, $x(t)$ is the input and $h(t)$ is the impulse response of the system.

If we took the Fourier Transform of $x(t)$ and $h(t)$, then this same problem is solved using an element by element product:

$$Y(\omega) = X(\omega)H(\omega)$$

Similarly we use convolution to multiply polynomials by convolving the coefficients. So we could use the FFT for polynomial multiplication as follows:

Consider the simple example of a product of $x^3 + 2x^2 - 4x +1$ with $5x^2+3$

This is solved by convolving the coefficients as:

$$[1, 2, -4, 1] \star [5, 0, 3] = [5, 10, -17, 11, -12, 3]$$

Resulting in $5x^5 + 10x^4 -17x^3+11x^2-12x+3$

Similarly this could be done in the frequency domain with FFT processing, but to do this properly the coefficients of the polynomial are loaded starting with the lowest power of x and we must zero-pad out to 6 samples that we would expect in the result:

FFT of $[1,-4, 2,1,0,0] = [0, -3+1.73j, 3+5.20j, 6, 3-5.2j, -3-1.73j]$

FFT of $[3,0,5,0,0,0] = [8, 0.5-4.33j, 0.5+4.33j, 8, 0.5-4.33j, 0.5+4.33j]$

Element by element product of the FFTs: $[0, 6+13.85j, -21+15.59j, 48, -21-15.59j, 6-13.95j]$

Inverse FFT of the product:  $[3, -12, 11, -17, 10, 5]$

Corresponding to $5x^5+10x^4-17x^3+11x^2-12x+3$

For something this small as in this example, it would be more efficient to do the convolution directly, but as the length of the polynomials grow, the FFT approach becomes significantly more efficient.


Why do we use scaling factors in DCT? Are they important? I don't understand the use of scale factors, why do we need them? Can we not use them or are they essential?

The wikipedia page on DCT II shows this:
https://en.wikipedia.org/wiki/Discrete_cosine_transform

Scale factors ensure the normality of each DCT vector. Since the vectors are already orthogonal, factors are needed whenever the orthonormality of the DCT is needed (often compression, denoising). For visualization, comparison, feature extraction, scaling factors may be less important.

For operation efficiency, especially when integer arithmetic is used, the scaling factor can be incorporated afterward, in the quantization stage when used in compression.


Forecasting with ARMA models, from a filter point of view ARMA
 models are afaik just filters with transfer function
$ {MA(z) \over AR(z)}  \equiv {FIR(z) \over IIR(z)} $ .
However forecasters of stock prices, market trends ...
seem to be mainly statisticians, with their own vocabulary and culture.
For example, "signal-to-noise ratio" is rarely mentioned;
for another, differencing must increase noise.
Can anyone suggest either

textbooks or introductory courses on ARMA forecasting
from a filter or signal processing point of view
websites with real time series and running code to ARMA-model them ?

(I'm interested in ARMA models for prediction, not in spectral analysis as such.
ARMA models may well be wrong for prediction from short, noisy data --
what will the economy do next year ? --
hence the need for real examples.)

Added: Some 40 years ago, R.W. Hamming wrote
in Digital Filters:

  ... We have a predicting filter without finding, or even talking about, the transfer function.
  Statisticians often do this,
  and neglect to examine the corresponding transfer function of the formula,
  which can often shed some light on the whole system.

ARMA divides the signal into two parts and that models the two parts.

Financial time series are corrupted by different types of correlated and uncorrelated noises with definite functions that allow modeling and others more difficult as per having to use aproximations. In addition financial or economic times series may exhibit long memory processes (Black noise) or may have other long range dependencies or other non stationary issues. It is for this reason that besides the MA part it is hard to tell if a series is stochastic or deterministically chaotic or a mix between these.

The two terms of ARMA include an autoregresive term and a moving average term. The deterministic part of time series is modeled with ARIMA, ARFIMA etc family of models and the stockastic part using ARCH and GARCH family of models (usually volatility is modeled this way as a large portion of the stockastic part is white noise which can not be forecasted but its variance can).

Diferenciation as you mention may eliminate some of the signal, however is necesary for regression in many cases as to make a time series weak sense stationary which is a requirement to have a meaningfull and valid regression.

Signal to noise ratio is hard to guess for any unknonw noise+signal composition. However there are some ways to aproximate this by determining what noises affect your signal and what weights of noise correspond to each signal as well.

I would sugest from a signal procesing point to google / take a look at "financial / non stationary time series signal procesing" papers. There are many methods and ways as to extract information from a signal.

I would sugest from a modeling point of view that you look at ARIMA or GARCH presentations from academia or papers from leading economic schools (LSE etc).

I will also sugest taking a look / google data preprocesing for regression (including the terms of unbalanced data set, stationarity, normality, linearity, aditivity and serial correlation, homocedasticity, spurious regression, cointegration, transformation of variables for financial time series).


Maintain constant total amplitude of sum of sines Given a set of N sine waves, each with their own amplitude and frequency, is there a straightforward way to select their phases such that the sum over all sines maintains a constant amplitude, given that the frequencies are spaced equally apart by df?

If all phases are 0, then at t=0 all the sines will constructively interfere to make a huge spike, then they will destructively interfere until they become commensurate again at t~1/df. My goal is to select a set of phases for the different sine waves which will leave the resultant sum with roughly constant peak-to-peak amplitude from t=0 to t=1/df. Picking all random phases does do a pretty good job, but ideally I'd like to find a way to find an optimal set of phases, or at the very least come up with a metric to quantify "how good" a given set of phase is at accomplishing my goal, which I can then optimize on.

edit: To clarify I'm not looking to have the sum of sines equal a constant value, but rather have the resulting periodic function have a roughly constant peak-to-peak amplitude. For instance, for a given set of frequencies all with equal amplitudes, if all phases are set to 0 then the resulting sum will look like:

If the phases are set uniformly randomly in the range [0 2pi] then the result is:
.

Clearly, in the random-phase case the peak-to-peak amplitude of the resulting oscillatory function is roughly constant over the given time series, while in the 0 phase case it spikes often at times where commensurate sines add constructively. I'm looking for a way to find an optimal set of phases which keeps the peak-to-peak amplitude roughly constant, without just generating a million different random combinations and picking the best one.

Any help would be greatly appreciated.

Hopefully this will be some help.

As you have stated it, this is a difficult problem to formulate and thus get an analytical solution.  First, are you talking continuous or discrete?  Makes a big difference.

Anyway, you want to formulate it as an optimization problem.  Something like "minimizing the RMS" can be expressed as an evaluation function.  "All the peaks having the same amplitude" means have to measure the value of the function with the condition applied to the derivative (=0).  So, your evaluation function (a function of your free variables) could be something like:

$$ F(\phi_0,\phi_1,\phi_2,...,\phi_{N-1}) = \sum_{p=1}^{P} \left( f(t_p) - \bar f \right)^2 $$

Where $f(t)$ is your function, evaluated when $f'(t)=0$, and $\bar f$ is the average height of all those peaks.

This evaluation function is going to be periodic in all dimensions.  There are going to be lots of local minima.  Your only recourse (as I can see it) is to do a coarse sweep of your whole domain, then a local optimization search on the best candidates.


What Distance Measure Is Proper for Image Patches? I have a bunch of image patches extracted from two images.

I want to find a proper measure of similarity between each patch from first images to patches from second image.

I have made vectors from each patch to calculate distance/similarity.
I have tried to use euclidean distance and closest vector by phase, but these two measures do not perform well when reconstructing first image by replacing each patch by the most similar patch of the second image.

Is there any well suited measure of similarity for small patches of two images?

Well, You're asking on one of the most researched topics.

There is no 1 "Right" answer.
Many models have been suggested and many of them works pretty well.

In the Non Local Means they use Euclidean Distance.
Later revisions on the idea used Mahalanobis Distance.

Others are using different norm, other use Learned Dictionaries.

Some use features or transformation of the data.

I think the best method depends on the image (Or set of images) and the usage (Denoising or something else).


Modified Tukey window applied in frequency domain In analysis of transient signals (ground motion time-series) in seismology, many operations are performed in frequency domain and it is common practice to apply a filter which is like a modified Tukey window in frequency domain to keep data band-limited before applying an inverse IFFT to go back to time-domain. Essentially, 4 frequencies f1, f2, f3 and f4 are specified and the filter is a combination of high pass and low pass taper. Here quoting from a web manual-
"The taper is unity between f2 and f3 and zero below f1 and above f4. Frequencies f1 and f2 specify the high-pass filter at low frequencies, while frequencies f3 and f4 specify the low-pass filter at high frequencies. Both f3 and f4 should be less than the Nyquist frequency: 0.5/DELTA. The filters applied between f1 and f2 and between f3 and f4 are quarter cycles of a cosine wave. To avoid ringing in the output time series, a suggested rule-of-thumb is f1 = f2/2 and f4 >= 2*f3".

Can someone please explain what is the reason behind the rule of thumb mentioned at the end ?
Is there any specific advantage of using such a filter compared to other filters like an acausal Butterworth filter ?

I have added links to 2 figures (for some reason, upload image is not working for me) which show time and frequency domain response of 4 filters (cosine tapers with f2=2f1 and f4=2f3 [black], and f2=1.25f1 and f4=1.25f3 [red] with f2=0.2 Hz and f3=0.9 Hz) and 2-pass bandpass Butterworth filters with 4 and 8 poles applied at frequencies f2 and f3.

Regarding the 2:1 frequency ratio recommendation for the transition band (f2:f1 or f4:f3):

The steeper the filter response (i.e. the smaller the frequency ratio between pass band and stop band limits), will necessarily produce a much longer impulse response and produce ringing artifacts in the time domain whenever there is a transient.

These filters, which are applied in the frequency domain, also seem to have linear phase response, which is very desirable for analysing this type of signal.

Butterworth filters have a compact and nice representation as difference equations (for filtering in the time domain), but do not have linear phase response, and as you see from your plots, have not very good transition between pass band and stop band.

(I don't understand how you plotted a symmetric impulse response (time domain) for Butterworth filters; their impulse response is causal and asymmetric in general.)


edge detection evaluation methods I'm working on some edge detection algorithms that I want to evaluate them. I want to know what are standard ways to evaluate them that check does them work better that old ones or no.

thank you for helping

I strongly agree with Ava. Without a clear formulation of your task, you cannot make a meaningful ranking of algorithms.

Think about your most important task and formulate clear goals. Here is an example of such an investigation:

"Detect only the outer edges of a shoe box in grey value images in a short amount of time".

Goal 1: Detect only the outer edges --> measure number of connected components, if 1, check if it is the contour of the parcel, if !=1, vary parameters and try again.

Goal 2: Be fast --> If goal 1 is reached, measure time needed to compute these edges.


variability among signals I have a set of $n$ low frequency time-varying signals $y(t)$.

I would like to express the variability within this set of signals using a single number. My current idea is to:

Compute the average signal $\bar y(t)$ for the set using

$$\bar y(t) = \dfrac{y_1(t) + y_2(t) + \ldots + y_n(t)}{n}  $$
Compute the 'standard deviation' of each individual signal with respect to the mean signal $\bar y$ using the Pearson correlation $r$

$$\bar r = \dfrac{ r(y_1(t), \bar y(t)) + r(y_2(t), \bar y(t)) + \ldots + r(y_n(t), \bar y(t))}{n}$$

This would then give me a sense of the variability within the set of signals.

My question is whether there are better, different, or more established ways of doing this?

As an example of what I mean, let's have a look at a simple sinusoid (though I see you're actually assuming double gamma functions; I may rework this later to use those).

This models the signals as:
$$
y_p(t) = \sin(\omega t + \phi_p)
$$
so the only difference between them is the phase $\phi_p = (p-1)\pi/6$.

The top plot is the mean signal $\bar{y}(t)$ and the bottom plot shows the Pearson coefficient between $\bar{y}$ and $y_p$ for $p=1\ldots 6$ (so $n=6$).

As you can see, the Pearson coefficient varies quite widely even in this simple example.

R Code Below

# 26438

Nsigs <- 6
T<- 1024
omega <- 2*pi*0.00982734982
t <- 0:(T-1)

y <- matrix(,nrow = T, ncol = Nsigs)
phi <- c(0, 1*pi/6, 2*pi/6, 3*pi/6, 4*pi/6, 5*pi/6 )

for (idx in 1:Nsigs)
{
  #phi[idx] <- runif(1)*2*pi
  y[,idx] <- sin(omega*t + phi[idx])
}

mean <- rowSums(y) / Nsigs

total <- matrix(,nrow = T, ncol = Nsigs+1)
total[,1] <- mean
total[,2:(Nsigs+1)] = y

pearson <- cor(total,use="complete.obs", method="pearson")

par(mfrow=c(2,1),pty="m")
plot(t,mean, type="l")
title("Mean value")

plot(1:6, pearson[1,2:7], pch=19, col="blue")
title("Pearson coefficient")


DSP Book with More Applications The process of my DSP class used the Oppenheimer DSP book as text. Although I can easily understand the math in the book, I found it lack examples and applications, hence I still cannot build a picture how people actually use those under different circumstances. Can someone recommend me a DSP book that contains more practical examples and applications more than "This example uses the mathematical formula above to derive some numbers from some numbers"? Thanks in advance.

Hopefully you'll get a bunch of answers here from the very general to the super specific. I'll put in my two cents here.
The recommendations I would make are from the field of radar and communication systems. These systems tend to exercise almost all aspects of signal processing:

Signal generation and mixing
Sampling, decimating/upsampling
Signal conditioning
General filter design and implementation
Detection theory
Target/observable tracking
Advanced techniques such as eigen-decomposition and super-resolution

Unfortunately books in the radar and communications fields are not free, but in my opinion are worth their asking price (as far as texts are concerned).
Bassem Mahafza: Radar System Analysis and Design using MATLAB, Chapman & Hall
This book is an awesome treatment of many of the aspects of a radar/comms system. As the title suggests, buying this book gives you access to all of the MATLAB functions provided by the author. As you work through the theory of the book, many sections have accompanying examples that you can play with allowing you to both see implementation details and visuals.
Merril I. Skolnik: Radar Handbook, Third Edition, McGraw-Hill
This book is a very large collection of different authors' work. It covers almost everything you can imagine about signal processing in a radar system. Theory, practicalities, and implementation specifics are all covered here. This book is to be used as a reference for applications and does not really contain examples of how to code certain things.
Harry L. Van Trees: Optimum Array Processing: Part IV of Detection, Estimation, and Modulation Theory, Wiley & Sons
This one is dense and I would consider to be the more hardcore of the three. This text exposes to you a bunch of different techniques in signal processing when you have access to an antenna array, which would not be possible otherwise. Applications include radar, astronomy, communications, seismology, medical diagnosis, etc.
I could list many more but these are the three that come to mind that are a good mix of theory, applications of that theory, and examples.


What is difference between along the edge and across the edge I am trying to implement a paper called structured tensor based image interpolation. In the paper I came across a concept called tangent direction and normal direction. What they mean is tangent direction is along the edge and normal direction is across the edge. Can anyone help me to find out what the tangent and normal direction the author is talking.
The paper link is below:-
https://arxiv.org/pdf/1402.5564.pdf

This is illustrated for a step edge. If you look at gray-scale values, one often considers that a 1D edge is (relatively) irregular across the edge (from gray to black), and regular along the edge (parallel to the border, the luminance value does not vary).

You can perceive it better in the 3D view, where the amplitude along with the gray value is represented in 3D.

In more generic cases, at the appropriate scales, the structure tensor is supposed to capture this behavior.


Impulse Response Function-Green's Function Assuming a input-output system
$$u(t) = c \cdot \big(h(t)\circledast m(t) \big)$$

where its impulse response function is:

$$
h(t) =
     \begin{cases}
       \frac{A}{\tau}te^{-\frac{t}{\tau}}, &\quad t\geq0 \\
       0, &\quad t<0 \\
     \end{cases}
$$

Given only these, how can one derive the differential equation below, whose Green's function is the above?:

$$
\ddot{u}(t)+\frac{2}{\tau}\dot{u}(t)+\frac{1}{\tau^2}u(t) = \frac{A}{\tau}m(t)
$$

Would, considering an inverse n-order L differential operator and performing the consequent integrals, be a possible way?

A simple way to derive the differential equation from the impulse response is to transform the latter to the frequency domain, rewrite the input/output relation, and then transform the resulting equation back to the time domain.

The Laplace transform of the given impulse response $h(t)$ is

$$\begin{align}
H(s) &=\frac{A}{\tau}\frac{1}{\left(s+\frac{1}{\tau}\right)^2} \\
&=\frac{A}{\tau}\frac{1}{s^2+\frac{2}{\tau}s+\frac{1}{\tau^2}} \\
&=\frac{U(s)}{M(s)}\tag{1}
\end{align}$$

where $U(s)$ and $M(s)$ are the Laplace transforms of the output $u(t)$ and input $m(t)$, respectively.

From $(1)$ we get

$$U(s)\left[s^2+\frac{2}{\tau}s+\frac{1}{\tau^2}\right]=\frac{A}{\tau}M(s)\tag{2}$$

Transforming $(2)$ back to the time domain finally gives the desired differential equation

$$\ddot{u}(t)+\frac{2}{\tau}\dot{u}(t)+\frac{1}{\tau^2}u(t)=\frac{A}{\tau}m(t)\tag{3}$$


Typical filter orders I am on my first DSP project where I want to implement lowpass/highpass filters for ECG application. I have started from MATLAB to give me orders and coefficients for my desired specs. I prefer to go with FIR filters but as you know they require much higher orders than IIR filters. Since this is my first project I have no clue what ballpark should I have in mind for a filter order that runs in a DSP.(Analog Devices BlackFin)

So I decided to frame this question as a general call for advice, that what filter orders are DSP engineers comfortable with choosing for implementation. Is 100's of order very high for a DSP? Do they go with less than order of 10? 50?

I understand this is very much dependent on filter specs and my application. But I still want to know if there is a range to have in mind.

I hope I have been able to convey my question.

my experience in audio is that for IIR filters, there is usually a cascaded Second-order Section (SOS) for every little feature (a bump or an edge in the frequency response deliberately placed there by the user).  often just a single second-order IIR is used for many little jobs where some frequency is given a little boost or cut.  we call those IIR filters a Parametric Equalizer.

however, for a fixed and sharp feature (like a brickwall of some sort), then the order of the IIR filter will be related to the sharpeness of that feature.  for a brickwall, there is a 6 dB per octave of steepness in the wall for every integer order of the filter.

for FIR, Fat's answer is good.  the order of the filter is 1 less than the number of taps of the FIR.  the minimum order again is related to the sharpness of some particular feature in the frequency response.  there are a couple of heuristic formula for initial guessing at what order you will need.  i will try to express them in a way to compare apples-to-apples.

using a Kaiser-windowed FIR design the filter order (one less than the number of taps) is about:

$$ N \ge \frac{(A-7.95)F_\text{s}}{14.357 \cdot \Delta f} $$

where $A$ is the amount of attenuation in dB of the brickwall, $\Delta f$ is the frequency width of the transition between the nominally un-attenuated and attenuated portions of the frequency response, and $F_\text{s}$ is the sample rate.

using Parks-McClellan will reduce $N$ a bit:

$$ N \ge \frac{(\tfrac12 A + \frac{4.343}{r_\text{p}} - 29.7)F_\text{s}}{14.602 \cdot \Delta f} $$

where the new symbol $r_\text{p}$ is the passband ripple in dB.

this appears to be twice as good as i was expecting.  so i am not sure i did this right, but i was being very careful with the equations from O&S about the P-McC alg.


Normalizing audio waveforms code implementation (Peak, RMS) I have some audio data (array of floats) which I use to plot a simple
waveform.
When plotted, the waveform doesn't max out at the edges.
No problem - the data just needs to be normalized. I iterate once to find the max, and then iterate again dividing each by the max. Plot again and everything looks great!
But wait videos which have a loud intro, or loud explosion, causes the rest of the waveform to still be tiny.
After some research, I come across RMS that is supposed to address this. I iterate through the samples and calculate the RMS, and again divide each sample by the RMS value. This results in considerable "clipping":

What is the best method to solve this?
Intuitively, it seems I might need to calculate a local max or average based on a moving window (rather than the entire set) but I'm not entirely sure. Help?
Note: The waveform is purely for visual purposes (the audio will not be played back to the user).

What is the best method to solve this?
  Intuitively, it seems I might need to calculate a local max or average based on a moving window (rather than the entire set) but I'm not entirely sure. Help?

What you are describing with your second of these two quoted bullet points is Automatic Gain Control, which eventually will take us to compression. What you might find useful in your case, i.e for visual purposes only, is a simple logarithmic transformation which is very close to what compression is about but without the complex controls it demands.

The most common way to solve such problems in the audio world is to use a compressor.

The compressor is monitoring the amplitude of the incoming signal. If that amplitude crosses a user-defined threshold it decreases the amplification in an attempt to keep the signal within certain limits. As you can understand, imposing such limitations to a signal without special care could result in annoying audible artifacts. For this purpose, the compressor has an additional set of controls called "Attack" and "Release". These are two parameters that control how fast does the compressor impose its limitation once the threshold has been crossed (Attack) and how fast does the compressor return to the initial amplification of 1:1 after the amplitude of the signal has returned below the predefined threshold. These three parameters, the threshold and attack and release times are what make the compressor "complex" in setting it up properly. Extreme values in these settings can result in more artifacts.

You may have "heard" a compressor working in radio shows. It is the effect that you hear when a speaker suddenly stops talking and after a very brief period of time the background sounds seem to rise in volume. When the speaker starts talking again, the background sounds seem to be diminishing abruptly again and the speaker's voice now dominates the dynamic range of the signal. That's the effect of a compressor.

The main idea here is that the compressor is making soft sounds louder and loud sounds softer.

  Note: The waveform is purely for visual purposes (the audio will not be played back to the user).

If you are not looking to actually apply dynamic range compression and you simply want to make sure that your viewer can see all details within a signal, you can apply a logarithmic transformation, which is very similar to how a compressor / limiter would operate but without the complexity of setting it up.

Given a signal $x(n)$, you can generate a signal $x_{transformed}$ with:

$$ x_{transformed}(n) = c \times log(1+x(n))$$

Where $c$ is a normalisation parameter. Due to its shape the $log(\cdot)$ function will boost small amplitude values and limit high amplitude values. To get a feel about how the waveform might look like, you might want to have a quick look at the way Audacity does its 'dB View' of waveforms.

You might find both the compression and logarithmic transformation concepts useful to what you are trying to achieve.

Hope this helps.


Transfer function of a nonhomogeneous difference equation Consider the following difference equation:

$y_k=\alpha y_{k-1}+\beta x_k$

The transfer function for this is given by:

$\displaystyle\frac{Y(z)}{X(z)}=\frac{\beta}{1-\alpha z^{-1}} = \frac{\beta}{z-\alpha}$

How does one calculate the transfer function for

$y_k=\alpha y_{k-1}+\beta x_k + \gamma$?

The last difference equation is not a linear system due to the addition of the constant $\gamma$, therefore it does not have a transfer function.


Complex gaussian random variable In my work i need to generate circularly symmetric complex gaussian random variables with non zero mean and certain variance in matlab.

I know the command for generating in case of zero mean , but how to do it for non-zero mean.

Answer and hint: Just add the mean to expression of zero mean random variable. Along with this scale the expression for unit variance gaussian with the standard deviation of the required variance. You would now have non zero mean gaussian RV with given variance