How to Extrapolate a 1D Signal? I have a signal of some length, say 1000 samples. I would like to extend this signal to 5000 samples, sampled at the same rate as the original (i.e., I want to predict what the signal would be if I continued to sample it for a longer period of time). The signal is composed of several sinusoidal components added together.

The method that first came to me was to take the entire FFT, and extend it, but this leaves a very strong discontinuity at frame 1001. I've also considered only using the part of the spectrum near the peaks, and while this seems to improve the signal somewhat, it doesn't seem to me that the phase is guaranteed to be correct. What is the best method for extending this signal?

Here's some MATLAB code showing an idealized method of what I want. Of course, I won't know beforehand that there are exactly 3 sinusoidal components, nor their exact phase and frequency. I want to make sure that the function is continuous, that there isn't a jump as we move to point 501,

vals = 1:50;
signal = 100+5*sin(vals/3.7+.3)+3*sin(vals/1.3+.1)+2*sin(vals/34.7+.7); % This is the measured signal
% Note, the real signal will have noise and not be known exactly.
output_vals = 1:200;
output_signal = 100+5*sin(output_vals/3.7+.3)+3*sin(output_vals/1.3+.1)+2*sin(output_vals/34.7+.7); % This is the output signal

figure;
plot(output_signal);
hold all;
plot(signal);

Basically, given the green line, I want to find the blue line.

I think linear predictive coding (otherwise known as an auto-regressive moving average) is what you are looking for.  LPC extrapolates a time series by first fitting a linear model to the time series, in which each sample is assumed to be a linear combination of previous samples.  After fitting this model to the existing time series, it can be run forward to extrapolate further values while maintaining a stationary(?) power spectrum.

Here is a little example in Matlab, using the lpc function to estimate the LPC coefficients.

N = 150;    % Order of LPC auto-regressive model
P = 500;    % Number of samples in the extrapolated time series
M = 150;    % Point at which to start predicting

t = 1:P;

x = 5*sin(t/3.7+.3)+3*sin(t/1.3+.1)+2*sin(t/34.7+.7); %This is the measured signal

a = lpc(x, N);

y = zeros(1, P);

% fill in the known part of the time series
y(1:M) = x(1:M);

% in reality, you would use `filter` instead of the for-loop
for ii=(M+1):P
    y(ii) = -sum(a(2:end) .* y((ii-1):-1:(ii-N)));
end

plot(t, x, t, y);
l = line(M*[1 1], get(gca, 'ylim'));
set(l, 'color', [0,0,0]);
legend('actual signal', 'extrapolated signal', 'start of extrapolation');

Of course, in real code you would use filter to implement the extrapolation, by using the LPC coefficients a as an IIR filter and pre-loading the known timeseries values into the filter state; something like this:

% Run the initial timeseries through the filter to get the filter state
[~, zf] = filter(-[0 a(2:end)], 1, x(1:M));

% Now use the filter as an IIR to extrapolate
y((M+1):P) = filter([0 0], -a, zeros(1, P-M), zf);

Here is the output:

It does a reasonable job, though the prediction dies off with time for some reason.

I don't actually know much about AR models and would also be curious to learn more.

--

EDIT: @china and @Emre are right, the Burg method appears to work much better than LPC.  Simply by changing lpc to arburg in the above code yields the following results:

The code is available here: https://gist.github.com/2843661


Computes velocity from sampled position My system is :

A step by step wheel returns its position each 150ms. I want to pilot an horizontal infinite scroll bar with the wheel. The scroll bar needs velocity as input to move.
Moreover the faster the wheel rotates the farthest the scroll must move. That means that if I turn the wheel one turn quickly, the bar must go further if I had turned the wheel one turn slowly.

How can I compute the right velocity to give as input to the scroll bar with sampled position of the wheel ?

You have to estimate a velocity based on a regularly sampled signal. The velocity is given by the derivative on a continuous time signal. The approximation of this value with discrete data is a long-lasting topic. @Juancho is correct with the proposal of the $2$-point difference scheme.

As the sampling is never perfect, and $x[i]-x[i-1]$ is not symmetric. So if you do not need fully causal results, there exists more symmetric derivative schemes, such as an average of forward and backward differences:  $\frac{1}{2}(x[i]-x[i-1])+(x[i+1]-x[i]) = \frac{x[i+1]-x[i-1]}{2}$. Other schemes can be found in Numerical differentiation formulae (see here for a discussion).

If you need to compute forward derivatives for real-time purposes, a paper (in French, but formulae are readable) is Estimation par maximum de vraisemblance de la dérivée d'un signal bruité, with cares about noise issues.


Ways to parallelize Sifting in EMD (Empirical Mode Decomposition) Sifting process, which is recursive in nature, is used in EMD (Empirical Mode Decomposition) to decompose the signal (x(t)) into Intrinsic Mode Functions (IMFs) and residue, which is a monotonic function. All the examples provided of Sifting is recursive and all parallel implementation of EMDs where GPU or GPGPU or FPGA are utilized (Example), only use the parallel nature of hardware to speedup the mode decomposition algorithm rather than parallelizing the Siting process instead.

Is it possible to parallelize the Sifting process? If Yes, then do you happen to have any ideas regarding this?

Is it possible to parallelize the Sifting process?

No.

The "problem" here is not the recursion itself but the fact that the data changes at each recursion step.

It is possible to parallelise divide and conquer type of algorithms that imply hierarchy as is the case of the split radix FFT or Quicksort for example. But in that case, the input sequence of values, let's call it $x[n]$ does not change with every step of recursion. Therefore, this ends up being an "exercise" in managing subsequences and the results of function application. Consequently, it is possible to resolve recursion by storing which indices are supposed to be processed at each step of recursion and storing them on a list, effectively transforming recursion to iteration.

Contrast this with empirical mode decomposition where the vector that is sent to the lower (or higher) levels of recursion is not $x[n]$ but $x[n]-q[n]$ where $q[n]$ is the mean waveform that is formed by combining the maxima and minima interpolated wavefoms.

In this way, the maxima and minima of the second level depend on the subtraction of the previous level. It is not simply a matter of finding the "first order" maxima and minima, then the "second order" maxima and minima, then the third order and so on and then recombining them at different levels of recursion because the first level subtraction has the potential to change which $x[n]$ become local maxima (or minima) at a subsequent step.

Therefore, empirical mode decomposition has an inherent sequential part that cannot be parallelised.

Hope this helps.


Window period(overlap) and FFT How does changing the window period (i.e the number of points overlap between two frames) affect the FFT results ?

Suppose that a time series signal was converted to frequency domain by FFT with window size =2048 and window overlap=1024 and the output was plotted as a spectrum. Now will the appearance of the spectrogram change if the window overlap is reduced to 512 ?

In addition to what others have already said, I'll try to answer it from a purely practical point of view (this is also a variant of the overlap-add technique).

If your FFT length is 2048, then an overlap of 1024 (50%) means that you have twice as many analysis (FFT) frames (as compared to the number of frames without any overlapping). A 512 samples overlap (75%) means 4 times as
many frames and so on.

The purpose of overlapping is primarily to reduce the effect of windowing. Most windowing functions (e.g Blackman, Hamming etc) are taper-shaped, which means that they drop to 0 (or close to 0) near the frame edges. This, of course, affects FFT results and
we may lose some important information (e.g. transients).

So to reduce this negative effect of the windowing, we use overlapping. The basic idea here is that we can average FFT results form overlapping frames and thus obtain a better frequency representation of
our time-domain signal. The actual frequency resolution is still the same as without overlapping.

As an example, let's say we use 1024 samples for the FFT and we have a 50% overlap (512 samples). We can average 3 overlapping frames to compute the final (averaged) magnitude values

M_avg[frame][i] = (M[frame-1][i]+M[frame][i]+M[frame+1]]i])/3 ; // i is just a bin index--in this case the bins range from 0 to 512

So, to compute the magnitudes for, say, frame #5, we average magnitudes from frames #4, #5# and #6 and divide by 3 (the number of averaged frames)

There are other use case for overlapping but the one described above is probably used the most.


Sampled signal reconstruction using matlab I have a wav file recorded from my smartphone's mic, and I want to reconstruct the sampled signal and plot the reconstructed signal. After some research and search, I was able to get the following code:

[y, fs, nb] = wavread('new.wav');
x=[1:1:length(y)];
z = y(:,[1]);
samp=transpose(z);
b = 48000/fs;

% Quantizing the Sampled Signal
    no = N/b;           % no. of samples
    z = ceil(5*(no));       % no. of bits required
    quan = zeros(1,z);      % constructing a zero array
    c = (abs(p_max) + abs(p_min))/32     % size of quantized level
    u = p_min;
    w = 1;
    quant = zeros(1,N);
    for (q = 1 : b : N)
        for (r = 1 : 32)
            y = r;
            if (samp(q)>=u & samp(q)<=(u+c))
                for (v = (w+4) : -1 : w) % this loop is for converting decimal into binary
                    quan(v) = mod(y,2);
                    y = (y/2)- (mod(y,2)/2);
                end
                 quant(q) = r;
            end
            u = u + c;
        end
        u = p_min;
        w = w + 5;
    end
    quan

    % Writing the Quantized Signal
    fid = fopen ('quant.bin','w');
    count = fwrite (fid,quant,'int8');

    % Reading the Quantized Signal
    fid = fopen ('quant.bin','r');
    [quant_r,count] = fread (fid,inf,'int8');

    % Reconstructing the Sample of the Quantized Signal
    for (n = 1 : N)
        u = p_min;
        for (ab = 1 : 32)
            if (quant_r(n)==ab)
                quant_r(n) = u;
            end
            u = u + c;
        end
    end
    figure
    stem (quant_r);
    title ('The Reconstructed Sample of Quantized Signal')
    grid on;

    % Reconstructing the Original Signal
    for (c1 = 1 : b : (N-b))
        if (b~=1)
            xyz = (quant_r(c1+b)-quant_r(c1))/b;
            for (c2 = (c1+1) : (c1+b-1))
                quant_r(c2) = quant_r(c2-1) + xyz;
            end
        end
    end
    figure
    quant_r = smooth (quant_r);
    plot (quant_r);
    title ('The Reconstructed Signal')
    grid on;

But it looks like I can't get my desired output, I'm completely new to signal processing, can someone tell me what is wrong with my code? Thank you!

Update#1: I just realized that I used a fairly high sample rate(44100hz), so I tried downsampling using the following code:

[piano,fs]=wavread('new.wav'); % loads piano
play=piano(:,1); % Renames the file as "play"

%t = linspace(0,time,length(play));          % Time vector
x = play;
y = decimate(x,25);

stem(x(1:30)), axis([0 30 -2 2])   % Original signal
title('Original Signal')
figure
stem(y(1:30))                        % Decimated signal
title('Decimated Signal')

%changes the sampling rate

fs1 = fs/2;
fs2 = fs/3;
fs3 = fs/4;
fs4 = fs*2;
fs5 = fs*3;
fs6 = fs*4;

wavwrite(y,fs3,'new1');

Is this what I should do? Thank you!

Update#2: I think my problem is this line:  b = 48000/fs; it should be sampled interval, I have no idea how to deal with this.

Update#3: I just found the following code that seems do the same thing, is this right? If so I don't need to implement it myself.

sub1=wavread('jamming.wav');
 magnSub1=abs(fft(sub1));
 phaseSub1=angle(fft(sub1));
 sub1_L=magnSub1.*exp(i*phaseSub1);
 sub1_L=ifft(sub1_L);
 q=(abs(sub1_L).*cos(angle(sub1_L)));
 plot(q);
 axis([44000,45000,-1,+1]);

I'm not exactly sure what you mean by

  I want to reconstruct the sampled signal

The term reconstruct has a special meaning in DSP and is related to converting a signal from discrete form to continuous using a DAC and a low-pass filter.

I'm going to guess that you are just looking to plot the signal in time domain, and if that is the case then your code would look like this:

[y, fs, nbits] = wavread('signal.wav')
sample_period = 1/fs;                % sample period is just the inverse of the sample rate
sample_time = size(y)/fs;            % total time represented by the samples
time = 0:sample_period:sample_time;  % create the time array
plot(time,y)


When might looking at rising edges perform better than cross-correlation? A simple question as the title indicates:

Let us assume that we have two signals, measured from two different sensors that are close together. (For this application, the sensors are measuring blood/biological activity at different positions on an animal body).

I am wondering when, (if at all), measuring the rising edges of both signals might offer a better estimate of their relative time delays, VS doing a simple cross-correlation?

When might we expect this to be true, if ever?

I ask because I would like to know if looking at rising edges may indeed be optimal in some sense, instead of it being something 'ad-hocy' and possibly mis-leading in lieu of cross-correlation.

Using the rising edge is better than cross-correlation when there is the physical equivalent of an inductor in the system, which will introduce inertia into the system.  That inertia will cause the output to rise gradually, then fall back down.  Even in that case the cross-correlation might be more useful, but if you wanted to know when the input starts to affect the output the rise would be more accurate.

The figure above is illustrative of what I'm talking about.  The inductor (or "inertia introducing object") resists any change in the current flow, so initially it resists the flow causing an exponential rise (as opposed to a sharp rise that matches the presumed non-ideal impulse at the input), and then opposes the current stopping when the input ends, causing the exponential fall at the output.

If the input looks more-or-less like an impulse function, it will correlate most strongly with the output at the output's peak, at time 100.  It's clear, though, that the input starts to affect the output at time 1.  In different situations, either measurement could be the more useful one.


Curved noise floor and ugly unexplained peaks after FFT It has been so long since I studied this at University, so I'll better ask here if I'm doing something really wrong here.

I'm using a USRP frontend for L1 band surveillance. I want a 20MHz bandwidth, but I found through testing that it's better to use a sampling frequency slightly higher than twice that bandwidth (I'm using 41MHz), so I can use a local oscillator to move the DC out of the spectrum. Probably there is some additional error of concept here, but this isn't the problem I'm facing now.

First thing to note: I'm using an FFT over a big chunk of data. Depending on the configuration this size may vary, but for example I'm currently using 8192 samples per FFT. Maybe this is too much and I should compose the whole FFT using many smaller ones, somehow?

In order to reduce noise, I'm performing 100 of these FFTs, summing their outputs and dividing by 100. Probably there's a better method, but so far ok.

After I apply the FFT (with FFTW library), I convert the output to dBW/Hz with the following formula (I'm assuming I have Vpeak):

10 * log10(psd_output / (2 * 50 Ohm)) - 60

Now, if I plot this vector, I get something like this:

The noise floor is getting some curvation (up to 5 dBs). I tried applying a hamming or a blackman-harris window to the FFT output with the hope or removing this effect, but it didn't work.

If I test with a jammer, I can see the peak, but I also see some other minor peaks which are unwanted (and that are not seen in an spectrum analyzer with the same configuration):

What am I missing? These should be simple operations, but I'm out of ideas at this point.

EDIT: Adding other plot

Does the 'same configuration of your analyzer' mean the 'same sampling rate'? Is that equal to 41 MHz? It looks like you faced a well-known problem with a Nyquist criterion practical implementation. I mean that some time ago in digital oscilloscope ads there was an exact Nyquist ratio between a sampling rate and channel bandwidth (especially for cheap ones). Now it is about 10. (You may want to look at any ads of digital oscilloscopes) I remember there were several articles with figures similar to those posted here. To check this hypothesis just check a sampling rate of your spectrum analyzer (you wrote that artificial pikes are not seen in a spectrum analyzer with the same configuration.) or try to use a higher sampling rate in your data processing procedure if possible. In order to get at least an intuitive idea of the reason lies in the foundation of the problem perhaps you want to check this: Why is it a bad idea to filter by zeroing out FFT bins?
I suppose the following paragraph from the cited answer works for you: 'So if your original FFT input data is a window on any data that is somewhat non-periodic in that window (e.g. most non-synchronously sampled "real world" signals), then those particular artefacts will be produced' This phenomenon is called resolution bias error, or more commonly, the picket fence effect http://www.azimadli.com/vibman/thepicketfenceeffect.htm. Here you can find an illustration http://www.mechanicalvibration.com/More_on_spectral_leakge.html


Impulse and step response of a system that contains cascaded second order sections I have to design an 8th-order IIR filter with the following transfer function magnitude response as a specification:

I use "ellip" function from matlab to calculate the unquantised coefficient set:

Then i implement the filter using a cascade of second-order sections using "tf2sos" matlab function:

And then i quantize these coefficients:

After all these steps i take the two matrices below:

I choose the filter representation for every second-order instance to be like below (Direct Form 2):

circle with "X" is a multiplication
circle with "Σ" is an addition
and the triangle is a delay element

After that i have done all of these i write some matlab code to model the structure above:

%Implementation of a second-order filter structure
function[y,state,pointer]=Second_Order(x,a,b,state,pointer)

%The filter routine when fed with an input sample produce one output
%sample each time it is called

%x=input value
%a=array containing feedback coefficients[1,a1,a2]
%b=array containing feed forward coefficients[b0,b1,b2]
%state=array(of length 2) containing old input values
%y=output value

%compute output value
y=b(1)*x+state(pointer+1);
pointer=rem(pointer+1,2); %update pointer in modulo form

%Update states
state(pointer+1)=b(2)*x-a(2)*y+state(pointer+1);
pointer=rem(pointer+1,2);

%Overwrite oldest sum with b(N-1).x
state(pointer+1)=b(3)*x-a(3)*y;
pointer=rem(pointer+1,2);%Increment pointer modulo-(N-1)

At the end, and the point that i don't know how to do it properly is how i can display the impulse and step response of the whole filter and not for every second order section.

For example if i write something like that, i think i take the impulse response of the first second order section. How i can take the impulse resonse of the whole system?

%Initialization
[b,a]=ellip(4,1.3,60,2*[0.2,0.3]);
[sos,q]=tf2sos(b,a);
sosq=Quantize('round',sos,2^6);

bf=sosq(1,1:3);
af=sosq(1,4:6);

%impulse invariant design
impsv=zeros([1,2]);%impulse response state variable
stpsv=zeros([1,2]);%step response state variable
len=20; % run length
impy=zeros([1,len]);%impulse response result
stpy=zeros([1,len]);% step response result
impp=0;
stpp=0;

% Simulation of filter operation

x=1;
for k=1:len
  [impy(k),impsv,impp]=Second_Order(x,af,bf,impsv,impp);
  [stpy(k),stpsv,stpp]=Second_Order(1,af,bf,stpsv,stpp);
  x=0;
end

Not sure I follow completely, but is this the question: how to obtain the impulse and step response of a filter that is realized as a cascade of second or sections? or are you asking for working code that implements this?

Your 8th order filter is decomposed into 4 second order sections. The ordering is not really important here. For any input $x$ to your 8th order filter you compute the output $y$ by first filtering the input with the first second order section $SOS_1$ to get output $s_1$, that is $s_1 = SOS_1*x$ (the $*$ means processing the $x$ with a second order difference equation with coefs from $SOS$). Then that output is made input to the second SOS, that is $s_2 = SOS_2*s_1$. So now you have $s_2 = SOS_2*(SOS_1*x)$. And then continue this pattern (output from one section is input to the next) through all sections.

If you understand this, then it should also be rather obvious that you obtain the impulse response or the step response by letting an impulse or a step waveform be input to your filter and the output $y$ contains the response.

I have not looked at your code. But try to start with a simpler case first just to verify your code and your understanding of the concept.


How is the topo / scalp map from an EEG Signal with ICA created? I am lately experimenting with an EEG and found this Library: MNE Python.
I am using an ICA to find and correct artifacts - this is working fine.

It's now possible to plot a topo / scalp map for the ICA (EEGLab and other also can do this).

Now I wonder how this is calculated?

I assume it takes the topographic information of the EEG sensors and combines them with the ICs for each channel. Is this correct?
Are there any sources explaining the process in general?

Thanks in advance

In general, the topographic map describes the distribution of electrical activity across the brain as measured in a number of well known "stations" (electrodes whose position is known) on the scalp.

This electrical activity is the combined result of millions of neurons firing below the skull. The strength of the electrical signal they create is measured across the REFerence electrode and one of the scalp electrodes. In a typical Electroencephalography (EEG) montage, you may have something like 16 to 20 channels (more in special cases).

The Inverse Problem in this case, that is, given the 20 (for example) signals of the EEG, calculate which neurons fired and how is ill-posed because it has a huge number of unknowns (millions of neurons) and only a handful of knowns (20 signals).

What can be done however, is to use Independent Component Analysis (ICA) to decompose the combined electrical activity of millions of neurons in the brain to a small number of independent (in a very specific mathematical sense) components.

In this case, the physical interpretation of this process is that the 20 channel activity is produced by 3,4,5,... independent "sources" within the skull. And, depending on the cognitive task that the subject was going through during data acquisition (whose performance is modulated by their mental state), these independent components might be able to provide some additional information. For example, maybe disease is reducing the number of sources or "changing their location" and so on.

OK, so, what we are really saying in this case is that we are recording the activity of a small number of components within the brain with 20 electrodes at the surface of the brain. We enter all these 20 time series in ICA, and ICA returns:

The Independent Components (the elementary waveforms whose linear combinations generate all 20 original recordings)
The "loadings" for each electrode. Or in other words, the mixing levels, how much of each of the elementary components is required to reconstruct the signal at each electrode.

Essentially now we have the signals of the sources and we are re-constructing the signals at the electrodes.

Since we are now in full control, let's think about what happens with just one elementary waveform, let's call it "IC1". If we now plot the loadings on the positions of the electrodes (and interpolate in between) we get a topogaphic map of the strength of the Independent Component across the surface of the brain and to an extent, we can localise its activity. The same can be done for the rest of the components and this is how these topographic maps are calculated.

In the case of EEG and assuming a "standard" montage, it is possible to generate these general images. In the case of Magnetoencephalography (MEG), you can get the actual positions of the electrodes from the manufacturer to generate these plots accurately.

For more information please see this,this and this papers while for a more detailed overview of the method please see this link. The papers are not "de-facto" papers, but give a good overview of the use of ICA in EEG, both for preprocessing but also for analysis.

Hope this helps.


GnuRadio: Clock Recovery MM: imu out of bounds I have created a flowgraph, based on an example of Michael Ossmann, which takes in a signal, and should output bits.

I need to use the clock recovery MM block, which I do not fully understand yet. However, after reading some blogposts, I am quite sure that I can leave most of the settings default, except for the Omega one. Here's my flowgraph:

When running this flowgraph, it gives me the following error:

thread[thread-per-block[5]:<block clock_recovery_mm_ff (9)>]: mmse_fir_interpolator_ff: imu out of bounds.

While searching, I stumbled upon this piece of code in the source of GnuRadio:

int imu = (int)rint(mu * NSTEPS);
  if((imu < 0) || (imu > NSTEPS)) {
    throw std::runtime_error("mmse_fir_interpolator_ff: imu out of bounds.\n");
  }

So, I suspect it is not due to my Omega setting (which might be wrong, I have to play with that setting), but that it is due to my Mu setting, which is just the default (0.5). However, I understand that Mu needs to be between 0 and 1, so I do not really understand what the problem is. Anyone who does?

Environment details:

GNU Radio Companion 3.7.7.1
Running a GNU Radio live DVD in a virtual machine (VirtualBox 4.2.12) on Windows 7.
Using Volk machine: ssse3_64

The M&M clock recovery block has a bug in it, such that it can't reliably handle more than 8 samples per symbol without inadvertently dropping samples.  That interpolating resampler alignment bug may also be responsible for generating mu values outside of [0.0, 1.0].

I'd strongly recommend, downsampling the input from 2.5ksamples/symbol down to 6 or 8 samples /symbol before the M&M block.


Inverse Discrete Time Fourier Transform of a given Function using the definition I would like to compute the Inverse Fourier Transform of a given function using the DTFT definition

$$H(e^{jω})=\frac{1}{(2−e^{jω})(2−e^{−jω})}$$

But I have to use the definition of the DTFT
$$H(e^{j\omega}) = \sum_{n = -\infty}^{n = \infty}{h(n) e^{j\omega n}}$$

My approach was that I assume that $x=e^{jω}$, which then gives the equation

$$H(x) = \frac{1}{(2-x)(2-x^{-1})}$$
Now I can apply partial fraction decomposition.

Is my approach correct?

Yes, you can do that. Because we substitute $e^{j\omega}$ for $z$ to get the frequency response of a discrete system, you can go back by substituting $z$ for $e^{j\omega}$. This yields

$$H(z) = \frac{1}{(2-z)(2-z^{-1})} = \frac{1}{5-2z^{-1}-2z} = \frac{z^{-1}}{-2+5z^{-1}-2z^{-2}}$$

Therefore the system you're trying to find is

$$-2y[n]+5y[n-1]-2y[n-2] = -x[n]$$

$$y[n] = \frac{-x[n]-5y[n-1]+2y[n-2]}{-2} = \frac{1}{2}x[n] + \frac{5}{2}y[n-1] - y[n-2]$$


What is the rule for manipulating the boundaries of a summation? When working with DTFTs or $\mathcal Z$-transforms, we sometime get summations that do not go from $n=0$ to $+\infty$.
For example, suppose we have the sequence $x(n) = -\alpha^n u(-n-1)$. To find the $\mathcal Z$-transform:
\begin{align}
X(z)&=\sum_{n=-\infty}^\infty x(n)z^{-n}\\
&=-\sum_{n=-\infty}^{-1} \alpha^nz^{-n}
\end{align}

I understand this, the step function(which was shifted and inverted) causes the sequence to exist for all values from $-\infty$ to $-1$, and is $0$ everywhere else. However, the next part of the solution confuses me:
$$=-\sum_{n=0}^{\infty} (\alpha^{-1}z)^{n+1}  $$
I am not sure how they managed to change the boundaries.

From my interpretation of the result , a flip of the boundaries $\sum_{n=-\infty}^{-1}$ to $\sum_{n=-1}^{\infty}$ for any sequence will cause the sign of $n$ to change. Is this true?
Then, to get from $\sum_{n=-1}^{\infty}$ to $\sum_{n=0}^{\infty}$,do you just need to add $1$ to wherever here is an $n$?
Should you always first flip the boundary (when necessary) and then only change the index afterwards?

Summations are not like integrals. When you calculate an integral, swapping the boundaries does change the sign of the result. However, when talking about summations:

$$ \sum_{n=-\infty}^{-1}a(n)=\sum_{n=-1}^{-\infty}a(n)$$

This holds because it is the same if we perform the sum of $a(-1)+a(-2)+...+a(-k)$ for any $k>0$ in any order. Namely, you can order those terms as you wish as long as all the terms are present. For infinite summations, this is true if the series is absolutely convergent.

In your example, they just made the change of variables $m=-n-1$ to get the summation from $0$ to $+\infty$. It may be confusing because they kept the original name of the index ($n$) instad of replacing it with the new one ($m$), but remember that they are just indexes or dummy variables, so in the end it doesn't matter. You can check that if $m=-n-1=-(n+1)$, then the boundaries change:

$$n=-1 \to m=-(-1+1)=0$$
$$n=-\infty \to m=-(-\infty)=\infty$$

where I used the fact that a finite number is negligible against infinity. So the summation can be written as:

$$\begin{align} -\sum_{n=-\infty}^{-1} \alpha^nz^{-n}&=-\sum_{m=\infty}^{0} \alpha^{-(m+1)}z^{-(-(m+1))}\\
&=-\sum_{m=0}^{\infty} \alpha^{-(m+1)}z^{-(-(m+1))}\\
&=-\sum_{m=0}^{\infty} \alpha^{-(m+1)}z^{m+1}\\
&=-\sum_{m=0}^{\infty} (\alpha^{-1}z)^{m+1}
\end{align}$$

And that's the same result you wrote in your question.


Costa's loop questions, design, tracking and limitations Why does Costa's loop PLL bandwidth need to be narrow? and why can't it track a certain frequency offset? and does that hold even if it has a filter with a pole at DC?

Costa's loop is related to traditional PLLs. One way of characterizing a PLL is by its loop bandwidth and damping factor. The damping factor determines the loop response (under-damped, over-damped or critically damped). The loop bandwidth determines how fast a PLL/Costas will achieve a "lock". The higher the loop bandwidth, the shorter the acquisition time. However, higher loop bandwidths result in higher tracking error (and hence higher locking threshold SNR - See graph below). Therefore, the loop bandwidth is a tradeoff between acquisition time and tracking error.  In some systems, the PLL/Costas start with a higher loop bandwidth before acquisition. After acquisition, the loop bandwidth is reduced to minimize tracking error.

The following terms associated with general PLL can answer your second question

Pull-in range: is the maximum initial frequency difference between the input
and VCO center frequencies both in positive and negative directions, for which
the PLL eventually achieves the phase-locked condition. The pull-in range is
related to the dynamics of the PLL
Lock-in range: is the frequency range over which the PLL achieves the
phase-locked condition without cycle slips

A little Googling can help you to find the specific formulas of the two terms above for the Costas Loop.

Regards,

Moses.


How can I graph a simulation of a FM 16 bits codification? I have to develop a software where the user inputs a 16 bits binary sequence, which is supposed to be the codification of a FM (or even AM) signal, then I should graph it according to Frequence / Time ....

What procedure I should take to treat such input?

It's hard to tell what you're asking here, but it sounds like you have a discrete-time frequency-modulated signal with 16 bits per sample. In order to convert this to a frequency-versus-time signal, you want to demodulate the FM. One way to do this is with an FM discriminator. A digital implementation could look like this (assuming a complex input signal $x[n]$ sampled at a rate of $f_s$ Hz):

$$
x_{d}[n] = x[n]x^*[n-1]
$$

$$
f[n] = \text{atan2}\left(\text{Im}[x_d[n]], \text{Re}[x_d[n]]\right) \frac{f_s}{2\pi}
$$

This simple implementation estimates the frequency of the signal $x[n]$ using a first-order phase difference. The phase of the product $x_d[n]$ is equal to the phase difference between successive samples $x[n]$ and $x[n-1]$. This phase difference can be used as an approximation of the frequency of the signal at time $n$. To extract the phase of $x_d[n]$, we use the atan2 function, then scale it appropriately so that $f[n]$ is measured in Hz. $f[n]$ is then the estimated frequency modulation input.

Note: As I said before, this technique assumes a complex input signal. If what you have is real-valued, you'll want to convert $x[n]$ to analytic signal format first.


What is the meaning of positive and negative shift in signal correlation? If $x(t)$ and $y(t)$ are real valued signal then from this page cross correlation between them is defined as

\begin{align}
R_{XY}(\tau) &= \int_{-\infty}^{\infty}x(t)y(t-\tau) dt\quad\text{+ve shift}\tag{1a}\\
R_{XY}(\tau) &= \int_{-\infty}^{\infty}x(t+\tau)y(t) dt\quad\text{-ve shift}\tag{1b}\\
\end{align}

I tried to solve a problem where $x(t) = e^{-at}u(t)$ and $y(t) = e^{-bt}u(t)$.

  From $(1a)$: Positive shift

\begin{align}
R_{XY}(\tau) &= \int_{-\infty}^{\infty}x(t)y(t-\tau)dt\\
&= \int_{-\infty}^{\infty}e^{-at}e^{-b(t-\tau)}u(t)u(t-\tau)dt\tag{2}
\end{align}
For $\tau\geq0$. Limits of integration will be from $\tau\leq t\leq\infty$. Hence from $(2)$ after simplification integration can be written as
\begin{align}
R_{XY}(\tau) &= \int_{\tau}^{\infty}e^{-at}e^{-b(t-\tau)}dt\\
&=e^{b\tau}\int_{\tau}^{\infty}e^{-(a+b)t}dt\\
&=\frac{e^{-a\tau}}{a+b}
\end{align}
For $\tau\leq0$. Limits of integration will be from $0\leq t\leq\infty$. Hence from $(2)$ after simplification integration can be written as
\begin{align}
R_{XY}(\tau) &= \int_{0}^{\infty}e^{-at}e^{-b(t-\tau)}dt\\
&=-\frac{e^{b\tau}}{a+b}\int_{0}^{\infty}e^{-(a+b)t}dt\\
&=\frac{e^{b\tau}}{a+b}
\end{align}
Therefore
\begin{align}
R_{XY}(\tau)=\begin{cases}
                    \frac{e^{-a\tau}}{a+b}; \tau\geq0\\
                    \frac{e^{b\tau}}{a+b}; \tau\leq0
             \end{cases}\tag{3}
\end{align}

  From $(1b)$: Negative shift

\begin{align}
R_{XY}(\tau) &= \int_{-\infty}^{\infty}x(t+\tau)y(t)dt\\
&= \int_{-\infty}^{\infty}e^{-a(t+\tau)}e^{-bt}u(t+\tau)u(t)dt\tag{4}
\end{align}
For $\tau\geq0$. Limits of integration will be from $0\leq t\leq\infty$. Hence from $(4)$ after simplification integration can be written as
\begin{align}
R_{XY}(\tau) &= \int_{0}^{\infty}e^{-a(t+\tau)}e^{-bt)}dt\\
&=e^{-a\tau}\int_{0}^{\infty}e^{-(a+b)t}dt\\
&=\frac{e^{-a\tau}}{a+b}
\end{align}

Corrected my analysis after @Fat32 answer

For $\tau\leq0$. Limits of integration will be from $-\tau\leq t\leq\infty$. Hence from $(4)$ after simplification integration can be written as
\begin{align}
R_{XY}(\tau) &= \int_{-\tau}^{\infty}e^{-a(t+\tau)}e^{-bt)}dt\\
&=e^{-a\tau}\int_{-\tau}^{\infty}e^{-(a+b)t}dt\\
&=\frac{e^{b\tau}}{a+b}
\end{align}
Therefore
\begin{align}
R_{XY}(\tau)=\begin{cases}
                    \frac{e^{-a\tau}}{a+b}; \tau\geq0\\
                    \frac{e^{b\tau}}{a+b}; \tau\leq0
             \end{cases}\tag{5}
\end{align}

Hence from $(3)$ and $(5)$ as can be seen the result match. If my analysis is wrong do correct me and if its correct then what does these negative and postive shift signifies?

On substituting $u = t - \tau \Rightarrow t = u + \tau \Rightarrow dt = du$ in $(1a)$, we get $(1b)$ as
\begin{align}
R_{XY}(\tau) &= \int_{-\infty}^{\infty}x(u+\tau)y(u) du = \int_{-\infty}^{\infty}x(t+\tau)y(t) dt\\
\end{align}

Hence

\begin{align}
R_{XY}(\tau) &= \int_{-\infty}^{\infty}x(t)y(t-\tau) du = \int_{-\infty}^{\infty}x(t+\tau)y(t) dt\\
\end{align}

Your Eq(3), for positive shift, is correct. And your Eq(5), for negative shift, should also be the same, but is wrong for $\tau < 0$ case.

The lower limit of the integral in Eq(4) should be $-\tau$ but you took it $\tau$. The reason is the step function $u(t+\tau)$ will have the nonzero range for $t > -\tau$.

You will get the same result with the correct limit.


Rayleigh fading with frequency selective fading channel I Read in one reference that, the Rayleigh distribution is commonly used to describe the statistical time varying nature of the received envelope of a flat fading signal, or the envelope of an individual multipath component. Can Rayleigh and Rician distributions be used with frequency selective fading channels? I need some explanation please

First it is useful to understand what a Rayleigh and Ricean distribution is, and then we can more easily see how that statistic can be used to describe certain fading channels.

A zero-mean complex signal that has an identical Gaussian distribution on both the in-phase and quadrature components will have a magnitude that is Rayleigh distributed. A useful explanation I came across that explained this well is an Agilent (HP) app-note on spectrum analyzer noise.

https://literature.cdn.keysight.com/litweb/pdf/5966-4008E.pdf

And I will paste the relevant figure below, thanks the HP/Agilent/Keysight. The first plot shows typical samples of a complex zero mean Gaussian signal. If we drew the histogram of the distribution on the real and imaginary axis, both would be the familiar bell curve of a Gaussian distribution.

If the complex signal with an identical Gaussian distribution on both the in-phase and quadrature components was not of zero-mean, such as I depict in the figure below, then the amplitude distribution of this process would be Ricean.

So how does this apply to fading channels? The central limit theorem comes into play here. The central limit theorem states that the sum of a large number of independent random variables will have a distribution that tends toward a Gaussian distribution regardless of the distribution of the individual processes.

In the case of fading channels we are seeing the result of many independent signal paths from the transmitter reaching the receiver at different delays. The different arrival times results in each signal being an independent random process. Often in terrestrial land-mobile communications there is no particularly strong direct-path, but a distribution of many reflected paths resulting in a zero mean Gaussian complex signal at the receiver (once the carrier is removed). To the extent that the signal actually has this process (which means to the extent that there are many indirect paths and no particularly single dominant direct path), then the AMPLITUDE or related power level of this received will be well modeled by a Rayleigh distribution.

When there is a dominant signal path present, particularly stronger than the others, with the presence of multiple reflected paths, then the amplitude of the received signal will follow a Ricean distribution. This is more common in satellite communications for example.

The above occurs whether we are dealing with flat fading (the amplitude across the entire signal bandwidth changes as one), or frequency selective fading (the amplitude at particular frequencies within the channel bandwidth is fading independent of the amplitude of other frequencies within the channel bandwidth.

Frequency Selective or Flat Fading

To be frequency selective fading or flat fading is independent of the distribution of the signal strength being Rayleigh. That only has to do with delay spread. (See figure below).

The Delay Spread of the channel is what determines the frequency selectivity of the fading (regardless if it is Rayleigh or not). Typical delay spreads for the macrocellular environment are shown in the figure below:


Convert 96 Khz to 48 Khz audio: is this simple downsampling method ok? I have an audio signal x sampled at 96 Khz. I have already lowpassed it with a cutoff somewhere between 20 Khz and 22 Khz, so there should be nothing left of frequency > 22 Khz.

I need to resample it to 48 Khz. Since we are in the specific case 96000=2*48000, is it ok to do it this way:

y[n] = x[2*n]

by just forgetting every two elements?

What are the drawbacks of this method, in comparison to more complex downsampling algorithms? (Since 96 Khz is an integer multiple of 48 Khz, I thought that there would be no risk of aliasing here).

Also, what would happen with a downsampling y[n] = x[n * 2] if I had not low pass filtered at 20 Khz? What could be the potential problems?

Since you have already properly low pass filtered the signal, then there is no risk in taking every other sample to complete the downsampling operation.  If you were to create a low pass filter that passed your spectrum of interest with no distortion, and rejected all energy in the alias frequency locations then taking every other sample would provide a perfect downsampling by a factor of 2 (there would be no distortion in you signal of interest). Of course, such a perfect filter is not feasible to implement but stated as such to show that the issue is not in taking every other sample but entirely in the filter design itself.


Spectral Envelope from MFCC I'm obtaining MFCC by using the function wave2mfcc from the Speech and Audio Processing library (SAP) from here http://mirlab.org/jang/matlab/toolbox/

The steps follows the procedure:

Take the Fourier transform of (a windowed excerpt of) a signal.
Map the powers of the spectrum obtained above onto the mel scale, using triangular overlapping windows.
Take the logs of the powers at each of the mel frequencies
Take the discrete cosine transform of the list of mel log powers, as if it were a signal

As I have previously read, the MFCC are the spectral envelope of the signal, but I don't know how to relate the coefficients with the envelope.

For example, if I have a signal 'example.wav'

[signal fs] = wavread('example.wav');
    figure(1); plot(log(abs(fft(signal)));
    mfcc = wave3mfcc(signal,fs);
    figure(2); plot(mfcc):

As I understand the MFCC should have the shape of the figure 1, but in my case it doesn't, maybe I'm missing some step.

Can anybody tell me please how to relate them and being able to graph?

It would be more accurate to say that the MFCC "capture" or "represent" the spectral envelope - in the sense that two signals with a similar spectral envelope will have a similar sequence of MFCC.

The MFCC coefficients themselves are not supposed to "look" like the spectral envelope when plotted. After all the spectral envelope is a "big" vector (as long as the size of an analysis window), while the MFCC are designed to be short (a handful of coefficients). The spectral envelope is highly redundant (this is a very smooth function, consecutive values are correlated), while the MFCC are, by design, very little correlated.

One could visualize the spectral envelope from the MFCC by inverting all steps in the process except the first one (note that you probably omitted one last step: truncate the vector and keep first 13 or 20 coefficients). The reverse operations would thus be: pad with zeros, compute inverse DCT, take exponential, and compute a matrix product with the matrix containing the triangular bin responses. This will yield a smoothed version of the spectrum that will show what the MFCC "see" of the spectrum.

Finally, unless "example.wav" is very short, there is something wrong in your code. Audio signals are analyzed in short overlapping frames rather than wholly. It's likely that you might instead want to look at spectrograms; and at the MFCC as a 2D matrix (imagesc).


Why sine wave leakage in FFT spectrum The input is three sine waves with different amplitudes and frequencies. After FFT the spectrum shows the correct characteristic of each wave without leakage error. Why does the other spectrum leakage when the input is multiplied by a Hamming window function? Thanks]1

Yes, we use windows to reduce spectral leakage and you show an example where using a window increases it.

There is more than one way to explain it.

If you express a DFT as a matrix vector product

$$ \mathbf{y}=\mathbf{Wx}
$$
where
$$
\mathbf{WW^H}=c\mathbf{I}
$$
where $c$ is a constant that depends on how you define your DFT. The point is that $\mathbf{W}$ is an orthogonal matrix. Your 3 original sine waves correspond to 3 rows of $\mathbf{W}$. Your 3 sines are perpendicular to all the other $N-3$ rows, so there is no leakage.

When you uses a window like a Hamming Window, it’s equivalent to multiplying each row element by element by the window. This modified matrix $\mathbf{\tilde{W}}$ isn’t orthogonal any more.
$$
\mathbf{\tilde{W}\tilde{W}^H}\ne c \mathbf{I}
$$
The resulting matrix is nearly equal to a constant times $\mathbf{I}$ but the off diagonals are nonzero.  Your 3 original sines now project on all the rows of $\mathbf{\tilde{W}}$. These nonzero projections are another way to say leakage.

The nonorthogonality of the widowed matrix is easy to verify in matlab using the dftmtx command.


Correct SNIR computation units I would like to know the correct units to make a SNIR computation.

I am using the following formula:

$$
\text{SNIR} = \frac{1}{\frac{1}{\text{SNR}} + \frac{1}{\text{CIR}}}
$$

I have SNR in dB and CIR in dBm.

Shall I convert SNR to dBm to get the SNIR in dBm, or the CIR to dB to get the SNIR in dB, or there is no need to make the units change?.

I want to use the SNIR to compute a BER which is $\approx \frac{1}{4\textrm{SNR}}$

Regards

You'll have to convert back into linear units before you can evaluate the SINR formula.

${\textrm{SNR}_{linear}=10^{\big(\frac{\textrm{SNR}_{dB}}{10}\big)}}$ and likewise for CIR...

${\textrm{SINR}_{dB}  = 10 \log_{10}(\frac{1}{\frac{1}{\textrm{SNR}}+\frac{1}{CIR}}})$

${BER = \frac{1}{4{\textrm{SINR}}_{linear}}}$


10 bit pixels convert to RGB I have a camera which is providing me raw, binary data with 10 bit data. It follows the Bayer pattern RGGB, and each pixel is saved using 16 bit chunks. I am able to demosaic the image almost successfully.

The problem

The RGB image looks is almost completely black. If I multiply by 64 the raw data, the RGB image looks better but rather green. The overall intensity is very low. The camera provider has a live view of the camera and it looks much better in both aspects.

My guess

Since 10bit are saved into 16bit chunks, the tools I am using assume 16bit depth (at least Matlab and OpenCV do for sure). In order to compensate for this, I multiply by 64 (6bit shift to the left) the raw data. But still is not enough. If I multiply by 128, some pixels saturate and then the image looks weird.

What I have tried

Using Matlab image processing toolbox, Python colour-demosaicing package, and OpenCV functions, I get the same result, but they do not match the live visualization.

Update
I analyzed the histogram on the vendor tool, and those from OpenCV. I am not sure, but interestingly there are no values above 255. However, on the raw data I record from their software, there are many values above 255, but nothing compared with values below 255. Thus, if I clip the image with max value of 255, the final result looks very similar!.

Questions

Is there anything else I can verify? I assume is a problem on my end, because the three tools give same results.
Is there a better way of handling the 10bit pixels to convert to 16bit?
Is there any sense on what I discovered on the update? It seems that the two extra bits are not used for anything but noise!

Any help is appreciated. I am very sorry for not providing data or images (company privacy).

Thanks!

The raw data you get from the camera is nearly linear in the colorspace spanned by the camera CFA filters. The 8-bit/sample monitor you have assumes the image is gamma-compressed, most likely within the sRGB colorspace. These two colorspaces are different, therefore you need to map between the two in some way. Understanding this, lets address each of your issues:

  The RGB image looks is almost completely black.

The maximum 10-bit value is 1023, the maximum 16-bit value is 65535. To Convert 10-bit linear to 16-bit linear you multiply by 65535./1023 and round to nearest integer. Alternatively it might be easier to work with floating point directly, so just divide by 1023 to map your data to the [0,1] range.

  the RGB image looks ... rather green

It is common for the camera sensors to be more sensitive to green, so the green response is higher than red and blue, giving an overall green tint. To get physically correct measurements of the incoming light within a known colorspace (say the ITU-R BT.709 primaries used by your monitor) you need to figure out the color transformation matrix from the camera colorspace to your preferred colorspace. There are different ways to accomplish this: get a color-calibration chart, contact the sensor manufacturer, or measure the exact spectral response using other means (like a monochromator).

That being said, if you don't seek for physically correct measurements but rather just for a good-enough picture, you might skip the above calibration step and instead do only a white-balancing step: find an ought-to-be gray object in the frame and divide each pixel by the RGB of that object, thus cancelling the green tint. This on its own, without a color transform, will not reproduce colors exactly though.

  The overall intensity is very low

This is due to the linear versus gamma discrepancy. You do all your math in linear color space with high bitness, and as a last step, convert it to 8-bit sRGB for display. See sRGB for the exact formula.


orthogonal family and pulse shaping filter I am a beginner in digital communications and I have a question on the pulse shaping filter.

For example, I have an orthogonal family $ \{s_1,...,s_N \}$. Can we preserve the orthogonality with a pulse shaping filter $g$ as a raised cosine filter, that is to say the family $ \{ g(n) \ast s_1(n),...,g(n) \ast s_N(n) \}$ is orthogonal?

In general for a pulse shaping filter and any orthogonal code family: No.

UPDATE: See @Engineer's answer which I confirmed by excluding the first row in the set of Walsh codes (the all 1 row) the remaining rows which have a balance of 1's and 0's are indeed all orthogonal.

*SECOND UPDATE: @MBaz has pointed out here that with Root-Raised Cosine filtering the entire code set will remain orthogonal: Orthogonal Codes for Band Limited Channel

A detailed mathematical proof along with the criteria where orthogonality can be maintained would be a much better answer, but since no one has submitted one I am at least providing here a simulation result that demonstrates a specific example of an orthogonal family where after pulse shaping with a raised cosine filter the orthogonality was not maintained. (The original signal can of course be extracted at the sampling locations with no distortion given the zero-ISI property of a raised-cosine filter, so that extracted signal would still be orthogonal after sampling, but the question was specifically on the orthogonality of the pulse shaped waveform itself ($g(n)*s_n(n)$)

The formal definition for orthogonal signals is that their inner product (also called the dot product) is zero.

Consider a simple case of the Hadarmard $H_4$ matrix, with each row forming a set of four orthogonal Walsh codes:

$\begin{bmatrix}1 & 1 & 1 & 1\\1 & -1 & 1 & -1\\1 & 1 & -1 & -1\\1 & -1 & -1 & 1\end{bmatrix}$

Before pulse shaping it is easy to confirm orthogonality, for example between row 1 and row 2 using the dot product of index by index product and sum:

$$(1\times1)+(1\times-1)+(1\times1)+(1\times-1) = 0$$

And likewise betweeen row 1 and row 4:
$$(1\times1)+(1\times-1)+(1\times-1)+(1\times1) = 0$$

The result of a simulation of this set with a raised cosine filter suggests that orthogonality is maintained in the case of the first two rows, but is not between row 1 and row 4. The resulting waveforms from the first two rows after upsampling by 10 and pulse shaping with a raised cosine filter are shown in the figure below (the sample locations for decoding the original signal are at samples 50, 60, 70 and 80):

The computed dot product for these two waveforms was $4.4E-16$. (Close enough to zero so orthogonal).

Row 1 vs Row 3 is similarly quite low at $9.1E-16$ (again orthogonal).

However
Row 1 vs Row 4 was $-0.465$ which is a very strong correlation and clearly no longer orthogonal.

The plot showing the pulse-shaped Row 1 and Row 4 is shown below.

The dot product is the sample by sample product and then the sum (or integration) of the result. Plot of the product of each of the waveform pairs gives further insight into the loss of orthogonality. Note how the product of the pulse-shaped row 1 and row 2 is very symmetric with similar positive and negative areas under the curve, while the product or the pulse shaped row 1 and row 4 has a significantly larger negative area, thus destroying orthogonality in that case:


CLAHE - Redistributing excess pixels among histogram bins In order to clip an image's histogram to a maximum value, according to the CLAHE algorithm, the excess pixels (above the clip limit) will be distributed equally among all the histogram bins (While keeping track none of the bins exceed the max limit of course).

Now, when working with (8,8) region size and distributing 13 pixels on 256 bins, the uniform added value (13/256) will be rounded to zero and so the total number of pixels in histogram will be reduced from (8,8) -> 64px to 64 - 13 = 51px.

That doesn't make sense since I'll be neglecting 13 pixels of the region.

Any tips on how to tackle this part of the algorithm?

you have to redistributed 13 pixels(# of clipped) to the histogram.

 // clip histogram
  194
  195             if (clipLimit_ > 0)
  196             {
  197                 // how many pixels were clipped
  198                 int clipped = 0;
  199                 for (int i = 0; i < histSize; ++i)
  200                 {
  201                     if (tileHist[i] > clipLimit_)
  202                     {
  203                         clipped += tileHist[i] - clipLimit_;
  204                         tileHist[i] = clipLimit_;
  205                     }
  206                 }
  207
  208                 // redistribute clipped pixels
  209                 int redistBatch = clipped / histSize;
  210                 int residual = clipped - redistBatch * histSize;
  211
  212                 for (int i = 0; i < histSize; ++i)
  213                     tileHist[i] += redistBatch;
  214
  215                 for (int i = 0; i < residual; ++i)
  216                     tileHist[i]++;
  217             }


How to calculate the quantum of the difference between two value of clock counter where the values are divided by 64? I have a counter in my embedded system with f = 64 Mhz.

let's say its read value is $cnt_i$
where $i$ is the index of the time I sampled it.
$$\Delta = cnt_{i+1} -cnt_i$$

let's say I take $\Delta^{'} = \Delta/64 = \Delta >>6$ where $>>$ is shift right operator.

I understand that the value of the actual counter now is equivalent to 1 Mhz but How can I show that.

$\Delta^{'} = \Delta/64 $ $[\Delta^{'}]_{sec} = \frac{\Delta}{64Mhz}*\frac{1}{64}$

but it is exactly the opposite, what do I miss?

In order to show that a counter is running at a specific frequency, you need to take two readings spaced at some known time period apart and then calculate how many ticks would occur in one second. Say that you take two measurements from a monotonically-increasing counter at times $t_1$ and $t_2$. Their difference is:

$$
\Delta C = C_{t_2} - C_{t_1}
$$

You can calculate the frequency of the counter from the difference as follows:

$$
f = \frac{\text{number of ticks}}{\text{time period}} = \frac{\Delta C}{\Delta t} = \frac{C_{t_2} - C_{t_1}}{t_2 - t_1}
$$

In your example, if you divide both counter values $C_{t_1}$ and $C_{t_2}$ by a factor of 64 (for instance, by right-shifting by 6 bits), then the above equation becomes:

$$
f = \frac{1}{64}\frac{C_{t_2} - C_{t_1}}{t_2 - t_1}
$$

Therefore, the action of dividing the counter's output value by 64 has the likewise action of dividing its output frequency by the same factor.


How to deal with "weird" phase plots in bode diagram when designing a controller I am trying to design a balance controller for a robot. With MATLAB simulink I arrived at the transfer function between the input and the pitch angle for the robot. I have plotted the bode and Nyquist plot below.

As you can see, the Nyquist plot looks very weird, and in the bode diagram the signals at low frequencies are attenuated. To fix this, I added a zero at the peak of the magnitude plot, so
$$\omega_i=10 \text{rad}/s \Rightarrow \tau_i = \frac{1}{\omega_i}=0.1\text{s}$$
This results in an I-term which looks like this: $C_i(s)= \frac{\tau_is+1}{\tau_is}=\frac{0.1s+1}{0.1s}$
When I add this to the controller the bode plot now becomes.

Now I also want to add a lead-term to prevent overshoots but the way I usually do that is to change the crossover frequency $\omega_c$ such that there is good phase margin, define some constant $\alpha$, and then find the lead term with the formula $\tau_d = \frac{1}{\omega_c \sqrt{\alpha}}$. Then I can make my lead-term with
$$C_d = \frac{\tau_d s + 1}{\alpha \tau_d s + 1} $$
But when I look at the current phase plot, it all around looks like there is no trouble regarding the phase margin. The phase is never close to $-180^{\circ}$ so I am not sure how to proceed with such a weird looking phase plot.
So my question is, how do I find my lead-term when the phase plot looks like this? Is there another way to do than I have proposed?
Edit
The open loop transfer function between the input and pitch angle is $G(s) = \frac{ 5.893e05 s^4 + 7.71e08 s^3 + 1.435e11 s^2 + 2.566e12 s + 7.683e09}{s^7 + 2418 s^6 + 1.317e06 s^5 + 1.977e08 s^4 + 1.312e10 s^3 + 1.796e11 s^2 - 9.658e11 s - 1.525e13}$
Using the "pole" command in MATLAB I noticed that there was a pole in the RHP which should cause some instability.

I'm a bit puzzled with your Bode plot, here's what I get with your transfer function and Matlab 2019b. My phase starts at almost -180 degrees while yours start at 0 degree. Is it possible that you've made a mistake somewhere in your analysis?

Your transfer function is a bit complex, you could simplify it by removing the fast poles and zeroes. For example, you have 2 fast poles at -1714 and -516 and a fast zero at -1088. You can remove the 2 fast poles and the fast zero to simplify your transfer function. This will make your analysis easier.
Secondly, I think the rule-of-thumb formulae you use for the cross-over frequency is only valid when your open loop system has a low-pass characteristic. Your system has a band-pass characteristic.


Is it possible to have high energy compaction in audio signal? We all know that the energy compactness in image can be achieved in high scale.But in case of audio signal is it possible to have such high energy compactness?

Depends on what you mean exactly by "energy compaction".

Most audio signals are "pink" in nature, so doing any type of frequency transformation will "compact" the energy towards the lower frequencies and most audio codecs use this property. Not sure you get "high" energy compaction: a typical FLAC compression rates are about 30% size reduction or so.

It's more difficult to compact a signal in time through some sort of linear transformation or invertible transformation. This would really depend on what exactly you are trying to do and what the general properties of your signal are.


Significance of sequency ordering in Walsh-Hadamard matrices So I've been doing some reading into Hadamard and related matrices (Slant, Walsh, , and have observed that there are three typical "orderings" of the rows based on various criterion. E.g. "Natural", "Dyadic", "Walsh" (or Sequency http://fourier.eng.hmc.edu/e161/lectures/wht/node3.html).

One of the books I'm reading is "Introduction to orthogonal transforms with applications in data processing and analysis" by Ruye Wang (of which a copy can be accessed freely here http://fourier.eng.hmc.edu/book/lectures/mybook.pdf )

In this book there is a chapter (8) dedicated to some of these matrices. In particular is the Hadamard matrix. It states a definition of sequency:

  "..the sequency s of each row, defined as the number of
  zero-crossings or sign changes in the row. Similar to frequency,
  sequency also measures the rate of changes or variations in a signal.
  However, sequency can measure non-periodic signals as well as periodic
  ones"

Later on in the discussion of making the Hadamard matrix into a seqeuency ordered Walsh-Hadamard matrix it states:

  "The rows (or columns) of the WHT matrix H, ... , are not
  arranged in the order of the sequencies, while it is desireable to
  arrange them according to the sequencies in a low-to-high order,
  similar to how the DFT coefficients are arrange."

In the section on the Haar transform, it makes a brief explanation comparing the Fourier, Walsh-Hadamard, and Haar transforms:

  "It is interesting to compare the Haar transform with other orthogonal
  transforms such as the Fourier, cosine, Walsh-Hadamard, and slant
  transforms discussed before. What all of these transforms, as well as
  the Haar transform, have in common is that their coefficients
  represent some types of detail contained in the signal, in terms of
  different frequencies (Fourier transform and cosine transform),
  sequencies (Walsh-Hadamard transform), or scales (Haar transform), in
  the sense that more detailed information is represented by
  coefficients for higher frequencies, sequencies, or scales"

One can observe an analogy between the fourier and sequency ordering ( Aittokallio et al, 2001 - Testing  for  Periodicity  in  Signals: An Application  to  Detect  Partial  Upper Airway  Obstruction  During  Sleep:

  Walsh  functions with even  and  odd  orders  are  called the cal and
  sal functions,  respectively, and  they correspond to  the cosine and
  sine  functions in Fourier  analysis.  One should note, that in
  contrast to harmonic functions the

TLDR

Now, i've been seeing the Walsh-Hadamard matrix in various papers regarding locality sensitive hashing, and approximate randomization, random projections etc. Returning to the Walsh-Hadamard, my confusion is that I'm trying to form an intuition behind the significance of using a sequency ordered Hadamard matrix as opposed to just a standard Hadamard matrix, in general, and in instances when the WH matrix is combined with other matricies to induce randomization (e.g. see Subsampled Randomized Transform, e.g. Tropp 2010 - IMPROVED  ANALYSIS  OF  THE
SUBSAMPLED  RANDOMIZED  HADAMARD  TRANSFORM)

There isn't much significance. Except for few very particular applications like auto-focusing where you might try to adjust the lens such that the higher sequency components has maximum magnitude.  The FFT of a particular sequency isn't bandwith limited. Therefore sequency has very limited meaning.
You should indeed be looking at the applications like random projections, locality sensitive hashing, associative memory and neural networks:
WHT applications


FIR filter - remove padded zeros I am using this class to filter EKG data.

My sampling rate is 167/s and the length of data is 1670 (10 seconds), when I apply the filter I get the correct results but looks like the number of samples i get back is bigger (I guess because of zero padding to become power of 2).

How can I get back only the samples I filtered? I tried to copy first 1670 samples from the filtered array but looks like it contains zeros as well.

Thanks

Given a linear phase FIR filter $h[n]$ of length $L$, if you use the conv or filter function to produce the output samples $y[n]$ for a given input $x[n]$, then there will be a constant group delay $N$ associated with the type of FIR filter being used. So this delay of $N$ samples must be taken into account when evaluating the output result in applications where data phase is critical such as in the case of using a Hilbert transform to produce an analytic signal...

Assuming that you have a filter with odd length $L$ which can be expressed as
$$ L = 2 N + 1$$ then the associated group delay will be the integer $N$; i.e., the information in the output data will be shifted $N$ sample to the right, hence you may select the output after the Nth sample as such:
$$ y_i[n] = y_d[n-N] $$.

In a computer environment such as MATLAB/Octave this translates into matrix indices such that:

xi = yd(N+1:end-N)

where $y_d$ was the matrix that holds the result of the raw convolution output;i.e.,

yd = conv(x,h)

if you have used filter function to produce the raw data than the slicing
should be done as:

xi = yd(N+1:end)

where $yd$ now is the matrix that holds the result of the raw filter output;i.e.,

yd = filter(h,1,x)

When the filter length $L$ is an even number, there happens a half-sample delay of information associated with the raw-data output of the computation and you should consider that as well. Therefore It's better to use an odd length FIR filter unless you have convincing reasons to do the otherwise.


How to implement the $\tt conv(x,h)$ in MATLAB without using loop? I'm trying to calculate convolution of two given vectors in MATLAB without using loop, and of course without the function conv itself, but I can't remove the last loop I've used in the code below. Any hint would be appreciated.

function [ convolve ] = myconv( x,y)
%This function calculate the convolution of the two inputs
lx=length(x);
ly=length(y);
lt=lx+ly;
x=[x zeros(1,ly)];
y=[y zeros(1,lx)];
convolve=zeros;
for n=2:lt
    convolve(n-1)=sum(x(1:n-1).*fliplr(y(1:n-1)));
end

what about this:

function [convolve] = myconv(x, y)
    L = length(x) + length(y) - 1
    convolve = ifft(fft(x, L) .* fft(y, L));
end

Using the convolution theorem.

Regarding deconv function, check this out:

x = randn(10,1);
y = randn(10,1);

L = length(x) + length(y) - 1;

z = ifft(fft(x, L) .* fft(y, L));
z2 = conv(x, y);
sum(abs(z-z2))

x2 = ifft(fft(z, L) ./ fft(y, L));
x3 = deconv(z, y);
x0 = [x; zeros(L-length(x),1)];
sum(abs(x0-x2))
sum(abs(x3-x))

Program output:

ans =

   1.1616e-14

ans =

   5.4956e-15

ans =

   2.3924e-12

The first sum calculates the deviation between conv and my proposal. The second sum shows difference between the Fourier-Domain decovonvoled sequence and original sequence (note that the output of the FFT has L samples, where the last L-len(x) samples are 0). The third sum compares the output of deconv and the original sequence. You see that they all match up to floating point errors.


What is the correct solution for Fourier transform of unit step signal? The unit step signal defined as

$$
 u[n]= \lbrace 1; n>=0; \\
               \qquad0; n<0 \rbrace
$$

has three possible solutions for its Fourier domain representation depending on the type of approach. These are as follows -

The widely followed approach (Oppenheim Textbook)- calculating the Fourier transform of the unit step function from the Fourier transform of the signum function.

$$
     F(u[n])=U(j\omega) = \pi\delta(\omega)+\frac{1}{j\omega}
$$

Fourier Transform calculated from the Z transform of the unit step function (Refer Proakis Textbook, Digital Signal Processing Algorithms and applications,  pages 267,268 section 4.2.8 )

$$
     U(j\omega) = \frac {e^{\frac{j\omega}{2}}}{2j\sin \frac{\omega}{2}};   \omega \neq 2\pi k; k=0,1,2,3...
$$

Fourier transform calculated by splitting into even and odd functions - followed in Proakis Textbook (Refer Proakis Textbook, Digital Signal Processing Algorithms and applications, page 618 section 8.1 )
$$
      U(j\omega) = \pi\delta(\omega)+\frac{1}{1-e^{-j\omega}}
$$

The 2nd representation can be ignored since it is not a well-behaved function. But the approaches followed by Proakis and Oppenheim are equally valid (they extend the Fourier transform to include impulses in the frequency domain)
But the confusion is that they provide different solutions.

Is there any mistake in my understanding? or am I missing any crucial point?? Kindly help me understand this and the correct form that can be used in all applications. (I found that the Oppenheim approach is used in deriving the Kramers-Kronig Relations and the Proakis approach used in the derivation of the Hilbert transform)

As Matt said, the second and third definition are the same except for the part with impulse. The impulse ($\pi\delta(\omega)$) account for the DC value of $u[n]$. Without that term (i.e., the second definition) it is actually the FT of $v[n]=\frac{1}{2}\operatorname{sgn}[n]$. We have $u[n]=v[n]+\frac{1}{2}$. And hence the FT of $u[n]$ has the additional term to account for the addition of $\frac{1}{2}$. Also, the discrete time FT (or DTFT) of $u[n]$ is properly written as $U(e^{j\omega})$.

The first definition, $U(j\omega)$ is the "continuous time" FT (or CTFT) of $u(t)$ (not $u[n]$) and hence different from the other two definitions.


What is the relationship between a digital audio signal and sound pressure level? I'm playing a pure tone on an iPad through a pair of "standard" earphones.  From a programming perspective, the output level is in the range of $0$ to $1$ where $1$ means "full volume".

I've observed using a sound pressure level meter that when playing a pure tone at $1000 \text{Hz}$ with output level in $0$ to $1$ with step size $0.05$, that the relationship between output level and sound pressure level is logarithmic.

However, I'm not sure if this is correct or an issue with my use of a sound field sound pressure level meter (instead of a coupler).

Should doubling the output level double the sound pressure level?  That is, should the relationship be linear?  It makes intuitive sense to me that as the output level or voltage is increased, the earphone speaker will vibrate more, and thus the sound pressure will also increase.

As you say, pressure will be proportional to output level, if the output level is equivalent to voltage. But sound pressure level (SPL) is by convention expressed in the logarithmic decibel (dB) scale, with 0 dB referenced to some standard sound pressure level. The relationship is:

$$\text{SPL} = 20\log_{10}\left(\frac{A}{A_\text{ref}}\right) \text{ dB,}$$

with a reference output level $A_\text{ref}$ that would give an SPL reading of 0 dB. For measures of power such as $A^2$ the formula changes slightly:

$$\text{SPL} = 10\log_{10}\left(\frac{A^2}{A^2_\text{ref}}\right) \text{ dB.}$$

With $A$ and $A_\text{ref}$ non-negative, the equations are equivalent with the square inside the logarithm or outside it as a factor $2$.


Digital Audio Mixing Fundamental. What will be the data rate of mixed audio? I have voice data at 128 kbps coming from three different sources. I'm looking for mixing all three voices together and sending it back to speakers. What methodologies are there for mixing digital audio bit streams and what will be the mixed voice data rate? will it be 128*3 or just 128 only?

When you add multiple audio sources, you must first make sure that they use the same sampling rate, and then same numerical format, so that addition can be realized.
So, if they are not at the same sampling rate already, then you have to resample them to the same rate. And also make sure (by type conversions) that all of them to use the same numerical data format (same type and precision, etc.)
When you say 128 kbps alone; it looks like either an compressed data stream such as mp3, or PCM raw data at 16 kHz, 8-bits per sample. Assuming the latter is the case, then the resultant stream will again be an 16 kHz, 8-bits per sample PCM data stream, having a bit rate of the same 128 kbps.


Phase shift of all-pass bi-quadratic filter Does anybody know how to derive this phase function for that all-pass filter structure? Many papers re-print them, and corresponding chart as well.

But nobody mentions that at normalized frequency 0.25, tan(pi/2) goes to +- infinity, and atan jumps from + to - pi/2 :

Is phase shift periodic here by pi, so we can move the right curve down to glue them at 0.25. What does that mean in real filter operation?

First of all, that formula for the phase has a sign error, at least if $\alpha$ is defined as in the diagram of the filter structure. Second, your phase plot shows a jump from $-\pi$ to $\pi$, which is no jump at all because a phase of $\pi$ is the same as a phase of $-\pi$. The jump occurs only because the arctangent function returns the principal value of the phase in the interval $[-\pi,\pi]$. You can get rid of the jump by using the phase unwrapping function unwrap in Matlab/Octave. But note that this is just a cosmetic thing, there is no phase jump in the physical sense of the word.

Now for the derivation of the formula for the phase. From the given filter structure we get the following equation for the $\mathcal{Z}$-transforms of the input and output signals:

$$Y(z)=X(z)z^{-2}+\alpha\left(X(z)-Y(z)z^{-2}\right)\tag{1}$$

From this the transfer function can be computed as

$$H(z)=\frac{Y(z)}{X(z)}=\frac{\alpha+z^{-2}}{1+\alpha z^{-2}}=\frac{\alpha z+z^{-1}}{z+\alpha z^{-1}}=\frac{B(z)}{B\left(\frac{1}{z}\right)}\tag{2}$$

with $B(z)=\alpha z+z^{-1}$.
$H(z)$ as given by $(2)$ is indeed an allpass transfer function with $|H(z)|=1$ for $|z|=1$. Note that the system described by $(2)$ is only stable for $-1<\alpha <1$.

On the unit circle $z=e^{j\omega}$, $H(z)$ can be written as

$$H(e^{j\omega})=\frac{B(e^{j\omega})}{B^*(e^{j\omega})}\tag{3}$$

where $^*$ denotes complex conjugation. Consequently, the phase of $H(e^{j\omega})$ is given by

$$\phi_H(\omega)=2\phi_B(\omega)\tag{4}$$

where $\phi_B(\omega)$ is the phase of the numerator $B(e^{j\omega})$:

$$\begin{align}\phi_B(\omega)&=\arg\big\{\alpha e^{j\omega}+e^{-j\omega}\big\}\\&=\arg\big\{(1+\alpha)\cos(\omega)-j(1-\alpha)\sin(\omega)\big\}\\&=-\arctan\left(\frac{1-\alpha}{1+\alpha}\tan(\omega)\right)\end{align}\tag{5}$$

The phase of $H(e^{j\omega})$ is thus given by

$$\phi_H(\omega)=-2\arctan\left(\frac{1-\alpha}{1+\alpha}\tan(\omega)\right)\tag{6}$$

which differs from the given formula in the sign of $\alpha$.


Why does DPSK constellation shrink when noise is added? I am generating a DQPSK waveform at complex base band with constellation points on the unit circle and when I add AWG noise the constellation shrinks.  The noisy constellation "blob" migrate inside the unit circle.  My question is why is this happening and is it supposed to happen or am I missing something?

My method of adding noise is:

Xnoise = X + sigma/2 * (randn(N) + 1j*randn(N))

where sigma is the overall standard deviation for the noise.  It's being divided by 2 since I'm using two separate randns for the real and imag parts and the two Gaussians added together give me the correct variance.

It looks like your circle got bigger. The constellation is fine.


Extremely large reconstruction errors in NMF First time in this stack...
I am doing some audio analysis. I have a spectrogram ($N\approx 33000, M=1024$) and I need to run an NMF algorithm on it. I am using the Scikit learn implementation. And the code on Github.

Immediately, what strikes me strange, is if you look at the fit_transform functions, the returned signal representation is called W, which makes me believe that the dictionary learned is our H. Aren't the naming conventions for these two components the opposite?
Second, my reconstruction error is in the thousands. Any suggestions on what I may be doing wrong?

My spectrogram is non negative, as required.

When you decompose $X \approx WH$, the columns of $W$ form the dictionary and the columns of $H$ are the time varying weights. In your example, if you choose 150 components, then $W$ should be $33k \times 150$ and $H$ should be $150 \times 1024$.

The error "in thousands" could mean one of two things:

It's practically "zero" if you normalize it by the Frobenius norm of the spectrogram $X$, in which case all is well, or,
Perhaps your sound signal is too rich to be completely representable by 150 components and you should experiment by increasing that number.

It is also worth noting that squared error is often not a good measure for audio perception. Two signals with a large squared difference may sound almost the same.


Is the frequency spectrum dependant on frequency, or on imaginary angular frequency? The frequency spectrum of a time domain signal x(t) can either be written as X(f) or $X(j\omega)$. But how is the later correct? I mean, the frequency spectrum is clearly dependant on the frequency, not on a $j\omega$ part, is it?

Basicaly it's just a convention and the Fourier transform is a function of $\omega$ .

The imaginary argument $j\omega$ can be understood in different perspectives.
In one convention, for the class of signals $x(t)$ which are absolutely summable (i.e., stable) , the continuous time Fourier Transform (CTFT) $X(j\omega)$ is considered to be derived from the existing Laplace transform $X(s)$, with $s$ being a complex number, by replacing $s$ with $j\omega$, which amounts to evaluating the Laplace transform on the imaginary axis ( where $s =j\omega$ )  of the complex plane. Hence in this interpretation, the argument $j\omega$ retains the conceptual and mathematical link between the Laplace transform $X(s)$ and the CTFT $X(j\omega)$ where we say $X(j\omega) = X(s)|_{s=j\omega}$

A similar argumentation also holds for the discrete time Fourier transform (DTFT) $X(e^{j\omega})$ as such: For the class of signals $x[n]$ which are absolutely summable (i.e., stable), the discrete time Fourier Transform $X(e^{j\omega})$ is considered to be derived from the existing $Z$-transform $X(z)$, where $z$ is a complex number, by replacing $z$ with $e^{j\omega}$, which amounts to evaluating the $Z$-transform on the unit circle ( where $z = e^{j\omega}$)  of the complex plane. Hence in this interpretation, the argument $e^{j\omega}$ retains the conceptual and mathematical link between the $Z$-Transform $X(z)$ and the DTFT $X(e^{j\omega})$ where we say $X(e^{j\omega}) = X(z)|_{z=e^{j\omega}}$

However, note that, there exist certain very important class of signals (such as the pure sine wave or the impulse response of an ideal Low-Pass filter) which are not absolutely summable (unstable) therefore whose formal Fourier transforms (as obtained by replacing $s$ or $z$ by the corresponding $j\omega$ or $e^{j\omega}$) do not exist. Those Fourier transforms are obtained from another perspective by utilizing the impulse function $\delta(\omega)$ or discontinuities in the frequency domain. In such a case the argument is misleading and there exist no coresponding $X(s)$ or $X(z)$ which you can obtain by replacing $j\omega$ or $e^{j\omega}$ with $s$ or $z$ in the Fourier transform $X(j\omega)$ or $X(e^{j\omega})$, hence in that case it's just a convention.


Apply Low Pass (Smoothing) Filter on a Set of Matrices and Reject Outliers Given a set of a $ 3 $ by $ 3 $ matrices $ {H}_{i} $.
Each matrix is an Homography matrix.
They are used to stabilize Video Stream.
Yet some of them are outliers which creates "Jumps" in the "Stabilized Video Stream".
I would like to apply Low Pass filter on matrices in a manner which "Low Pass" the homography (Meaning, the Mapping).
Moreover, I'd like to reject Outliers, Matrices which compared to other rotates / scales the frame in a completely different way.
The problem I don't know any metrics which deals with that property of a $ 3 $ by $ 3 $ matrix.

Any idea?

The method used was a parametric representation of the matrix $ H $.

For intuition, thing of rotation matrix in $ \mathbb{R}^{2} $, namely $ H \in \mathbb{R}^{2 \times 2} $.
Since we have series of those (In case of video, a series in time) we can define $ H \left( t \right) $.

As a rotation matrix it can be parameterized by a single parameter $ \theta $ which is the rotation value. Since we are dealing with a time series then we have $ \theta \left( t \right) $.

Now, this is just a simple Time Series and regular methods of smoothing and outlier rejection can be applied.

Another idea would be to find the best base to approximate all matrices (Which with respect to that base the can be diagonalized) and then work on the "Singular / EigenValues".


Preserve specific information in compressed sensing I have a signal that isn't perfectly sparse and I would like to apply compressed sensing on it for lossy compression. However I would like to preserve specific section of the signal so that this "important" section can be perfectly reconstructed.
Is it possible? How?
(I hope this is ok to ask this kind of question - I'm not very deep in the theory of CS, just playing with code samples after reading brief overviews).

As a result of the discussion, let me provide an answer here.
Say you want to perform "linear compression" àla Compressed Sensing to a vector, i.e., reduce a length-$N$ vector $\mathbf x$ to a length-$M<N$ vector $\mathbf y$ via a linear transformation $\mathbf \Phi$ such that $\mathbf y = \mathbf \Phi \mathbf x$. Further assume that $\mathbf x$ may be partitioned into $\mathbf x = [\mathbf x_1^T, \mathbf x_2^T, \mathbf x_3^T]^T$, where $\mathbf x_i$ are length $N_i$ vectors and $N = N_1 + N_2 + N_3$. Further, $\mathbf x_2$ is a part of your signal you do not want to compress so that the actual compression is only applied to $\mathbf x_1$ and $\mathbf x_3$.
The most general way of achieving this is by structuring $\mathbf \Phi$ as $$\mathbf \Phi = \begin{bmatrix} \mathbf \Phi_{11} & \mathbf 0 & \mathbf\Phi_{13} \\
\mathbf 0 & \mathbf I_{N_2} &\mathbf 0 \\
\mathbf \Phi_{31} & \mathbf 0 & \mathbf\Phi_{33}
\end{bmatrix} $$ where $\mathbf \Phi_{ij}$ is $M_i \times N_j$ for $i, j \in \{1,3\}$ and $M = M_1 + N_2 + M_3$.
Applying this $\mathbf \Phi$ we can see that $\mathbf y$ will be partioned as well, i.e., $\mathbf y = [\mathbf y_1^T, \mathbf y_2^T, \mathbf y_3^T]^T$ where $\mathbf y_1 = \mathbf \Phi_{11} \mathbf x_1 + \mathbf \Phi_{13}\mathbf x_3$, $\mathbf y_3 = \mathbf \Phi_{31} \mathbf x_1 + \mathbf \Phi_{33}\mathbf x_3$ and $\mathbf y_2 = \mathbf x_2$, as desired.


pyo: Retrieve amplitudes of audio signal frequenciy ranges I'm playing around with audio signal processing for the first time and just to start somewhere, I decided to use the Python library pyo. My goal is to analyze an audio input signal from the microphone to detect different types of strikes on a drum and to emit an event (or call a callback, whatever) in that case. My first step is to define some frequency ranges, I want to observe and then to have a function called periodically and pass in the peak amplitudes within those defined frequency ranges.

I found that there is something like FFT/DFT and my naive idea is, that I can setup an FFT object with my input signal and create an instance of the Pattern class, that calls a function periodically. In that function I can then call the FFT object in a way, so that it will return me a list with 1024 items. The "whole" frequency spectrum is sliced into 1024 slices and each item in the list contains an amplitude value for that specific frequency range slice at the given time.

I tried the following code and now I don't know how to interpretate the values provided by the FFT.get() method. My question is: Are my expectations wrong or how do I use the FFT class to get my list?

from pyo import *

server = Server(nchnls=1, buffersize=128, duplex=1, winhost="asio").boot().start()

class FrequencyRangeAmplitudeProvider:

    __frequency_ranges = ()
    __listener = lambda x: x

    __fft = None
    __pattern = None

    def __init__(self, input, frequency_ranges, listener):
        self.__frequency_ranges = frequency_ranges
        self.__listener = listener

        self.__fft = FFT(input, size=1024, overlaps=4, wintype=0)
        self.__pattern = Pattern(function=self.__extract_and_emit, time=0.06).play()

    def __extract_and_emit(self):

        #print(self.__fft.get('real', True))
        #print(self.__fft.get('imag', True))
        #print(self.__fft.get('bin', True))

        my_list = []  # <-- How to get the data, I'd need?

        peak_amplitudes = self.__extract_peak_amplitude_in_frequency_ranges(my_list, self.__frequency_ranges)
        self.__listener(peak_amplitudes)

    @staticmethod
    def __extract_peak_amplitude_in_frequency_ranges(amplitudes, frequency_ranges):
        return (('A', -2.34), ('B', -.65))

frequency_ranges_to_observe = (('A', 430, 450), ('B', 870, 890))

def do_something_with_peak_amlitudes_in_frequency_ranges(amplitudes):
    print(amplitudes)

provider = FrequencyRangeAmplitudeProvider(
    Input(),
    frequency_ranges_to_observe,
    do_something_with_peak_amlitudes_in_frequency_ranges)

server.gui(locals())
```

pyo is not really the issue here at all. The question still stands if you were to substitute it with a million other packages of similar functionality.

What is at the centre of this question is "How to build a classifier to distinguish between different types of hits on a drum". Once you have that sorted, you can then implement it with whatever package you like, including pyo.

One of the simplest classifiers you can use for this task is the Nearest Neighbour Classifier.

The end effect of this classifier is not too different than your current thinking: I will split the frequencies of the sound and then try to work with if-then-else constructs, the final outcome of which would be the classification. Every if-then-else, puts a binary decision boundary on the outcome. For example, if you monitor the audio volume then you can create a rule that tells you whether there is something playing or not (if the volume is above a threshold, then something is playing). The nearest neighbour classifier over spectral features (of the FFT), does exactly the same thing but puts $N_{DFT}$ if-then-else on a given frame (it is not exactly if-then-else, but the end-effect is the same).

To get it to work, you would need a few representative recordings of each drum shot. For example, 20 different rimshots, 20 different straight hits. If you pass each one through the Discrete Fourier Transform (DFT) (of which the FFT is a quick implementation of) this will give you $c_{k,i,N_{DFT}}$ coefficients. That would be one set of $N_{DFT}$ coefficients for the $i^{th}$ recording of the $k^{th}$ type of drum hit.

Find the centroids (means of each DFT coefficient) that correspond to each drum hit and this will give you $k$ points in feature space that represent what each drum hit tends to sound like (obviously, here $k=2$).

To figure out what sort of drum hit has just been recorded (a drum hit you have not encountered during the "training" phase), take the waveform of it, pass it through DFT, take its DFT coefficients and work out which centroid they are closer to.

Doing this in real time might be a bit more challenging as you would have to establish an audio "loop" that continuously records frames from the audio interface, transforms them and passes them over to the classifier to make a decision.

Eventually, you will probably also need a way to describe a "third class" for sounds that are not drum hits and in which case you would like your classifier to remain "silent".

Overall, the system will also require some tuning to bring it as close as possible to a real time identification of drum hits, presumably, from a clean audio source, but it would be a start.

Hope this helps.


Using the traditional instruments for measuring high frequency signal All we know that we can not use the probe for measuring the high-frequency signal (microwave signals). I still can not understand the reason. I have read that the reason because of the wavelength size of the microwave signals is comparable to the electronic components, but this does not make any sense to me as the wavelength of the 1GHz is about 20 cm which is already comparable to the size of the electronic components.

Could any of you please explain to me in more details why we can not use the probe for measuring the high-frequency signal.

Not sure to understand what you are asking, Microwaves wavelengths range between 300 MHz (1 m) and 300 GHz (1 mm) which is indeed of the same order of magnitude as the probe.

1GHz corresponds to a wavelengths of about 30 cm  and is also in the Microwave range.

Are you saying that it is possible to measure a frequency of 1GHz with a probe while it is usually said that Microwaves can't be measured with it?


EmguCV: Detect a ball in the frame I am having a VERY hard time capturing when the ball comes into frame in this video

http://www.youtube.com/watch?v=x1V6_OHrXIo

I have no idea what to set the parameters too, and trying random values doesn't seem to help at all, I can't see the pattern to know what effect the numbers are having.  This is what I'm attempting.

Background subtraction, or color detection won't work, as they could easily change.  I need a more generalized solution.

var laneSixBedRect = new Rectangle(new Point(240, 185), new Size(240, 150));
var laneSixBed = imgSource.GetSubRect(laneSixBedRect);

TrackShot(laneSixBed);

private static void TrackShot(Image<Bgr, byte> laneBed)
{
    var laneBedGray = laneBed.Convert<Gray, byte>();

    var circles = laneBedGray.HoughCircles(
        new Gray(150),  //cannyThreshold
        new Gray(75),   //circleAccumulatorThreshold
        2.0,            //Resolution of the accumulator used to detect centers of the circles
        20.0,           //min distance
        5,              //min radius
        50              //max radius
        )[0];           //Get the circles from the first channel

    foreach (CircleF circle in circles)
        laneBed.Draw(circle, new Bgr(Color.Brown), 2);
}

Edit: This is my currentTrackShot method, it's running about 70% as far as accuracy, still some false positives, and some missed positives.

    private static void TrackShot(Rectangle laneRect, Image<Bgr, byte> prevImage, Image<Bgr, byte> currImage, Image<Bgr, byte> background)
    {
        var laneBedBack = background.GetSubRect(laneRect);
        var laneBedPrev = prevImage.GetSubRect(laneRect);
        var laneBedCurr = currImage.GetSubRect(laneRect);

        var laneBedBackGray = laneBedBack.Convert<Gray, byte>();
        //TODO: Should these be Smoothed?
        var laneBedPrevGray = laneBedPrev.Convert<Gray, byte>().SmoothGaussian(5);
        var laneBedCurrGray = laneBedCurr.Convert<Gray, byte>().SmoothGaussian(5);

        //TODO: Uncomment to see difference between frames, nearly black when no motion
        //var diff = laneBedPrevGray.AbsDiff(laneBedCurrGray);
        //diff.Convert<Bgr, byte>().CopyTo(laneBedCurr);

        //TODO: Improve these values
        var prevFeatures = laneBedPrevGray.GoodFeaturesToTrack(100, .01, .1, 3)[0];

        var returnFeatures = new PointF[1];
        byte[] status;
        float[] trackError;
        OpticalFlow.PyrLK(laneBedPrevGray, laneBedCurrGray, prevFeatures, new Size(15, 15), 5, new MCvTermCriteria(5), out returnFeatures, out status, out trackError);

        //TODO: Now how do I refine this???
        for (var i = 0; i < returnFeatures.Length; i++)
        {
            var prevPoint = prevFeatures[i];
            var currPoint = returnFeatures[i];
            var state = status[i];
            var error = trackError[i];

            var line = new LineSegment2D(new Point((int)prevPoint.X, (int)prevPoint.Y), new Point((int)currPoint.X, (int)currPoint.Y));
            if (state == 1 && line.Length > 15 && error < 10)
            {
                laneBedCurr.Draw(line, new Bgr(Color.Green), 3);
                Console.WriteLine("E: {0}", error);
                Console.WriteLine("L: {0}", line.Length);
            }
        }

        viewerModified.Image = currImage;
    }

EDIT2: My ulimate goal is to be able to draw a line along the trajectory the ball took and store that trajectory in some fashion that can be compared to others.

Remove things you don't want

Since the camera is static, you might want to use a background remover first. I found that the standard one provided with OpenCV works pretty well. I create it like this in the Android OpenCV SDK (you can play with the parameters) :

backgroundSubtractor = new BackgroundSubtractorMOG(3, 4, 0.8);

Then, apply it to each image in your video (no need to re-instantiate it unless you move the camera, it will learn from the sequence of images).

Computing the motion

First, using the Hough Transform for fast-moving shapes isn't very useful. The Hough transform requires you to be able to detect the shape's boundary with some accuracy, which obviously you can't, thanks to the motion blur.

I would suggest looking into algorithms for computing the optical flow. I don't know how close EmguCV is to the original OpenCV, but you could try calcOpticalFlowPyrLK (computes the flow for a set of points in the image, usually obtained using the goodFeaturesToTrack method) or calcOpticalFlowFarneback (computes a dense flow for the whole image, but more CPU-hungry).

Additionally, there have been some people working on CUDA/GLSL versions of Lucas-Kanade (the first one) if you really want to have a fast, dense flow algorithm, but I haven't seen a lot of useful implementations.

See the OpenCV documentation here : http://opencv.willowgarage.com/documentation/cpp/motion_analysis_and_object_tracking.html#calcOpticalFlowPyrLK


Incorporating Radial Distortion in Imaging Model I'm currently developing a multi-view stereo system but I'm confused as to how all the standard equations in Structure from Motion are affected in the presence of radial and decentering distortion. For example, the standard imaging model is that of

a) World to Camera Coordinates --> translate to camera center and rotate to match orientation of camera (R,t)

b) Perspective Projection --> x = fX/Z and y = fY/Z

c) Conversion to Pixel Coordinates --> scale by no. of pixels/m and add principal point
(b and c can be combined into matrix K)

As a result, we can build a projection matrix P = K[R|t]. Thus a point X in world coordinates can be represented in pixel coordinates as u = PX.

We can also create a fundamental matrix F between two images such that a pair of corresponding point (f on image I1 and f' on image I2) satisfy the epipolar constraint;

x'Fx = 0.

The matrix F depends upon the matrix P. But, my question is that if we incorporate radial and decentering distortion, the matrix P can no longer represent the imaging model. How do we get an epipolar "curve" representation ?

Moreover, for things like triagulation back to 3D coordinates from 2D point.. how do we do it ?

If you know the radial distortion parameters (e. g. by calibrating the cameras), then you should simply compensate for the distortion before you do structure from motion.  You can either undistort the images before you do anything else, or you can undistort the coordinates of the points that you are trying to match.


What's the difference between SISO and SIMO in using matched Filters Suppose we are transmitting a signal $x$ in $OFDM$ with $N$ subcarriers, After performing Parallel/Serial conversion and then N-points
 FFT, we suppose to have signal $s = IFFT^Hx$, Then adding the Cyclic prefix $P$ will create the transmitted signal $M$ such that:
$M = s +  P$. Before transmitting that signal, let's use Passband with matched filter (Or rrcpulse filter) whose sampling frequency is
$Fs$ and period $Ts = 1/Fs$, and carrier frequency $Fc$.

In SISO (Singl input singl output) system, it's clear that we should use matched filters with same parameters of that matched Filter
which was used in transmitter side.

My question, suppose that we are using SIMO (Singl input Multi output) system, four copies of signal should be recieved with
total length 4-times the transmitted signal.
My question, In that case, Should we use same matched Filter parameters for the received signal of SIMO system? How should we deal
with received SIMO signal? Or I think the period Ts should equal to $Ts/4$  ??  What's paramters should be changed?..

thank you

I don't know why we can't use RRC filter with OFDM ?!! Yes, of course, we do can. The only disadvantage for that is the complexity. however the gain will be in spectral efficiency improvement and bandwidth efficiency because of avoiding the needs of guard interval which is used in the regular OFDM. In Engineering, we should always ask, why and why not? why should we use it and why can't we use it?

So, Firstly, it's important to know that ISI in context of OFDM is different about ISI in the other types of modulations. As known,  the ISI happens when the symbol, in time domain, has a tail to the time of next symbol, but  that's not the same case for OFDM. In term of OFDM, that issue can be solves in two steps, because of the orthogonality, (and remember the use of FFT and iFFT), and second we insert the guard interval like CP which create circular convolution to eliminate the ISI. Be careful, in multi-path channels and large delay spread, that ISI will be exist in OFDM !! so we should present new solutions to deal with. "That exactly what the example reference,you provided above, tries to deal with that issue too".

But the question you will ask now, Why regular OFDM doesn't use RRC filter?  In reality, As I mentioned above, Regular OFDM is using Rectangular filter instead of RRC filter for minimal complexity, but the cost of that low complexity is a.) inserting CP which reduce the bandwidth efficiency b.) spectral efficiency reduction by sidelobes.

These two issues in regular OFDM, CP insertion and low spectral efficiency, can be solved by using RRC-OFDM, this RRC-OFDM doesn't use CP and the spectral efficiency is improved very well there. That technique was first presented in 2001, and then patented as US7065150 in 2006.

Coming back to the question, Should we use one matched filter in received SIMO signal or Multi filters. That depends on the simulation you are doing. you can use only one, but be careful how implement it. Or try to filter each copy of signal separately and then collect them back again and continue your work.

thnx