How to get info (mean,std) out of normplot in matlab? Lets say we have the following data and then we have a normplot:

load census;
h=normplot(cdate);
[mean std]=normfit(cdate);

I know i can get the mean(1890) and std(62.0484) of the data using normfit fxn, but thats not what i need.

As you can see in the plot, the straight red line is the ideal for the data.

i want to know the mean and std of the line.

I think we can do it from the slope of the line.. using handles of the figure probably with h(3).Xdata and Ydata.., but somehow i cant figure it out.

Can you please help me out.

Thanks

Just made a few modifications to the normplot program to give the reference data as output as well. xexp is the reference data and yperexp is the reference probability
(reference is the red line data). You can copy paste the whole code below and run it

function [h,xexp,yperexp] = nplot(x)

%   Copyright 1993-2010 The MathWorks, Inc.

if size(x,1)==1
x = x';
end
[n, m] = size(x);

[sx,originds] = sort(x);

minx  = min(sx(:));
maxx  = max(sx(:));
if isnan(minx) % Data all NaNs, setting arbitrary limits.

minx = 0;
maxx = 1;
end
range = maxx-minx;

if range>0
minxaxis  = minx-0.025*range;
maxxaxis  = maxx+0.025*range;
else
minxaxis  = minx - 1;
maxxaxis  = maxx + 1;
end

% Use the same Y vector if all columns have the same count.
if (~any(isnan(x(:))))
eprob = ((1:n)' - 0.5)./n;
else
nvec = sum(~isnan(x));
eprob = repmat((1:n)', 1, m);
eprob = (eprob-.5) ./ repmat(nvec, n, 1);
eprob(isnan(sx)) = NaN;
end

y  = norminv(eprob,0,1);

minyaxis  = norminv(0.25 ./n,0,1);
maxyaxis  = norminv((n-0.25) ./n,0,1);

p     = [0.001 0.003 0.01 0.02 0.05 0.10 0.25 0.5...
     0.75 0.90 0.95 0.98 0.99 0.997 0.999];

label = {'0.001','0.003', '0.01','0.02','0.05','0.10','0.25','0.50', ...
     '0.75','0.90','0.95','0.98','0.99','0.997', '0.999'};

tick  = norminv(p,0,1);

q1x = prctile(x,25);
q3x = prctile(x,75);
q1y = prctile(y,25);
q3y = prctile(y,75);
qx = [q1x; q3x];
qy = [q1y; q3y];

dx = q3x - q1x;
dy = q3y - q1y;
slope = dy./dx;
centerx = (q1x + q3x)/2;
centery = (q1y + q3y)/2;
maxx = max(x);
minx = min(x);
maxy = centery + slope.*(maxx - centerx);
miny = centery - slope.*(centerx - minx);
yinter = centery - slope.*(centerx);

mx = [minx; maxx];
my = [miny; maxy];

% Plot data and corresponding reference lines in the same color,
% following the default color order.  Plot reference line first,
% followed by the data, so that data will be on top of reference line.
newplot();
hrefends = line(mx,my,'LineStyle','-.','Marker','none');
hrefmid = line(qx,qy,'LineStyle','-','Marker','none');
hdat = line(sx,y,'LineStyle','none','Marker','+');
if m==1
    set(hdat,'Color','b');
    set([hrefends,hrefmid],'Color','r');
end
if nargout>0
    h = [hdat;hrefmid;hrefends];

% Compute position to display.
yper = normcdf(y,0,1);
yperexp = normcdf(polyval([slope,yinter],sx));
xexp = polyval([1/slope,-yinter/slope],y);

end

set(gca,'YTick',tick,'YTickLabel',label);
set(gca,'YLim',[minyaxis maxyaxis],'XLim',[minxaxis maxxaxis]);
xlabel(getString(message('stats:probplot:Data')));
ylabel(getString(message('stats:probplot:Probability')));
title(getString(message('stats:probplot:NormalProbPlot')));

grid on;
box off;

for i=1:m
    % Set custom data cursor on data.
    hB = hggetbehavior(hdat(i),'datacursor');
    set(hB,'UpdateFcn',{@normplotDatatipCallback,slope(i),yinter(i)});
    % Disable datacursor on reference lines.
    hB = hggetbehavior(hrefends(i),'datacursor');
    set(hB,'Enable',false);
    hB = hggetbehavior(hrefmid(i),'datacursor');
    set(hB,'Enable',false);
    if m>1
        setappdata(hdat(i),'group',i);
    end
    setappdata(hdat(i),'originds',originds(:,i));
end

function datatipTxt = normplotDatatipCallback(obj,evt,slope,yinter)

target = get(evt,'Target');
ind = get(evt,'DataIndex');
pos = get(evt,'Position');

x = pos(1);
y = pos(2);

% Compute position to display.
yper = normcdf(y,0,1);
yperexp = normcdf(polyval([slope,yinter],x));
xexp = polyval([1/slope,-yinter/slope],y);

% Get the original row number of the selected point.
originds = getappdata(target,'originds');
origind = originds(ind);

% Get the group number, which is set if more than one.
group = getappdata(target,'group');

% Generate text to display.
datatipTxt = {
sprintf('%s: %s',getString(message('stats:probplot:Data')),num2str(x)),...
sprintf('%s: %s',getString(message('stats:probplot:Probability')),num2str(yper)),...
''
};
datatipTxt{end+1} = sprintf('%s: %s',...
getString(message('stats:probplot:Observation')),num2str(origind));
if ~isempty(group)
datatipTxt{end+1} = sprintf('%s: %s',...
    getString(message('stats:probplot:Group')),num2str(group));
end

datatipTxt{end+1} = '';
datatipTxt{end+1} = sprintf('%s %s: %s',...
getString(message('stats:probplot:ReferenceLine')),...
getString(message('stats:probplot:Data')),...
num2str(xexp));
datatipTxt{end+1} = sprintf('%s %s: %s',...
getString(message('stats:probplot:ReferenceLine')),...
getString(message('stats:probplot:Probability')),...
num2str(yperexp));


does aliasing occur always if i sample a vibration in real world applications? I was reading about the aliasing effect and nyquist. I understand that aliasing effect occurs if the sampling rate is lower than twice the maximum frequency in the signal I want to sample.

So I was wondering if aliasing is always in signal measurements. Example: If I have a vibrations sensor that has a max sample rate of 8kHz -> It can reconstruct signals till 4kHZ perfectly right?
But what about frequencies which occur also in the measurements with much higher frequencies? If there is a vibration with 16 kHz and because I sample with much lower rate, there could be aliasing in my measurement and false signals I don't know about? is this possible? That means i could have peaks in my measurements which are because of aliasing?

That is why, before sampling, a (steep) lowpass filter with cutoff frequency $f_c \leq \frac{f_s}{2}$ shall be applied. Thus, the amount of aliasing will be insignificant.


Using MATLAB `xcorr()` to Find a Signal Hidden Behind Noise I am trying to develop some code to find a transient signal in a data set I have.  Before I get there however I am trying to write code using simulated data.

I first am creating a signal with a specified frequency and then adding white noise to it.

I then iterate through a loop where I create a reference signal at an incremental frequency.  I take the FFT of both signals and the I xcorr the results.

With this methodology I get descent results when I use a SnR ratio of greater than -15 dB, but after that my results are no good.  Is there anyway I can optimize this to find signals that are hidden deeper than -15 dB SnR?

My Code:

t = linspace(0,16*pi,1000);
y = cos(2*t);
z = awgn(y,-10);
fftz = fft(z);

subplot(2,1,1);
plot(t,z)
subplot(2,1,2);
plot(t,y)

freq = linspace(0.1,10,500);

for i=1:500
    reference_signal = cos(freq(i)*t);
    fftr = fft(reference_signal);
    Q(i) = max(abs(xcorr(fftz,fftr)));
end

subplot(3,1,1);
plot(t,y)
subplot(3,1,2);
plot(t,z)
subplot(3,1,3);
plot(freq,Q)

The output

You're basically doing a bank of hypothesis to find your signal using Matched Filter.

Though you use a slightly different method.

First of all, you should leave the signal in the time domain and calculate the cross correlation or their multiplication at the frequency domain.
Yet, since your signal doesn't have unknown phase (Or delay) multiplication will do (It covers the cross correlation at time zero) -> which is equivalent of what you do at the frequency domain.

Since you're using AWGN the performance depends on the Effective SNR.
The effective SNR is a function of the energy of the signal, variance of the noise and the number of samples.

If you set the signal and the variance of the noise, just use more samples (Of the Signal -> Equivalent of faster sampling rate).


Downsapled and blurred image I have an image with size $256*256$ then i get the Fourier transform of the image, If i only take the central region of the Fourier transform whose size is $64*64$ and get the inverse FT ,What is the effect on the output image$(64*64)$?

Down-sampling, Blurring, or both.

Have you fftshifted your results or not? The answers will be different:

If you fftshift your $256$ x $256$ 2-dimensional Fourier transform result, retain only the $64$ x $64 $ central co-efficients, inverse fftshift, and then inverse fourier transform the result, the effect will be a blurring of your original image.

This is because you have culled our your images' high frequency components, hence, the effect will be a blurring.

If however you have not fftshifted, then you are doing the opposite, and will actually be doing a high pass filtering. (Your image will certainly not be blurred).

By the way, hard thresholding fourier co-efficients like this is certainly used, but will give you ringing effects. However, those effects might be negligible and/or removable based on your application.


Intuitively, what is fourier series representation of a signal? Also intuitively what is frequency response? I know these formulae and can solve problems mechanically but i never got the core concept. Please help me out with this.

A fourier series essentially breaks apart a periodic signal to represent it as an infinite sum of sine waves that are in that signal. So a fourier transform of a signal breaks apart the signal to look at all the frequencies that are in that signal (essentially looking at all the ingredients that make up that signal, ingredients being the frequencies).

A frequency response is the frequencies of an output signal compared to the input signal. Essentially, a measurement of how well the output can replicate the input.


Best colour to detect I am using emgucv with opencv.

I have an object I wish to detect based on colour.
I convert the image to Hsl colour space and perform a filter on it.

The object I am detecting is perceived as yellow so I cover that range and make sure I am using/filtering good saturation.

This has been working well.

But I have now upgrade my usb camera and I have changed the lighting in the room.

Now running the same code I can it is ‘missing’ a lot of the object.

So, I inspect the hue of the pixels that have been missed out.

I can now see that the hues of those pixels have ‘shifted’ to the red range.  This incidentally is more noticeable when the object is closer to the camera.

So my question is this;

Is yellow an unreliable colour to detect compared to say red, green, blue etc?

Is the object not truly yellow?  Though when the object is further away it does not appear to be a problem,

Is colour segmentation unreliable?

Should I consider using a different colour space?

I will post my existing code if people require but I do not think it is required?

Thanks

A standard color camera, along with its software driver, is suppose to reflect the color that the eye would see. If it isn't doing that, then you have a badly designed camera or driver. In that case, either get a good camera+driver combo, or if the damage is linear you might be able to repair it by first multiplying the RGB with a 3x3 matrix that you have designed. Start with a unity matrix and adjust the elements until the output looks most like it does when you are looking right at the object (and assuming there is negligible color damage caused by the monitor you use to do this). Or you could even programmatically create the matrix by comparing the camera output with an ideal camera output. Once that is taken care of then you can pass the output to the algorithm that recognizes the object.

As far as the color space used for recognizing the object, a simple algorithm might do best keying off of the hue and saturation, so HSL might be best for that. But your real goal is to key off the intensities of all three colors (RGB) but still ignore luminance. So I suggest you use normalized RGB. That is, first scale the RGB so the total energy is some constant (say, unity). Specifically: ScaledRGB = RGB / VectorMagnitude(RGB). Once it's been normalized, you simply measure how far that normalized RGB "vector" is from the target color.

Other than the above points, there shouldn't be an issue with the particular color you are using, except that it should be most different from the particular background environment.

EDIT: answering subsequent questions...

If the scene lighting changes, you must first measure the change by noting the coloration of parts of the scene or the whole scene that always remain in view. For example, if you know of a mostly white object that remains in view and remains stationary, you can use that color vector as the normalization divisor rather than the magnitude of the pixel as described earlier (that is, divide a pixel's red by the white object's red, etc.). If you don't know of any white objects that remain in view and remains stationary, then take the average RGB of the entire scene (including the target object, if present, since you don't know where it is, yet. And also ideally compensate for known coloration of these background objects, if possible) and use that RGB as the divisor in your normalization (pixel's red divided by average red, etc.). This normalizes the scene for lighting and camera oddities. Then you can search for the object without regard of changes in lighting, etc..


What is the most computationally efficient way to generate 100 sine waves layered on top of each other? Right now I simply do this (pseudocode)

for freq in freqs:
    for i in len(samples):
        samples[i] += sin(...)

Is there a better way? This is for outputting to a 16-bit WAV file, using C on an ARM7 device. The frequencies are actually the centers of FFT bins for a particular bin size.

The most computationally efficient approach is to set up the data you want in the frequency domain by putting in the tones as "spikes" in the appropriate frequency bins, and then inverse FFT'ing it to get the time domain tones.  This approach will likely be a couple of orders of magnitude or so faster than your approach.

There are some downsides to doing it this way.  For tones whose frequencies are integer multiple of the sample frequency divided by the number of samples it is easy- just put a spike in at the appropriate frequency bin.  For all other frequencies, though, you need to put some energy into the neighboring bins as well.  I would have to think about how exactly to figure out how to do this.  I'm sure it can be done, but it's not a no-brainer.

Also, if you need to just produce a block of samples with the tones this approach is ideal.  If you need to create a continuous stream, though, care would have to be taken to get the phases right from block to block.  Again, not a no-brainer.

All that being said, I am pretty sure that this is the most computationally efficient approach.


Correct way of zero padding in time domain I am very new to signal processing and want to learn the correct way of zero padding for 'n' even and odd input signals.

For example,

$N=6\;\&\;M=32$

$x(n)=[a, b, c, d, e, f]\;\;\;n=0,1,\dots,N-1$

$x_{zeros}(m)=[a, b, c, d, e, f, zeros(M-N)]\;\;\;m=0,1,\dots,M-1$

-or-

$x_{zeros}(m)=[left\;zeros\dots a, b, c, d, e, f,\dots right\;zeros]\;\;\;m=0,1,\dots,M-1$

-or-

$x_{zeros}(m)=[a, b, c,\dots zeros\;in\;the\;middle\dots d, e, f]\;\;\;m=0,1,\dots,M-1$

What to consider when the input is $\mathrm{Real}$ or $\mathrm{Complex}$ of even and odd lengths?

Zero padding changes the bin size of the DFT resulting a finer frequency resolution. It is done by appending zeros to the end of the signal because we are only artificially increasing the length to decrease the bin size, $\frac{2\pi k}{N}$. An $N$ point DFT is just evaluating the DTFT on certain frequencies. If you were to append some zeros to the front of the signal then you are doing a time shift, and a shift in time domain is modulation in the frequency domain: $x[n-k]=X(\omega)e^{-j\omega k}$.


Get spectrogram matrix without specgram in MatLab I have a signal in Malab, and I need a matrix, probably transposed specgram matrix without using specgram tool. Is it possible?

Yes, the easiest way to do this is using reshape and FFT, since reshape will give you a matrix and the FFT will operate along a dimension of the matrix. I've applied a 1024 rectangular window on a data set of 65536 points, sampling frequency Fs (1024), and signal frequency (100).

x = randn(1, 65536) + cos(2*pi*100/1024*(0:65536-1));
z = reshape(x, 1024, []);
y = 20*log10(abs(fft(z)));
imagesc(y)

This example doesn't have any overlap, but all you would need to do is shift z then concatenate it with z prior to taking the FFT. This all depends on your overlap percentage.

z = [z; circhift(z, [0 -1])];

This is not complete given the first column is now at the end, so you would have to omit that column, but you get the idea.


Disagreement between Analytic Convolution and FFT Convolution I've been working on a small project, and as a part of this I need a function that can compute the structure function of the input data. In the process of writing this part/testing it I found the results were actually a poor match to what I would expect based on an analytic solution. Upon some further investigation, I've found this mismatch applies even to 2D convolutions and I have been completely stuck trying to debug this.

The test case that I am running is two rectangular functions convolved with one another, my expected result is a 2D triangle function (a pyramid essentially). But the difference between the analytic output and the FFT convolution output is large enough that I am concerned something is wrong. Furthermore, when I compute a direct convolution of the two inputs using scipy.convolve I get yet another output, different from the other two in it's own unique way!

To demonstrate the differences I'm seeing, I've plotted: my FFT Convolution's output, the analytic solution, the difference between them, and the difference between scipy.convolve's output and the analytic solution. My question then boils down to a couple things:

How large of a difference is expected between the FFT convolution and the analytic solution?

I would have expected something small (of order 1e-15 or less) but the peak difference I am seeing is 0.12. It also follows a very particular pattern, the vast majority of the output is close enough (difference of ~1.0e-16).

Is there a key step or concept that I am missing here?

I've been taking the FFT of each input, multiplying them together, and then taking the inverse FFT (with some scaling based on pixel size).

The code that produces the figure can be found in this relatively simple jupyter notebook on my github.

After playing with the notebook for a while, I figured out the following:

The result from scipy is shifted by one pixel in each axis (probably due to its definition of the "same" mode). Or, alternatively, both the analytic result and pyfftw results are shifted the opposite way.

After working around that, the result is exactly (within error) the same as the pyfftw result. I did this as follows (though this drops some data on the edges):

B_padded = np.pad(B, ((0, 1), (0, 1)), 'edge')[1:, 1:]
C1 = scipy.signal.convolve(A, B_padded, 'same', 'direct') * delta**2

The analytic solution (which peaks at 4.0) is not quite right for your input data (31x31 pixels at 1.0, edges at 0.5, and corners at 0.25). The peak value (where the image and the window align perfectly) should be:

box_bulk = 1.0**2 * (31 * 31)
box_edges = 0.5**2 * 31 * 4
box_corners = 0.25**2 * 4
peak = (box_bulk + box_edges + box_corners) / 1024 * 4  # ~= 3.876

which is the answer produced by both scipy and pyfftw.

Gist (also includes plotting tweaks).

(disclaimer: friend of asker)


FFT bin interpolation for low frequency with DC offset I have a noisy signal in I/Q format with a peak frequency that constantly changes between 0.5Hz and 4Hz. The sample rate is 64Hz. I want an output showing the current peak frequency with a maximum lag of about 1 second.

This is my first attempt to solve this problem with FFT, and the question is: What is required to make the frequency estimate accuracy around 0.01 Hz?

The idea is to use the FFT to both do the filtering of unwanted higher frequency noise and do the estimation of the 0.5-4Hz signal peak frequency.

With one second, we have 64 complex samples go into the FFT and 64 complex magnitudes coming out. They will then represent frequencies from -32Hz to +31Hz. We look at the 6 bins representing 0Hz to 6Hz and want to interpolate here to find the peak.

The algorithm must run on a low power CPU, so computational power is quite limited and power consumption is an important parameter.

The algorithm looks like this in Scilab:

sample_rate = 64;
secs = 1;
nr = sample_rate * secs;

f_signal = 0.45; // Hz - The signals peak frequency

// Construct the ideal complex signal vector
BT = complex( cos(f_signal*2*%pi*((1:nr)-1)/nr) , sin(f_signal*2*%pi*((1:nr)-1)/nr) );

// Remove any DC offset (or not?)
// BTavg = sum(BT) / nr;
// BT = BT - BTavg;

// Make a time domain plot
figure(1);
clf();
plot(real(BT), 'red');
plot(imag(BT), 'blue');

// Perform FFT (and shift bins around to have -32Hz to +31Hz)
BF = fftshift(fft(BT));

// Create the associated frequency vector f(1) = frequency of bin 1
f = sample_rate * (-(nr)/2:(nr-2)/2)/nr;

// Make a frequency domain plot
figure(2);
clf();
plot(f(1:nr), abs(BF(1:nr)), 'blue')

// Find max peak in the first 6 bins
pmax = 0;
imax = 0;
for i=(33:38)
    p = abs(BF(i));
    printf("i = %d f = %f p = %f \n", i, f(i), p);
    if p > pmax then
        pmax = p;
        imax = i;
    end
end

// Interpolate to find better estimate of peak frequency
y1 = abs(BF(imax-1));
y2 = abs(BF(imax));
y3 = abs(BF(imax+1));
// d = (y3 - y1) / (2 * (2 * y2 - y1 - y3)); // Quadradic method
d = (y3 - y1) / (y1 + y2 + y3); // Barycentric method
fmax = f(imax) + d * (f(imax+1) - f(imax));

printf("Estimated = %.2f Hz Actual = %.2f Hz Error = %.2f Hz\n", fmax, f_signal, abs(fmax-f_signal));

The input plot looks like this:

The output from the FFT (shifted) looks like this:

The console output showing the power from the 6 bins as well as the final estimated frequency and error is here:

  i = 33 f = 0.000000 p = 44.717018
  i = 34 f = 1.000000 p = 36.588120
  i = 35 f = 2.000000 p = 12.993838
  i = 36 f = 3.000000 p = 7.911240
  i = 37 f = 4.000000 p = 5.696681
  i = 38 f = 5.000000 p = 4.459184
  Estimated = 0.24 Hz Actual = 0.45 Hz Error = 0.21 Hz

As can be seen the error is way bigger than desired for this specific test frequency of 0.45 Hz. By changing the frequency it's clear that the error near zero only for near integer frequencies.

A caveat is that the input signal also has a varying DC offset, so ideally the DC offset should be "filtered" as well.

I am sure there are ways of saving power by only calculating the 6 bins we are interested in, but that optimization only makes sense if we can get the required accuracy from the algorithm in the first place.

So again: What is required to make the error on the estimates more like 0.01 Hz?

FFTs don't work well for sinusoidal frequency estimation if there is a nearby strong interfering signal.  Normally one needs at least 2 FFT result bins of separation to just start to resolve 2 frequency peaks.  For frequencies around or lower than FFT result bin +-2, any significant  DC offset will constitute that strong nearby interference (in bin 0) preventing resolving either frequency as an independent peak.

My guess is that, if the S/N ratio is low enough, a parametric estimator will produce a more accurate estimation, given that you have more samples than unknowns (3 sinusoidal parameters + 1 DC offset parameter).

Another possibility is to try to remove the DC offset by averaging or low pass filtering over a much longer window of samples, perhaps several hundred to thousands.

Added:
Strictly real valued sinusoids very near DC (or Fs/2) will also see interference from their own conjugate images in the spectrum.  Any input gain imbalance in the 2 quadrature components for a complex sinusoidal signal will contain some of that same conjugate interference as well.


Sampling: question regarding sampling period what exactly happens when I increase the sampling period ?

Providing some Matlab plots so that I can be more clear.

Using the signal $x(t) = \sin(10\pi t)$ and periods $T_s=0.02 ,\; 0.05 ,\; 0.1$

What happens exactly to the signal when the $T_s = 0.1s$ ?

You sample at times $nT_s$, so the values of the sampled signal are

$$x_d[n]=x(nT_s)=\sin(10\pi nT_s)\big |_{T_s=0.1}=\sin(n\pi)=?$$

Can you take it from here?


Difference between diversity gain and coding gain I am reading about space-time codes. They break down further into 2 sub categories called space-time trellis codes and space-time block codes. It says that trellis codes offer both diversity and coding gain. However, block codes only offer diversity gain. I don't quite understand this.

I understand the general forms of diversity gains such as time, frequency, site, etc. However, how exactly is space-time codes different from space-diversity or even MIMO?

Further more, what exactly is the criteria for coding gain? My intuition is that if we are actually decoding the original message which was encoded using parity bits, etc, then this is coding gain. This confusion carries over when I am trying to understand the distinctions between trellis and block codes.

MIMO refers to any system with multiple antennas, where signals are transmitted over the same bandwidth and at the same time.

Space-time coding (block or trellis, or a different one) is a specific way to transmit the data symbols. In space-time block coding, symbols are repeated over time and space in a specific pattern, which creates (potentially) both transmit and space diversity. They don't offer an intrinsic coding gain because the information is not encoded.

Trellis codes, on the other hand, allow for a coding gain because not all symbol sequences are permissible. This is similar to non-MIMO trellis modulation, where previous transmitted symbols define a sub-set of the constellation that can be used for the following symbol. This is a form of encoding, with the main benefit of neither increasing the data rate nor the required bandwidth.


What are the Even and Odd Components of a 2D Signal A signal $x$ can be decomposed into
its even, $x_e$, and odd, $x_o$, parts.
By,

$$
x_e(t) = \frac{1}{2} [x(t) + x(-t)]
$$

$$
x_o(t) = \frac{1}{2} [x(t) - x(-t)]
$$

What would be the equivalent for a multidimensional signal, such as a image signal?

The even part of a 2D image $x[n,m]$ is defined as:

$$x_e[n,m] = \frac{x[n,m] + x[-n,-m] }{2} $$

and the odd part is
$$x_o[n,m] = \frac{x[n,m] - x[-n,-m] }{2} $$.


Finding Argmax in MATLAB I have a Gaussian random variable $z$, whose mean $\mu$ and variance $\sigma^2$ are functions of two variables $x$ and $y$. Is there any way in MATLAB so that I can find the values of $x$ and $y$ for which $P(z)$ = $N(z$;$\mu$$(x,y)$,$\sigma$$(x,y))$ is maximised for a given value of $z$? Thanks in advance.

Try this:

[argvalue, argmin] = min(x);
[argvalue, argmax] = max(x);

Also, have a look to this question.


What exactly is a connected componnent in an image? The following is a paragraph from "Digital image processing, 4th edition, Gonzalez and Woods".

Let S represent a subset of pixels in an image. Two pixels p and q are said to be
connected in S if there exists a path between them consisting entirely of pixels in S.
For any pixel p in S, the set of pixels that are connected to it in S is called a connected
component of S. If it only has one component, and that component is connected,
then S is called a connected set.

I do not understand the last sentence: How can S have one component and that component is not connected??
Another question: In a binary image, can S be a connected component and have both 0 and 1 pixels?
By the way, the concept of component is not define precisely in the book.
Update:
In the book, a region is defined as a connected set. The following figure is an example of two adjacent regions (So, $R_i$ and $R_j$ are both connected sets.)

Thanks.

Let S represent a subset of pixels in an image.

This is typically a set of pixels with the same value, for example a group of pixels with the value 1, surrounded by pixels with a value 0.

Two pixels p and q are said to be connected in S if there exists a path between them consisting entirely of pixels in S.

A path is a set of steps that take you from pixel p to pixel q, and where each step is a jump from a pixel to a neighboring pixel. Thus, two pixels are connected if you can travel from neighbor to neighbor starting at p until you reach q.

If it only has one component, and that component is connected, then S is called a connected set.

If all pixels in the set S are connected to all other pixels through a path, then S is a connected component.
Thus a connected component is the very intuitive concept of a continuous region of equal-valued pixels.


Unsupervised learning algorithms to detect anomaly in waves I have a sample of graphs (more than 10000...).
that look like in the image below:

I am searching for an unsupervised learning algorithms that can help me to detect anomalous observations.

Here what I suggest for beginning: for every observation I have a collection of points $(x,y)$. With this collection, I find Fourier series with regression (I compute coefficients with the base $\{1,\sin(x),\cos(x),\sin(2x),\cos(2x)\dots\}$).
Now I have a set of coefficients instead of waves.

Somebody have an idea how to detect anomaly?

Try Nonparametric Statistics.
You may divide your random trajectories into parts and compare their statistic characteristics (mean and variance). Smth like this (but not exactly):

http://www.cscjournals.org/manuscript/Journals/SPIJ/volume7/Issue1/SPIJ-210.pdf


Licence plate detection with different backgrounds I just studied the Algorithm of ANPR (Automatic Number Plate Recognition) from the chp 5 of book Packtpub.Mastering.OpenCV.with.Practical.Computer.Vision.Projects in which they detect the vehicle licence plate , but their are some restriction in that algorithm , i call it restriction because its specific for one country and can not be used in my country , well in that algorithm they use the following algorithms to detect the licence plate area

Sobel filter
Threshold operation (otsu)
Close Morphological operation
Mask of one filled area
Svm classifier

But the thing they used is the background color of licence plate as their keypoint for the recognition and detection , but in my area licence plate have different backgrounds , most common are black , white and green , green is not much common but rarely used , black background used white number and white contain black numbers , they use floodfill algorithm for white background , but what should i use for 2 different backgrounds ? Size of the licence plate and size of character matter in any case ?

Here are the sample images of my regional vehicle's liecence plate

Here are a few ideas:

First convert from color to grayscale. It looks like you have fairly good contrast already. There are various methods to perform this conversion; choose the simplest at first: gray = (red + green + blue)/3. Quite often you don't need anything better than that. For some applications, using just the green color plane is sufficient.
If you want the algorithm to operate on black text, but if some images have white or light text on a dark background, reverse the polarity of the grayscale image by swapping grayscale 255 for 0, 254 for 1, and so on until all images have black text on a light-color background. You'll need an algorithm to detect when this is necessary, OR you can run an algorithm on the original grayscale image and on the negative (reversed brightness) image and choose the better result. Keep in mind that the "close" morphological operation has a different meaning for light-on-dark text versus dark-on-light text.
Investigate some adaptive thresholding methods. Using a global Otsu threshold can be problematic--if possible, you should first have some sort of method to detect regions that are likely to be license plates, and then threshold and binarize those regions only. Thresholding an entire image will likely cause trouble.
Look into the Stroke Width Transform and improvements upon it. The SWT or equivalent can help either with segmentation, or with identifying likely license plate regions, or possibly even with providing edges useful for identification.
Since in your region there are only a few color combinations, you might be able to create some method that would could a template for every possible color combination and then try to find a match at each of several image scales. The license plate might be 400 x 250 pixels, or 100 x 65 pixels, or whatever.
For multiscale operations, consider holding an image pyramid in memory. In some cases you can perform an initial match in the smallest image of the pyramid and work your way down. With some modification you could create an image pyramid has a scaling factor of less than 2 from level to level.
Pull the sample images into Photoshop or GIMP and play around with the filters. You may be surprised to find a filter that is helpful.
Whenever possible, work on raw, uncompressed images rather than on JPEGs or other compressed images.
You may not have much of a choice in hardware, but given the fuzzy look of your "LEA 1956" sample, it'd be nice to have more resolution (more pixels per image) and/or sharper optics and/or a raw image instead of a JPEG.

Try to generate the best possible images that you can with pre-processing before your SVM goes to work. Thresholding is often the weak link here. Thresholding and binarization is a common approach, but when you see how noisy a binary image can be you should think very hard whether that is the kind of image you actually want to process (no matter how many people may claim that this is the standard method and hence the best).

Even if you're already starting with a framework and/or algorithms, I recommend doing a little background reading to gain an understanding of both human reading and the many approaches to OCR:

Character Recognition Systems by Cheriet and others is a good reference (though a few years old now) of the various techniques for Optical Character Recognition.

Reading in the Brain by Stanislas Dehaene gives a good high-level overview of human reading, and covers the range of scripts. It's a relatively fast read, and has references to academic works that are useful.


Why does OFDM use cyclic-prefix while QPSK uses root-raised cosine filters? My understanding is that BPSK/QPSK/... are often implemented with root-raised cosine (RRC) filters to reduce ISI and bandwidth.

OFDM can be understood to be BPSK/QPSK/... implemented on multiple orthogonal carriers, with some additional complications like cyclic prefixes, pilot tones, and implementation using FFT/IFFT.

From what I've read, it seems OFDM doesn't usually use RRC filters.

Why is it that BPSK/QPSK/... use RCC while OFDM doesn't use RRC and instead uses a cyclic-prefix to combat ISI?  In other words, why isn't OFDM implemented the same way as BPSK/QPSK/... just with multiple orthogonal carriers?

Please let me know what I'm missing! :)  Thanks.

ISI, or intersymbol interference, means different things in the context of PSK and OFDM signals.  In PSK signals the symbols almost always have tails that extend, in the time-domain, into the times of other symbols.  This is what they mean by "intersymbol interference".

Unfortunately they have to do this to reduce the bandwidth of the signal.  They eliminate the negative effects of intersymbol interference by very carefully shaping the symbols so that when the receiver filters them through a matched filter, the other symbols have zero output at the sampling times of the other symbols.

With OFDM signals there is no overlap in the time-domain between symbols, so there is no intersymbol interference in the PSK sense of the term.  What they mean by ISI in the context of ISI is when time-delayed multi-path signals cause the time-delayed symbols to overlap with the following symbol in the non-delayed signal.  The cyclic prefix is added to the symbols to give the receiver enough un-corrupted data to recover the signal.

For instance, let's say that there are 64 sub-carriers and when we FFT them we get 64 samples of time-domain data.  Let's also say that we think the worst multi-path we will generally see will be a time delay of 7 samples.  If we make the cyclic-prefix 8 samples long then we should always have at least 65 "clean" samples, which gives us enough to recover the symbol.


What does it take to design a 4th order band-pass filter with -24dB/Octave? I would like to implement in software the same isolators that the Pioneer EFX-500 features:

The manual says the following about the isolators:

High-performance 3-band Isolator Function
[Fourth-order filter (–24 dB/ oct.)]
Level control is possible between +6 dB to –inf. and a separation of up to –24 dB/oct. for each band ensure accurate sound reproduction.

The thing is, I can program in C# or C++ but I have no experience at all in audio processing so I am a bit at a loss on how to get started.

I have found the winfilter software which can generate C code I can reuse but I am not sure about what settings to choose:

Based on the specifications I've put above, can you suggest which settings should be used ?

After experimenting with winfilter, generating the following filters mimics the EFX isolators, amazingly :)

2 * Butterworth LP 2nd order @ 300Hz
2 * Butterworth BP 2nd order @ 300Hz/6000Hz
2 * Butterworth HP 2nd order @ 6000Hz

So each signal passes through two filters of the same type, for some reasons I don't know the cause of, a 4th order BP in winfilter simply doesn't work and produces garbage instead.

As a reference, here's what really helped me out in matching the filters properties, white noise recorded for each of the bands then guesses until it matches:

Surprisingly, I'd say the above filters match the EFX about 99%, the only differences are at the edges where the winfilter ones really trim down the frequencies much more.


Analysis Received Signal at receiver side I'm trying to understand something which is really confusing me and I made my best to understand it but still confused.
If I have transmitted and receiver. Assume I want to transmit S(t) data and there's two multipath ways L1 , L2 in order to arrive the receiver.
Why when I have constructive interference then it's good for me ? in the paths/multipath there's might be destructive interference or constructive interference...my question why constructive interference is good for me..assume I get 2S(t) -constructive interference- in the received signal on the receiver ..why it's good for me ? what does it mean that I have 2S(t) and not S(t) ?
Im confused on the term constructive interference and how it's related to my received signal ..?
My confusion exactly is that-
I received 2X although I transmitted X , so 2X isn't X ..and we say that we received a good signal? if we received 2X which isn't equal to transmitted Data X so it's not the same data so it's considered distortion and it's not a good received signal.
thanks alot

Constructive interference just increases the gain of the signal at the receiver, which is designed (using automatic gain control) to adjust the signal level over a very wide range to that desired for reception. So the issue the OP refers to is of no real consequence as long as the received signal is in the (relatively wide) range of the receiver between maximum and minimum signal levels allowed.
Multipath distortion can come in different forms based on how fast the multipath is changing (such as in a mobile environment) resulting in "fast" or "slow" fading, and the relative distance between the paths compared the symbol rate, resulting in "flat" or "frequency selective" fading.  These are the real issues with multipath distortion which are corrected with equalizers in the receiver.
Please see this post which expands on these different fading models:
Rayleigh fading with frequency selective fading channel


time for a signal after a sum I have a signal $s(t)$ of duration $T$ with limited bandwidth. In matlab, I use a sampling frequency $fs$ and I have $L=T*fs$. After, I have $r=ones(size(s))+s$. To have a limited bandwidth, I filter with a raised cosine filter: $rr=filter(h,1,r)$ with $h$ the raised cosine filter.

My question: What is the duration for the signal $rr$? Personnally I would have said $T$. How to prove that?

I have a signal s(t) of duration T with limited bandwidth.

No, you don't.

A signal cannot be finite in both time and frequency at the same time. So the whole question about length in both domains comes down to how exactly you define it. Typically you need a criteria like "less than -100dB", "99.99% of the entire energy" or something like this.

These criteria are specific to your application. Maybe -40dB is plenty, maybe you need more than -200dB. Once you have a criteria defined you can apply it to your signals and determine the "lengths", but the result will always be tied to an application specific definition.


How to filter very low amplitude staircase signal from a background noise I have a staircase signal of a specific pattern that is drowned in a background noise. The ratio of signal to noise amplitude is nearly $10^{-4}$. I have attached the image of both with this post. I am an absolute noob in this field so please suggest me a direction on which path I should start to hunt for the solution.

I have tried to do fft on both of the signal and looked for any distinct signature to filter the background noise. Also I have tried with matched filtering and basic convolution techniques. But none of these seem to work. Your advice will be much appreciated. Thanks.

Linear filtering won't work here. If anything works at all then it might be something like total variation denoising, which is a non-linear technique for removing noise while preserving edges. This article is a very good read, and it also shows a Matlab implementation of the algorithm.


redundancy of sin and cos waves with real data I have the following question. Isn't it true that when applying a fourier transform to a real function (i.e. computing a characteristic function for a density), we only ever need one of the two waves: sin or cosin, to capture it's behavior.

I did a numerical exercise of moving back and forth between a probability density and a characteristic function, and every time it seems like either of the two waves is enough. Here are the waves:
Here are the original real data histogram (yellow), estimated smooth density (red) and the same density reconstructed separately from the sin (green) and cosine (blue) waves of it's own characteristic function.

Maybe for the complex input, the two waves are really necessary, but there is clearly some redundancy when the input is real.

Has anyone seen this before?

Both sin and cos waves are needed. Sure, you can find a signal (somewhere between a sine and a cosine, which seems as if it can be evaluated using either a sine or a cosine but in the general case you need both waves. And FT needs both waves anyway because internally it works with complex numbers regardless of the input  you pass to it (real input is usually interpreted as real0, img0, real1, img1...)

Suppose you have a pure sine wave and want to estimate its spectral content. Which wave would you use, a sine or a cosine?

It should be evident that only sine waves can accurately estimate the frequency and amplitude of your tested sine signal (provided that they're in phase and the frequency resolution is narrow enough to reduce spectral leakage effetcs). If you now inspect the DFT/FFT results, you'll see that only a single, imaginary coefficient (corresponding to the sine wave being tested) contributes to the overall amplitude/magnitude.  Likewise, cosine waves are necessary to accurately estimate the cosine components of your signal. Since most real signals are neither sin or cos waves, you need to use both to estimate them properly.

In addition to that, the DC component can only be represented by using a cosine wave (since sin(DC) = 0)

I personally don't think that you need to fully understand the math behind FT to realize that both waves are necessary.


Is this a valid z-transform of $ \frac{1}{n^2}u[n] $? I was looking over some DSP and came across the following signal: $x[n] = 1/n$. So I wondered whether it had a z-transform but I soon realized that it fails to meet two conditions described here:

  ...(1) the signal starts at a finite time and (2) it is asymptotically
  exponentially bounded...

So I modified it to meet those requirements: $x[n] = \frac{1}{n^2} u[n] $

(I assume $n>0$ is a condition or perhaps $u[n]$ should be changed to $u[n-1]$.)

Considering that signal, I believe it satisfies both requirements however I don't see anything like it in a table of common z-transforms, but it can be rewritten as:

$ x[n] = \frac{1}{1^2}*\delta[n-1] + \frac{1}{2^2}*\delta[n-2] + \frac{1}{3^2}*\delta[n-3]\ +\  ...$

Which I believe transforms to $ X(z) = \frac{1}{z^1} + \frac{1}{(2z)^2} + \frac{1}{3^2z^3} +\ ...  $

Does this transformation make sense? Is it the best approach? And if it's correct, is there a more compact way to express it?

What you've got is definitely the z-transform of that signal (by definition).  For a more compact expression, notice that the polylogarithm is exactly the infinite sum you've shown (just replace $ z \rightarrow 1/z$ in the Mathworld expression). Unfortunately I don't think there's any nicer way to write it; sometimes there just isn't.


Group Delay of a system Say I have an LTI system whose impulse response h[n] is:

$$
    h[n] = \delta[n] + \delta[n-4]
$$

I want to find the group delay of this system. I know the group delay ($grd$) is defined as:
$$
  grd = \frac{-d}{dw}(\angle H(e^{jw}))
$$

So I get $\angle H(e^{jw}) = -4w$ and the $grd = 4$ but I believe this is incorrect.

For a more complicated example:

$$
h[n] = -\delta[n+1] + \delta[n] + 2*\delta[n-1] + 2*\delta[n-2] + \delta[n-3] - \delta[n-4]
$$

What would my $\angle H(e^{jw})$ be? Thanks for helping me understand something fundamental that I am missing about group delay.

The z domain transfer function  of the system is the z transform of the system impulse response, so start by taking the Z transform of h[n] ...

$$H[z] = -z^1 + 1 + 2z^{-1} + 2z^{-2}  + z^{-3} -z^{-4}$$

You may be able to message this into a nicer form, but that isn't necessary.

Next, to get the the frequency response, replace z with $e^{jw}$

So this yields $-e^{jw} + 1 + 2e^{-jw} + 2e^{-2jw}  + e^{-3jw} -e^{-4jw}$ which is a complex function with both a phase and a magnitude.

Find the phase of this expression and you are done.  I don't have time to wrestle this into Magnitude*phase form, so I'll go to your first example which is a simple case.

$$H[z] =  1 + z^{-4} $$
or
$$H[z] =  z^{-2}(z^2 + z^{-2}) $$
then you have
$$H[e^{jw}] =  e^{-j2w}(e^{j2w} + e^{-j2w}) = 2cos(2w)e^{-j2w} $$

Your phase is the argument of the complex exponential $e^{-j2w}$, so...

$\angle H(e^{jw})$ is the expression $-2w$


How to deal with overlap in STFT when considering the highest energy among multiple channels I have a 4 channels audio file. As a pre-processing step, it goes through an STFT with $T$ frames and $M$ frequencies. The output is divided into magnitude and phase per channel, so for $C$ channels, we have a $T\times M\times 2\cdot C$ tensor.

I wish to focus on the channel with the highest energy per frame $T$. If the STFT had no overlap, I could simply take the highest energy channel per frame and concatenate to a single channel audio file. The $T$ frames, however, are produced by 50% overlapping windows.

In my case for example:

Audio length = 2880000

number of channels, $C$ = 4

Window Length = 1920

hop size = 960

number of FFT bins = 2048

number of frequencies, $M$ = 1024

Hence the output size:

$T = \frac{2880000}{960} + 1$ = 3001

How do I produce a single channel audio file, with the highest energy per window?

For reproduction:
I am using librosa.core.stft for the STFT  with a Hanning window.

The obvious approach to this would be a crossfade between channels. This leads to problems, just as obvious, as you would get nasty comb filter effects at 2000Hz and above. So this is a no-go.

A solution just as simple would be to use just one channel, but apply a gain at each window, so that the energy is aligned with the one of the loudest channel at that point in time, i.e. the gain would be
$$a_i = \frac{\max\limits_{c}(E_i^c)}{E_i^1}$$
for $c$ as channel index, $i$ as time index and assuming you use channel $1$.
This way, you avoid any phase problems whatsoever and still have a signal, that suffices your energy requirement. Also, there is no non-linear processing and thus, no artifacts resulting from it, which would most certainly occur, if you follow some approach with time aligning the signals from different channels to avoid comb filter effects.

EDIT: I know that the overall SNR will not be improved by this approach, but I think the advantages are worth it, since SNR differences should be quite small, as distance between microphones is also small.


Confused about how Echo filter diagram works The picture above is the block diagram of an echo filter. My question is, how does the delay, and gain block affect the output signal? I think that when the input signal is outputted with a delay from when the user and a smaller amplitude.

The output is the sum of the input and a scaled and delayed version of the input

$$y[n] = 2 \cdot (x[n] + 0.8 \cdot x[n-4000])$$

The multiplier of 2 seems a little suspicious. Divide by 2 (or multiply by 0.5) would make more sense.


Bandwidth of a bandpass signal If the Fourier transform of an aperiodic continuous time signal has signal components between the minimum frequency w1 and the maximum frequency w2, but not all the frequencies between w1 and w2, is the bandwidth said to be w2-w1? Or does it require that there be signal components at every frequency in the bandwidth?

The answer to this question could be tricky.
For a simple answer you may take the bandwidth of a bandpass signal as $\omega_2 - \omega_1$ (highest minus the lowest frequencies), and don't care about what happens in between the two. Hence your bandwidth allocation is not fully efficient, and a sampling based on this bandwidth will include redundant information of those empty bands.
Thus, for a more rigorous answer about finding the absolute minimum number of samples to represent a bandpass signal whose spectrum may include empty bands (gaps) between the minimum and maximum fequencies, you would take only the nonzero intervals into count.
Yet when the signal is to be transmitted over a physical channel without frequency allocation control, then you would again consider the minimum and maximum frequencies as the bandwidth, unless cross channel spectrum ovarlapping is allowed.
So the answer depends on your intentions.


Validity of differentiation property of Fourier transform we know the differentiation property of Fourier transform says that, if $$x(t)\longleftrightarrow X(j\omega)$$ then $$\dfrac{d}{dt}x(t)\longleftrightarrow j\omega X(j\omega)$$
We know that we can use this property to find Fourier transform of signum function $sgt(t)$, but cannot be used to find Fourier transform of unit step function $u(t)$.

My argument with intuition is that, we can apply this property only if the average value of the signal is zero. That is, $$\lim_{T_1\to \infty}\frac{1}{2T_1}\int_{-T_1}^{T_1}x(t)dt=0$$

Is this correct or any other validation exists?

Thanks in advance.

You're right about your intuition, but note that the differentiation property is always valid. This can be seen as follows:

If $x(t)$ is a real-valued function which is constant everywhere except at $t=0$ where it has a discontinuity such that $x(0^+)-x(0^-)=1$, then you have

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

In the frequency domain, this is equivalent to

$$j\omega X(j\omega)=1\tag{2}$$

From $(2)$ we can conclude that

$$X(j\omega)=\frac{1}{j\omega}+c\delta(\omega),\quad c\in\mathbb{R}\tag{3}$$

because $j\omega\delta(\omega)=0$.

Note that the real part of $X(j\omega)$ corresponds to the even part of $x(t)$, and the imaginary part corresponds to the odd part of $x(t)$. For the given $x(t)$, the even part is just its DC value. So if $x_0$ is the DC value of $x(t)$, i.e.,

$$x_0=\lim_{T\to\infty}\frac{1}{2T}\int_{-T}^{T}x(t)dt\tag{4}$$

then we must have

$$\mathcal{F}\{x_0\}=2\pi x_0\delta(\omega)=c\delta(\omega)\tag{5}$$

from which we obtain

$$c=2\pi x_0\tag{6}$$

Consequently, the Fourier transform of $x(t)$ is given by

$$X(j\omega)=\frac{1}{j\omega}+2\pi x_0\delta(\omega)\tag{7}$$

For $x(t)=\frac12\textrm{sgn}(t)$ we have $x_0=0$ and consequently,

$$\mathcal{F}\left\{\frac12\textrm{sgn}(t)\right\}=\frac{1}{j\omega}\tag{8}$$

and for $x(t)=u(t)$ we have $x_0=\frac12$ and

$$\mathcal{F}\left\{u(t)\right\}=\frac{1}{j\omega}+\pi\delta(\omega)\tag{9}$$

Also take a look at this and this question, which deal with the same problem in discrete time.


How to match a piece of very short audio based on key and pitch to find a piece in a large database that is most similar to it? For a large amount of musical audio tracks I am doing the following with every single track:

I am slicing the track into very small pieces whereas each piece represents a note onset event. So the piece begins at a note onset event and ends before the next event starts. The result is having hundreds of very small, typically sub-2-second slices for each track.

Now I want to take one of these slices - a totally random pick - and search through all the remaining slices to find the slice that is most similar to the picked one.

Most similar means that it "sounds" like the picked one. More specific this means that it is having the same key, pitch - not timbre and not rhythmically (since rhythm is cut off by slicing at note onsets)

The Problems

To find the slice that is most similar to / that most "sounds-like" a randomly chosen slice I have to determine its "sound". I have been reading about MFCCs but I am not sure if this will help determine the sound, which in this case will be key, pitch.

Next, when having a way to determine the "sound" in numerical way I need to find a way how to compare these results. There are things like euclidian distance or consine-similarity.

Last but not least the slices are of different length.

What I have done so far

Slicing based on note onset-events is done with the help of the onset-feature extraction methods of the librosa python library. This function returns a list of onset-events represented as time-stamps. Each time-stamp is used as a "cut-mark" to slice the audio tracks.

I have been playing with the fastDTW python library which does a dynamic time warping analyzation, compares the result with using an euclidian distance function and fed it with slices. But I am not sure if the result is the distance in MFCCs

The question

To summarize the above into a single question:

How to match a piece of very short audio based on key and pitch to find a piece in a large database that is most similar to it?

well, the "sound" of something depends on a lot of different parameters, which are not all well understood, nor all innumerated.

could be

loudness contour.
brightness contour.
pitch contour.
degree of periodicity vs. inharmonicity.
more parameters that might affect this nebulous thing we call "timbre" (like formant locations or even/odd harmonic ratio).

good luck.


How to deal with signal not equally spaced in time when performing FFT? I wonder to know what is the best way to handle not equally spaced in time signal when performing FFT ?

I guess it depends on the signal itself. I work with signal with about 1000 - 5000 samples and they can contain between 0 to 20% missing samples. And I am interested in both main frequencies AND amplitudes.

Currently I do a basic linear interpolation to "fill" missing samples. I heard somewhere that we could actually "fill" missing samples with NaN (I didn't try yet).

What do you think about these two options ? Advantages and disadvantages ?

It sounds like you are dealing with 'missing' data. Usually in cases like this, you can perform a simple polynomial regression on your data first, and then FFT.

In this way, you will fill in the gaps of the missing data based on the model polynomial you pick. Once you do that, you can perform an FFT. I would try this method first, before considering anything more involved. If you put an image of your data as well with the missing samples, that would go a long way.


What is the ideal FIR length, given a specific pulse length? I am trying to ascertain what the 'ideal' FIR filter length should be, given the pulse length $T_p$ of a windowed sinusoid in noise that I seek to filter.

As parameters into an FIR filter that I design, I have:

$F_c = 15 \text{ kHz}$, The center frequency. (This is the carrier frequency of the signal). I know this.
Since this is a BPF FIR, I specify the pass band as $F_{c} - \frac{1}{T_p}$ to $F_{c} + \frac{1}{T_p}$. This is because the bandwidth of the windowed sinusoid is $\frac{2}{T_p}$
The last parameter that I do not know how exactly to specity, is the length of this FIR... this is where I am lost. What is the ideal length here, (if any?)... Should it just be the length of the pulse (in samples of course), thereby making it something akin to a matched filter? Does this mean I have no further gains in increasing filter length?

As further context, I am seeking this 'ideal' length, should it exist, because I am trying to filter out as much noise as possible, but also try my best to retain the sharp transients. This is what led me to ask, is there an ideal filter length to start from. For example, in the following plot below, I have filtered a noisy version of my signal, with filters of length 11 (red) and 171 (black) respectively. They are shown below:

As you can see, while the black result is 'smoother', you can see that it is also more 'smudged' in as far as its transients go. In contrast, the red still retains some noise, but transients are not as affected.

The plot below shows the spectrums of the above filters:

TLDR: So, is there an 'ideal' length for FIR filters, in so far as that further increasing the filter length does not buy you any more noise immunity, but might actually smear your transients even more than needed?

EDIT:

I have added two new images. The first one has filter of length 11, (red), filter of length 171, (black), and filter of length 901, (blue). The thick-blue is the spectrum of the data.

Here is the corresponding results for the filter of length 11, (red), and the new filter of length 901, (black).

I'm converting some of my comments from above into this answer:

If you are trying to detect the presence of some pulse shape $p(t)$ on the AWGN channel, then the optimum detector uses a filter that is matched to the pulse shape $p(t)$ (aptly named the matched filter); this maximizes the signal-to-noise ratio (SNR) at the output of the filter and thus provides the best detection statistic. This approach is equivalent to a sliding cross-correlation of the pulse shape with the observed signal, expressed mathematically as follows:

$$
d(t) = x(t) * p(t)
$$

where $x(t)$ is the observed signal and $d(t)$ is the resulting detection statistic.

The main problem, therefore, consists of selection of an appropriate threshold that can be used to determine where the pulses of interest occur in $x(t)$. Specifically, one would indicate a detection when $d(t) > T$, where $T$ is a threshold used to balance between two opposing performance metrics: probability of detection $P_d$ and probability of false alarm $P_{fa}$. Here is a link to a previous answer where I talk about the tradeoff a bit more. Target values for these metrics would be chosen according to the requirements for your specific application.

For this case, we can come up with general expressions for $P_d$ and $P_{fa}$ pretty easily:

Probability of false alarm: A "false alarm" indicates a case where the detector reports the presence of the target pulse $p(t)$ when it is in fact not present. Since we've defined the channel to be AWGN, that means that for a false alarm, the input signal $x(t)$ is a white Gaussian noise (WGN) process. Without loss of generality, we will assume that noise to be zero-mean with variance $\sigma^2$.

To determine $P_{fa}$, we are concerned with what the signal at the output of the matched filter looks like. Recall that $d(t)$ is defined as:

$$
d(t) = x(t) * p(t) = \int_{-\infty}^{\infty} x(\tau) p(t-\tau) d\tau
$$

Assuming that $x(t)$ is a white Gaussian noise process with variance $\sigma^2$, it can be shown that $d(t)$ will also be Gaussian, with a variance equal to:

$$
\sigma_d^2 = \sigma^2 \int_{-\infty}^{\infty} |p(t)|^2 dt
$$

That is, the variance at the output of the matched filter is just scaled by the total energy of the pulse waveform $p(t)$. The detection statistic $d(t)$ is then just a Gaussian process with variance $\sigma_d^2$. The probability of a false alarm at any given time $t$ is equal to the probability that the detection statistic exceeds the threshold $T$. Using the properties of the Gaussian distribution, we can write this as:

$$
\begin{align}
P_{fa} &= P(d(t) > T\ | \text{ no signal present}) \\
&= 1 - F_d(T) \\
&= Q\left(\frac{T}{\sigma_d}\right)
\end{align}
$$

where $F_d(d)$ is the cumulative distribution function (CDF) of the Gaussian distribution and $Q(x)$ is the Q-function.
Probability of detection: This case differs from the false-alarm case in that the signal of interest is present. Specifically, we examine the situation where the matched filter is perfectly aligned with the pulse of interest. The same noise component that we analyzed previously is present, but the autocorrelation of the desired pulse shape gives it a non-zero mean. This mean is equal to the total energy of the pulse waveform:

$$
\mathbb{E}\left(d(t)\right) = \int_{-\infty}^{\infty} |p(t)|^2 dt
$$

The probability of detection, therefore, is the probability that the detection statistic exceeds threshold:

$$
\begin{align}
P_{d} &= P(d(t) > T\ | \text{ signal present}) \\
&= 1 - F_d(T) \\
&= Q\left(\frac{T - m_d}{\sigma_d}\right)
\end{align}
$$

The design process would then look like this:

Select an operating range for your detector, defining a minimum signal-to-noise ratio (or equivalently, a maximum noise variance $\sigma^2$) at which you will operate.
Assuming the worst-case conditions (i.e. the maximum noise level), select a threshold $T$ that meets either your required $P_d$ or $P_{fa}$ (whichever is more important to you).
Plug the resulting value for $T$ into the other equation to determine your predicted performance metrics.

This is a pretty high-level treatment of the problem, and if you get down to business trying to build something that works practically, you'll run into some other details of note:

One thing that may be relevant for your problem of detections sinusoids is that it's likely that the received pulses will be at some unknown starting phase $\phi$ relative to your template pulse shape $p(t)$. You would then observe a reduction in the correlation peak based on the amount of phase offset, which will wreak havoc with your detector's performance. If that is the case, a noncoherent detector is a better approach:

$$
d_1(t) = x(t) * p(t)
$$
$$
d_2(t) = x(t) * p_Q(t)
$$
$$
d(t) = \sqrt{d_1^2(t) + d_2^2(t)}
$$

where $p_Q(t)$ is a 90-degree shifted (or quadrature) version of the template sinusoidal pulse. The statistics of this case are slightly different and are left as an exercise for the interested reader.
The above treatment assumes that the pulse is received in $x(t)$ at the same power level as the template $p(t)$, which is almost certain to be untrue. Two ways to approach this complication come to mind: either use some sort of automatic gain control (AGC) process in order to steer the received power level to what you would expect, or you can make $T$ an adaptive threshold that adjusts itself relative to the observed signal (e.g. you could try to estimate the background noise variance $\sigma^2$ and then set $T$ appropriately).


Pitch Detection: avoiding frequency doubling / halving I am working on a tuning app. I have so far tried 3 different libraries with a number of different algorithms. However, when I compare to other tuner apps, I seem to be getting frequency doubling on certain notes with certain instruments. Sometimes frequency halving depending on other factors. I have tried adjusting the buffer size (size of array of samples) and also the sampling frequency.

Algorithms include: Auto-Correlation, YINFFT, Dynamic Wavelet.

Is there sometimes a need to filter the signal before hitting the algorithm?

so Adamski, what you have here is what we call "the Octave Problem" and is well-known and oft discussed with people who work on the pitch detection problem.

first, some definitions.  given a quasi-periodic signal, $x(t)$, where

$$ x(t) \approx x(t+T) $$

in some local region where $t \approx t_0$, we normally define "pitch", measured in octaves as the base-2 logarithm of the local fundamental frequency, $f_0$, of a note relative to some reference frequency $f_R$.  "pitch" is always about an interval, which means it's expressed relative to something.

$$ p = p_R + \log_2 \left( \frac{f_0}{f_R} \right) $$

the local fundamental frequency is the reciprocal of the local period

$$ f_0 \triangleq \frac{1}{T}$$

so if $p = p_R$ it's because the fundamental frequency $f_0$ is equal to the reference frequency $f_R$.  sometimes we use A440 as the reference, most often i use MIDI note 0 as the reference and measure pitch in semitones which are $\frac{1}{12}$ octave.  in the MIDI convention, middle C is MIDI note 60 and A440 is MIDI note 69.  then the definition of pitch (relative to MIDI note 0) is

$$ \begin{align}
p & = 69 + 12 \log_2 \left( \frac{f_0}{440 \text{Hz}} \right) \\
& = 0 + 12 \log_2 \left( \frac{f_0}{8.1758 \text{Hz}} \right) \\
\end{align} $$

so when $f_0$ = 8.1758 Hz (5 octaves below middle C), the pitch is "0" and when $f_0$ is an octave above that, the pitch is "12".

alright, so the basic issue regarding the octave problem is this:

(1)  given that $x(t)$ is fully periodic with period $T$, that is

$$ x(t) = x(t+T)  \quad \quad \forall t $$

then the same is also true that $x(t)$ is periodic with period $2T$ or $3T$ or any integer multiple of $T$.  that is easy to show.  the same is nearly true for quasi-periodic $x(t)$, so there are many candidates of possible periods $T$ to choose from, which one do we choose?  usually the candidate period $T$ that is smaller than all the other candidate periods is the one we choose.  in other words, A440 (A just above middle C) sounds like A440 with period $\frac{1}{440}$ second, but mathematically it is also periodic with period $\frac{2}{440}=\frac{1}{220}$ second, but we don't hear it as A220 (the A just below middle C).

(2)  now here's the rub:  suppose we add just a teeny bit of a 220 Hz waveform to our A440.  but this 220 Hz waveform is, say, 60 dB lower in level than the A440.  what are we going to hear and, mathematically, what is the period of this waveform?

so, your assignment for the day is to think about that and, maybe tomorrow i'll come back and tell you maybe one or two "secrets" (they belong to me, so i think i can tell them, but lacking contracting with me to develop an alg for you, i won't be too explicit) about how to deal with this, but you first have to understand the nature of the Octave Problem which you get a hint of by answering for yourself the question in (2).

stay tuned.


Using continuous verses discrete wavelet transform in digital applications I am familiar with much of the mathematical background behind wavelets. However when implementing algorithms on a computer with wavelets I am less certain about whether I should be using continuous or discrete wavelets. In all reality everything on a computer is discrete of course, so it seems obvious that discrete wavelets are the right choice for digital signal processing. However according to wikipedia it is the continuous wavelet transform that is primarily used in (digital) image compression as well as a large number of other digital data processing activities. What are the pros and cons to consider when deciding whether to use an (approximate) continuous wavelet transform instead of an (exact) discrete wavelet transform for digital image or signal processing?

P.S. (Checking an assumption here) I am assuming continuous wavelet transforms are used in digital processing by simply taking the value of the continuous wavelet at evenly spaced points and using the resulting sequence for wavelet computations. Is this correct?

P.P.S. Usually wikipedia is pretty precise about mathematics, so I am assuming that the applications in the article on Continuous Wavelet Transforms are in fact applications of the Continuous Wavelet Transform. Certainly it mentions some that are specifically CWT so there is clearly some use of CWT in digital applications.

As Mohammad stated already the terms Continuous Wavelet Transforms (CWT) and Discrete Wavelet Transforms (DWT) are a little bit misleading. They relate approximately as (Continuous) Fourier Transform (the math. integral transform) to DFT (Discrete Fourier Transform).

In order to understand the details it is good to see the historical context. The wavelet transform was originally introduced in geophysics by Morlet, and was basically a Gabor transform with a Window that grows and shrinks together with the selected scale/frequency. Later Daubchies (a physicist-ett from Belgium) realized that by choosing special orthogonal wavelet bases the infinitely redundant CWT can be critically sampled on a dyadic grid. From the resulting DWT the corresponding full CWT can be obtained by convolving the DWT with the reproducing kernel of the respective wavelet. The reproducing kernel is the CWT of the wavelet itself.

Daubchies findings gave a big boost to the wavelet theory in the early 80ies. The next big result was that the DWT can be computed very efficiently (this is sometimes called FWT [fast WT] as well) by using techniques from the theory of filterbanks, namely quadrature mirror filters (QMF) together with downsampling filterbanks. By constructing special QMFs the corresponding DWT can be computed via filtering and downsampling, which is the state-of-the-art algorithm to compute DWTs today. You do not need the scaling function to compute the DWT, it is just an implementation detail that FWT process.

Concerning the application side the CWT is the more ideal candidate for signal or time series analysis due to its more fine grained resolution and is usually chosen in most tasks (e.g. singularity detection). The DWT is more of interest in the context of fast non redundant transforms. The DWT has a very good energy compactification and is thus a good candidate for lossy compressions and signal transmissions.

Hope that clarified things.


laser speckle noise reduction I'm trying to process an 1D signal containing a waveform of 100000 pixel.
As it was acquired with a laser scanning system, and therefore it is corrupted by some annoying speckle noise.
I tried to prefilter it with a well-known type of non-linear image processing filter such as Lee and Kuan filters, obviously modified for the 1D nature of the signal.

http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2784081/ (see citation 5,7)

However, the resulting signal is identical to the result of a "naive" mean filter with the same support size (e.g. 15 px).

Is it an "expected" behaviour due to the 1D processing or am I probably writing wrong code?
Can you link me some working code for these algorithms or some other speckle reducing ones?
I found some source code on codeforge, but I think the Lee and Kuan algorithm are not correctly implemented.

Thank you very much in advance.

EDIT:
added some images of the signal:

1) original

2) mean 1x15

3) mean 1x15 with superimposed kuan filter 1x15 with input parameters  nm = 1.0 and nv = 0.1

4) mean 1x15 with superimposed kuan filter 1x15 with input parameters  nm = 1.0 and nv = 0.2

5) mean 1x15 with superimposed kuan filter 1x15 with input parameters  nm = 1.0 and nv = 0.3

We can see in 3,4,5 that there're no differences by mean and kuan: the mean signal is exactly coinciding with the kuan filter one.

Note that at the pixel 50569, where there is an exampe of the speckle I want to reduce, the local mean is 121.666 and local variance is 16.755. Therefore it can be found that an estimation of the noise variance is 0.13771, then I tried the parameters according to said observation.
Further, note that at pixel 50652 the signal is lower and therefore de speckle effect is lower too. There, the local mean is 40.933 and local variance is 4.195 and the ratio is 0.102, therefore in the order of the ratio on an high peak.

Thank you again for your attention

I seems to me like wavelets can offer a neat solution to your problem for noise filtering. Just incase you are not aware a wavelet transform uses two filters (also called mother and father wavelet). These two filters result in a sub-signal (also called sub-band) which is of half the size of the original signal. What I mean to say is say your signal length is 100. Upon wavelet transform you get 2 signals each of size 50. One of these signal is a approximation of original signal and the another one is detail signal. Usually the noise component is present in the detail section.

So as to get rid of the noise, what is usually done is to take a wavelet transform. Implementation of it is available in matlab (dwt() ). This will result in two outcomes cA and cD.

[cA,cD] = dwt(X,'haar')

Then I would suggest you to do usual noise removal algorithms on cD and with cA and cD' (ie. filtered cD) reconstruct the original signal.

cD_1 = someProcessing( cD )

X = idwt(cA, cD_1,'haar')

Please have a look at the following links to get more information on noise filtering with wavelets
1. PhD Thesis titled "Signal Recovery and Noise Reduction with Wavelets" download from HERE.


Why is doing fft on a non-stationary signal a problem? Why is it a problem to do frequency analysis on a non-stationary signal?

what makes the frequency interpretation incorrect?

The frequency interpretation on a non-stationary signal is not "incorrect".  For example, it still works for doing convolution/filtering.  However, it becomes opaque with regard to reflecting features interesting for analysis.

As a simple example, a morse signal of fixed carrying frequency is basically signal times carrier in the time domain and since the carrier is represented by two pulses in the frequency domain, the total frequency representation is an "echoed" variant of the transform of the signal.

However, the signal is interesting as a sequence in time.  Its Fourier transform is an obfuscation rather than a recovery of information.

For many non-stationary signals, the information you are looking for is contained in the sequence of parameter changes over time: an overall Fourier transformation hides those further under the surface than they were before.

There are "number theoretic transforms" also obeying the convolution theorem of the Fourier transform and based on powers of unit roots.  They have FFT-like algorithms and data flow and can be used for convolving number sequences.  They are utterly useless for any form of analysis since their "transforms" invariably look like pseudo-random noise.

In this case, we are talking about something which is "correct" and thus applicable in given situations but which is not useful as an analysis tool.  To a lesser degree, that is the situation when using a Fourier transform encompassing all of a non-stationary signal.


Slow Data rate of GPS data GPS satellites send 50 BPS of info apart from 1 MBPS PRN. What is the basis for choosing 50 BPS. Why can it not be 100 or 1500? The reason I am asking this is Almanac download time is 12.5 mins which is quite long.

The low data rate comes from the need to transmit the data reliably over a noisy channel with a low power transmitter.

A GPS receiver has to be capable of picking out a signal under conditions that other radio systems would just call "total loss of signal."
In a typical two-way radio (analog police radio and the like) the signal is usually unusable if the received signal strength gets below -120dBm.  Digital signals usually get lost long before you get down that low.

GPS receivers work in ranges where a good signal is -155dBW.  In order to fish the data out of the noise, they use a faster signal (1023 million "chips" per second for the civilian codes) that is then modulated with the actual data.  You modulate the data on the chip signal, and modulate that onto the transmit frequency.  You use correlation on the chip signal to recover it, and demodulate the data from that.

The end result works like a very narrow band filter that only lets through data at a low rate, but very clearly.

The chip modulation scheme also has the advantage that you can transmit a bunch of them on a single carrier frequency, but be able to separate them cleanly on the receive side.  This is known as Code Domain Multiple Access (CDMA) and is the reason that the US GPS system only needs a single transmit/receive frequency despite having some 30 satellites.  Glonass uses 15 carrier frequencies (each used twice) to cover the globe with 30 satellites.  It makes the transmitters and the receivers just that bit more complicated.

I would suggest reading the Wikipedia GPS article.  It includes a lot of details I've glossed over in trying to keep this short.


Fourier Transform of ECG signal in Python I have an ECG signal which I am analyzing using Python, as opposed to the mainstream MATLAB. So, I have digital form ECG in .dat file with .hea (header file). Below is the Fourier transform

The problem, as you can see, that it is not the correct Fourier transform. Below is my code.

import wfdb
import matplotlib.pyplot as plt
import numpy as np
from scipy import signal
from scipy import fftpack

cutoff = float(input("Cutoff: "))
cutoff = cutoff/(360.)
noiseAmp = float(input("Noise Amplitude: "))
programed_to = 5000

# load the aami3a ECG signal using special package wfdb for python.
record = wfdb.rdrecord('aami3a', sampto=programed_to, physical=False) # digital data, hence physical=False

# 2D numpy array storing the actual signal data. record is an object of wfdb which has many attributes.
# one of those attributes fs (sampling frequency) and d_signal which stores the amplitude at various times.
original = record.d_signal
# define time of 5000 length
t = np.arange(0, (1/record.fs)*programed_to, 1/record.fs)
# define a sinusoidal noise.
noise = noiseAmp*np.cos(2*np.pi*50*t)

# add noise to the ecg signal
record.d_signal = record.d_signal.flatten() + noise
# make the signal again a 2D array.
record.d_signal = np.reshape(record.d_signal, (-1, 1))

# define a fifth order butterworth lowpass filter with cutoff/fs/2 = 0.11 (tunable).
b, a = signal.butter(5, cutoff, 'lowpass', False)
# apply this low pass filter on the noisy ecg.
output_signal = signal.filtfilt(b, a, record.d_signal.flatten())

# plotting the FFT of the signal
# fft of the original signal.
org_fft = np.fft.fft(original)
freq = np.fft.fftfreq(original.size, 1/720.)
f1 = plt.figure()
A1 = f1.add_subplot(211)
A1.plot(freq, np.abs(org_fft))
plt.show()

I think I have taken Fourier transform correctly. Please help.

Ok, so a little trick and the transform would be correct. The problem is on last 6th line. Instead of np.fft.fft(original), it should be np.fft.fft(np.array(origitnal).flatten()). I think that the wrong results were because original was a 2D numpy array.


Algorithm for finding length and end points of a curly human hair I'm trying to do some analysis of a curly human hair on an iPhone app. I've been using the openCV library.

Particularly I want to be able to estimate the length of the hair, and locate its end points (a measure of curliness is stretched end-to-end length vs. relaxed end-to-end length).

My approach so far has been to use a canny edge detector to find the hair. Then to dilate the result. Finally to skeletonise the image and consider the ends to be those pixels that only have one neighbour. Since the resulting image would hopefully be a one pixel thick line, I could just count the pixels to estimate the length (and maybe adjust for diagonal movements).

The trouble is when the hair is so curly it crosses itself. This makes the skeletonisation add in extra branches and rings that mess up the results.

Is there a better approach for very curly hairs?

Some examples:

Note that once you obtain the skeleton, it is very hard to reverse back to separate the connected components that should not be connected.

The problem is that your original image contrast is too low. I would operate an open morphological operation on your raw image to remove the background, hence increase the contrast of hair.

Your raw image (reverse each pixel by 255-I):

Background ( Open operation with the 20*20 all ones kernel, in order to remove the dark regions, along with the small, thin bright regions, only keep the large bright regions, which are viewed as the noisy background):

Subtract the background from the raw image, and binarize the output, you have:

Next I would do a dilation with the kernel:

[0 0 1;...
 0 1 0;...
 1 0 0]

Remove the small isolated object, you have:

Finally you extract the skeleton:

Looks better than your result.

The most important step is the binarization part. The better you can obtain after thresholding, the easier you will implement the skeleton extraction. And the contrast enhancement is the key operation during binarization process.


Eigenfunctions of discrete-time LTI systems Given a discrete-time signal x[n], how does one determine if it is an eigenfunction of a stable, discrete-time LTI system?

For example, consider

$5^nu[n]$
$2^nu[-n-1]$

Intuitively, second one seems to be an eigenfunction of an LTI system while first one does not. Since the second sequence is absolutely summable while first one is not, is this the criteria for deciding?

i think the only family of functions that are eigenfunctions to LTI systems the sole exponential.

for continuous-time, if it's LTI:

$$ y(t) = \int\limits_{-\infty}^{+\infty} h(u) \ x(t-u) \ du $$

the eigenfunction is $x(t) = A e^{st}$ and the output is

$$ y(t) = H(s) \ x(t) $$

where

$$ H(s) = \int\limits_{-\infty}^{+\infty} h(t) \ e^{-st} \ dt $$

for discrete-time LTI systems, it is similar

$$ y[n] = \sum\limits_{m=-\infty}^{+\infty} h[m] \ x[n-m] $$

the eigenfunction is $x[n] = A z^n$ and the output is

$$ y[n] = H(z) \ x[n] $$

where

$$ H(z) = \sum\limits_{n=-\infty}^{+\infty} h[n] \ z^{-n} $$

often $s=j\Omega$ or $z=e^{j\omega}$ but it wouldn't have to be that.


Comparision between SSIM and MAD Image Quality Assessment Algorithms I have been working with Most Apparent Distortion(MAD) tool to evaluate the quality of images. I have read a paper that compares SSIM, PSNR, FSIM, etc. with MAD. I am uncertain about some calculations done in the paper.

(https://arxiv.org/pdf/2005.01338.pdf)

The output of MAD gives the distance between two images between [0, infinity). How do you convert this value in the similarity of range [0,1], so that it can be compared with other metrics (ssim, ms-ssim)?

A natural practice to convert, in a monotonous way, data in $[0\,\infty)$ to $[0\,1)$ is an increasing function, like those called sigmoid functions. They behave like $x\mapsto x$ close to $0$, and flatten to $1$ when $x\to \infty$. Classical examples are:

$$  x\mapsto \frac{x}{1+x} $$
$$  x\mapsto \frac{x}{\sqrt{1+x^2}} $$
$$  x\mapsto \frac{2}{\pi}\arctan{x} $$

There are many more versions, since they are used a lot in statsitics, articial neural networks, etc.


In cyclic codes, a code vector can be obtained from another by simple circular shift. Why then the other method? In cyclic codes, a code vector can be obtained from another by simple circular shift.  Why then the other method is taught that is a bit round about?  Take the code polynomial. Multiply it by the polynomial with power equal to the circular shift amount and then take the modulo of another polynomial and then convert this resulting polynomial to the codevector to find that the code vector obtained by circular shift and this polynomial method are the same.  Apart from this interesting property is there any advantage of the polynomial method in implementation or in any other way?

While it's true that a circular shift of a codeword is another codeword, it's not true that you can generate all codewords by shifting one codeword. Consider the all-zeros or all-ones codewords; they can't be generated by shifting any of the other codewords.

For very large codes, it is important to have simple algorithms to generate a codeword based on the input. Even if the algorithm seems to operate on a "roundabout" way for small codes, try to implement your idea using a code with $(n,k)=(100,50)$, just to give one example. You'll find that the "roundabout" way is actually much easier to implement!


How to calculate an All-Pass filter quantization coefficients Given

$$H(z) = \frac{z^{-1}-0.8}{1-0.8z^{-1}}\text,$$

does the filter stay an All-Pass filter after quantization of coefficients?

If so, how can I prove it mathematically?

Hint:

Write the transfer with a general coefficient $a$ (instead of the value $0.8$) in the numerator as well as in the denominator.
Compute the magnitude of the numerator for $z=e^{j\omega}$.
Compute the magnitude of the denominator for $z=e^{j\omega}$.
Show that both magnitudes are equal, regardless of the value of the coefficient $a$, hence $H(e^{j\omega})=1$, i.e., the system is an allpass, even after coefficient quantization.


How can I compare two pictures from the camera and tell if there are enough differences to detect movement? I want to use my phone as a CCTV system to control a room in my house, and display an alert when something is moving.

For now, what I have succeeded in doing is to grab a picture every minute and upload them via a PHP script to my server. Now, I would like to compare the current picture and the picture 1 minute ago and detect if someone entered the room. So basically, I would need to compare differences of pixels on the picture (but also taking in account that a cloud may just say hello and changed brightness during one minute)

Does anyone have a clue on how to achieve that or some doc to read?

It seems to me what you are looking for is background subtraction technique. With noisy images and changing lighting conditions it could be nontrivial. The current state of the art technique for this is the low-rank matrix representation, but it require not two but many (~dozen) images. Some heavy duty math follow:
Each image considered as  a vector of pixels, vectors combined into matrix and this matrix is decomposed into low-rank matrix and remnant. Low-rank matrix columns are backgrounds and remnant is moving objects. There are some open-sourced implementations, but only for factorization itself, not complete image pipeline IIRC

Here is a paper and code for matrix factorization
http://www.ece.rice.edu/~aew2/sparcs.html

Here is overview from CS blog and link to other code:

https://sites.google.com/site/igorcarron2/matrixfactorizations

http://nuit-blanche.blogspot.com/search/label/MF

Survey of some other techniques:
http://www.vis.uni-stuttgart.de/uploads/tx_vispublications/Brutzer2011-2.pdf


How to restore an image from its filter output? I'm wondering if there is a way that I can reconstruct an image, given the output vector obtained from the convolution with a certain filter.

For example I convolve my image with a filter from MR8 filter bank, so I got a feature vector v. Which is :

v = img ** filter

so can I define a function that takes v and filter as an input, and restores the img with a certain error ?

img_restored = reconstruct(filter,v);

What have I tried

I searched for deconvolution, I found two things :

Wiener Deconvolution : But it only works for de-blurring (I guess)
deconv() function of MATLAB, but it takes two vectors as inputs, and gives one -complex- number as output.

Edit

Another approach could be storing all feature vectors for pixels (centers of sliding windows) and try to restore only a single pixel, then obtain an image by concatenating these pixels.

Thanks for any help !

This type of problem is called deconvolution.
There are several approaches to tackle it, including Wiener filtering (another name for Wiener deconvolution), but there are many of them.
Wiener filtering is not reserved to deblurring: it is to cancel linear degradation operators.
As remarked by Aaron, it is very often an ill-posed problem and additional knowledge is required to find a good solution.
A classical approach very popular today is to use a regularized inverse problem approach.
The problem is stated as the sum of two terms $E(x) = d(Ax, y) + \lambda r(x)$:

$x$ is the variable;
$A$ is the degradation operator (the forward convolution);
$d(Ax,y)$ is a distance (for example, Euclidean distance) between the degradation of the current solution with respect to the observation $y$. It is called a data term;
$r(x)$ is a prior or regularization term that penalizes solutions $x$ that have a bad shape (usually, solutions $x$ that are not regular enough);
$\lambda$ is some weight that balances the effect of the two constraints.

Then, a good solution $x^\star$ is found by minimizing the functional $E(x)$.
See for example, have a look at Sparse Approximation via Iterative Thresholding for an application of this approach to sparse approximations.


Suitable sampling rate for triangle wave I know this is a fairly simple question but I can't convince myself of the right answer.

For a 5KHz triangle wave, being sample at Fs sampl/second, what's a suitable choice for Fs such that an accurate reconstruction of the cont. time wave can be generated?

My first thought was to satisfy Nyquist, so sample at twice at the wave rate. But, is there something I'm missing here? Like for the square wave, where technically the sample rate is infinite for a perfect reconstruction, is there something similar here?

Thanks!

The 5kHz is the repetition frequency and it's not the maximum frequency of your signal or its bandwidth and your signal isn't band-limited but it has small frequency content over very high frequencies which lead to the Fat32 answer.

but there are other sampling and reconstruction scheme for non band-limited signals like Finite Rate of Innovation (FRI) which could perfectly reconstruct certain families of non band-limited signals like piecewise polynomials (which your signal exactly fits in) with finite sampling rate (in your case 6 sample per each repetition).


Processes: Orthogonal, Uncorrelated, Statistically Independent How are they all related? You can define them as:

Orthogonal Processes: $E[XY] = 0$
Uncorrelated Processes: $E[XY] = E[(X - \mu_x)(Y - \mu_y)] = 0$
Statistically Independent Processes: $E[XY] = E[X] \cdot E[Y]$

If two processes are orthogonal:

they are also uncorrelated
they are not necessary independent

If two processes are uncorrelated:

they are not necessary orthogonal
they are not necessary independent

If two processes are independent:

they are uncorrelated
they are orthogonal

Is that correct? I'm not sure.

You got some definitions wrong. It's correct that orthogonality means that $E[XY]=0$. Uncorrelated means that $X-\mu_X$ and $Y-\mu_Y$ are orthogonal, i.e., $E[(X-\mu_X)(Y-\mu_Y)]=0$. If you work that out you should arrive at the equivalent condition $E[XY]=\mu_X\mu_Y$ for uncorrelatedness (not for independence!). Consequently, if at least one of the two RVs $X$ and $Y$ have a zero mean, then orthogonality implies uncorrelatedness and vice versa.

Statistical independence means that the joint PDF of two random variables can be written as the product of the individual PDFs:

$$f_{XY}(x,y)=f_X(x)f_Y(y)\tag{1}$$

Independence implies uncorrelatedness, but the opposite is generally not true. If $X$ and $Y$ are jointly Gaussian, then independence and uncorrelatedness are equivalent. Consequently, in the special case that $X$ and $Y$ are jointly Gaussian and at least one of them has a mean of zero, then orthogonality, uncorrelatedness, and independence are all equivalent.


Shannon's entropy for non-equiprobable symbols I am trying to understand the concept of Shannon's entropy and deciding the codelength. In first case, b is an array of 5 symbols. The symbol set is $\{1,2,...,8\}$. In general, there could be any integer value between 1 and 8 in b. For this data array, Shannon's entropy = NaN.

clear all
b = [1,3,2,6,1];
p_1 = sum(b==1)/length(b);
p_2 = sum(b==2)/length(b);
p_3 = sum(b==3)/length(b);
p_4 = sum(b==4)/length(b);
p_5 = sum(b==5)/length(b);
p_6 = sum(b==6)/length(b);
p_7 = sum(b==7)/length(b);
p_8 = sum(b==8)/length(b);

ShEntropy =  -p_1 * log2(p_1) - (p_2) * log2(p_2) - p_3 * log2(p_3) -p_4 * log2(p_4) -p_5 * log2(p_5) -p_6 * log2(p_6)...
    -p_7 * log2(p_7) -p_8 * log2(p_8)
%codelength
L = max(- log2(p_1), -log2(p_2), -log2(p_3), -log2(p_4), -log2(p_5), -log2(p_6), -log2(p_7), -log2(p_8))

For the otehr case, I am considering a binary arrays of symbols, d.

 d = [0,0,0];
    p1 = sum(d==1)/length(d);
    Shannonentropy_binary = -p1 * log2(p1) - (1 - p1) * log2(1 - p1)

I can see that when all the symbols are occuring then only there is a legal value for entropy, else it is infinite. But, in real life we cannot ensure that all the symbols would occur in a signal or some other data. So, my questions are

(1)for non-equiprobable sysbols, how can I calculate entropy without encoutering infinity value. Where am I going wrong?

(2) I need help in determining the L, the codeword length or the minimum number of bits required to encode a message. The optimal codeword length / block size should be $HN$ for $N$ number of messages. But, I need to know the value of $L$ which is the block size as I don't know how many number of messages, $N$ would be transmitted. Is there a way to find the exact value of $L$ without using the knowledge of $N$? I had remembered reading somewhere that the block size would be $\max\{-log2(p_i)\} $. What is the correct way?

You need to distinguish between the entropy of your underlying distribution and the "entropy" of the realization.

Consider a set $\mathcal{X}\subset \mathbb{N}$. This set defines, which symbols are allowed to be contained in a random sequence $X=(x_1, \dots, x_N)$. In this sequence, each element $x_i\in\mathcal{X}$ is drawn from the distribution over $\mathcal{X}$ (i.e. there is a function $p(x)$ which gives the probability of occurence for each symbol $x\in\mathcal{X}$). As a constraint, you have $\sum_{x\in\mathcal{X}}p(x)=1$ (i.e. it is certain that one symbol out of $\mathcal{X}$ is chosen for each $x_i$.

Now, the entropy of the random variable for each symbol is defined by $$H=\sum_{x\in\mathcal{X}} -p(x)\log(p(x)).$$
Here, we use the convention $(p(x)=0) \Rightarrow (p(x)\log(x)=0)$, i.e. if the probability for one symbol to occur is 0, it does not influence the overall entropy. This is inline with the limit

$$\lim_{p\rightarrow +0}p\log p=0.$$

So, now to encode a given sequence $b$ (which is drawn from the distribution above), you take each element $b_i$ from the sequence $b$ and associate a bit-word of a length which is (roughly) proportional to its "realization entropy" (sorry, I forgot the common term for this) $p(b_i)\log(p(b_i))$.

So, on average, the encoding for each symbol has $H$ bits. Hence, to encode $N$ symbols, you'd need (on average) $NH$ bits.


Difference between power spectral density, spectral power and power ratios What 'exactly' is power spectral density for discrete signal?

I was always under the assumption that taking the Fourier transform of the signal, and then the ratio of desired frequency range magnitude over the entire frequency range gives the power ratio for that frequency range which is the same as power spectral density. Is that wrong?

Reading a student paper got me confused as it says to compute PSD and then 'absolute and relative spectral powers in desired bands' as well. Are they different? If yes, how does one compute it?

I have no idea what your calculation of power spectral density gives since
I cannot understand it.

If a signal $x(t)$ has Fourier transform $X(f)$, its power spectral density
is $|X(f)|^2 = S_X(f)$.  The absolute spectral power in the band of frequencies from $f_0$ Hz to $f_1$ Hz is
the total power in that band of frequencies, that is, the total power
delivered at the output of an ideal (unit gain) bandpass filter
that passes all frequencies from $f_0$ Hz to $f_1$ Hz and stops everything
else.  Thus,
$$\text{Absolute Spectral Power in Band}
= \int_{-f_1}^{-f_0} S_X(f)\,\mathrm df + \int_{f_0}^{f_1} S_X(f)\,\mathrm df.$$
The relative spectral power measures the ratio of the total power in the band
(i.e., absolute spectral power)
to the total power in the signal.  Thus,
$$\text{Relative Spectral Power in Band}
=  \frac{\displaystyle\int_{-f_1}^{-f_0} S_X(f)\,\mathrm df
+ \int_{f_0}^{f_1} S_X(f)\,\mathrm df}{\displaystyle\int_{-\infty}^{\infty} S_X(f)\,\mathrm df}.$$