Create Matched Filter Reference Signal when Scaling Unknown I am working on a problem which raised what should be a simple question, but I'm stuck on it. The objective is to detect a shifted and scaled version of a known signal in the presence of Additive White Gaussian Noise (AWGN).

Think of a radar system. In a radar system, you know the transmitted signal. The returned system is delayed in time and scaled. I understand that a time-delay, only shifts the location of the matched filter's peak output. No problems there. But if the scaling is unknown, then how does one construct the reference signal?

Using a reference signal with the an amplitude that doesn't match the amplitude of the signal in the noisy data will mean that a threshold tests based on that reference signal will be incorrect.

I'm guessing that in a radar system, the noise power is approximately known, and so the power of the received signal is estimated by taking the power of the received signal and subtracting the noise power. That value could then be used to scale the reference signal. Still, that would only be and approximation, and it would still effect the probability of detection and false alarm in some way. Any help is appreciated. Thanks!

I think the problem is not as bad as you suspect it is.

I wasn't around at the time, but from what I've read, early radar systems essentially connected the matched filter's output to an oscilloscope, and a trained operator would look at the phosphor and decide, from experience and intuition, when the signal raised above the noise ("the grass") indicating a reflection.

In practice, you never know in advance the amplitude of the pulse you're looking for. There is attenuation in the channel that can change over time, even on very short time scales. Sometimes you don't know how much power the transmitter is using.

There are at least two solutions to this problem. One is to use automatic gain control in the receiver, so that you know in advance the power or the peak amplitude of the pulse going into the matched filter. The second is to use a training sequence or a pilot signal, so that the receiver can calibrate the thresholds appropriately.

To add a bit of detail: to simplify the math, assume the signal $p(t)$ is time-symmetric and has energy 1. The transmitter emits signal $kp(t)$, with $k$ unkown at the receiver. The receiver correlates the received signal with a matched filter with impuse response $lp(t)$. When the transmitted signal is present at the receiver, the matched filter's output is $$r=kl+n,$$ which is a Gaussian random variable with mean $kl$ and variance $N_0/2$.

When nothing but noise is present at the receiver, the matched filter's output is just $$r=n,$$ which is a Gaussian random variable with mean 0 and variance $N_0/2$. If the receiver can estimate $l$, then the problem reduces to standard hypothesis testing (see for example these notes).

I would suggest first estimation of the gain/loss by running a window with the size of your matched filter that will be compared to the value of the "filter" with the same window - this will average the noise impact as well (square root of the window size). Use the scale as the window moves to scale your signal and then run the matched filter.


Correct way to implement Discrete and Stationary Wavelet Transform? Can someone please confirm if this is indeed the correct way to compute DWT:

load image
load Haar wavelet

for i = 1: Level
  periodize edges of image
  convolve image with Low pass Haar wavelet filter and store result in G
  convolve image with Low pass Haar wavelet filter and store result in H

  convolve G with low pass Haar to calculate approximation co-eff
  convolve G with high pass Haar to calculate horizontal co-eff

  convolve H with low pass Haar to calculate vertical co-eff
  convolve H with high pass Haar to calculate diagonal co-eff

  i = i + 1
  downsample approximation by 2
end

In case of SWT, instead of down-sampling the approximation, Haar filter co-efficients are up-sampled twice at each level and the same process is repeated. Am I correct here?

In the Haar wavelet case, I just used the following NumPy:

# assume x is passed in as a variable
haar_transform = zeros(N)
i = arange(N/2, dtype=int)
haar_tranform[i]      = x[2*i] + x[2*i + 1]
haar_transform[i+N/2] = x[2*i] - x[2*i + 1]

I believe the mathematical description includes some scaling by $1/\sqrt{2}$.

I've learned that while academic papers may describe something very complicated there's a simple way to compute it.


Sampling rate vs sampling time of FFT I am having some trouble understanding the sampling time and frequency used in the computing of the FFT, and I was wondering if someone could make it clear for me :

The Sampling time: the time in which I took my sample for example 5 minutes, its the difference between the time at which i started taking measures and the time when I was done.
The sampling rate: how many samples do I save per second, which means if I took $N$ sample within $13\textrm{ s}$, my Sampling rate will be $f_s=N/13$.

Is this sampling rate the same as the frequency used when calculating the FFT when creating the frequency axis of the FFT?

The sampling time is the time interval between successive samples, also called the sampling interval or the sampling period, and denoted $T$.

The sampling rate is the number of samples per second. It is the reciprocal of the sampling time, i.e. $1/T$, also called the sampling frequency, and denoted $F_s$.

The frequency axis for the FFT is linked to the number $N$ of points in the DFT and the sampling rate $F_s$. It is defined as $f=k\cdot\frac{F_s}{N}$. With $k$ going up to $N$

When you calculate the Discrete Fourier Transform of a signal, if you have N samples as an input, you calculate two half-sized arrays which are N/2+1 values long.

It implies that when sampling at Fs, you cannot represent frequencies that are higher than Fs/2.

So the frequency axis resulting of a FFT cannot be higher than Fs/2.


From a physics perspective, why does D/A quantization error result in a noise floor? For the last week or so I have been trying to understand how quantization error results in the noise floor outside of a mathematical perspective and I haven't really had any luck finding a source that discussed quantization noise without using equations to show where quantization error comes from.  I certainly have nothing against using math proofs but I was hoping to understand why quantization noise existed rather than how to calculate it.  I think I’ve been able to work through it but I was hoping someone would be able to double check my work or point me in the right direction if I’m missing something.

So here is how I figure it so far:
Because the sample rate is allowing us to have a perfect recreation of the frequency of the wave the bit depth only determines the range of values available for the amplitude of a wave being converted to a digital signal.

When the original signal gets rounded to the nearest bit it creates slight changes in the amplitude which result in a slightly altered waveform at each sample.  This new waveform is basically the same as the sum of our original signal and a new frequency.  So if we were to duplicate this in an analog signal we could do it by adding additional frequencies at very small amplitudes to our signal 44.1 thousand times a second (or whatever the sample rate is).  The lower the bit depth the bigger these changes to our original waveform will be, and so the larger the amplitude of the added waves will be raising the “noise” floor.

If you sampled a waveform that allowed for a slower sample rate, say sampling a 10 Hz wave 20 times a second.  Instead white noise as a result of quantization error would you hear a quick succession of random frequencies paired with your original wave (which would obviously be to low to hear)?

If you can follow my logic up to this point I have one additional question because in my tests changing the phase of my sampled analog wave doesn't seem to change the frequencies of the quantization error like I thought it would. Am I making a mistake in my test?

For the last week or so I have been trying to understand how quantization error results in the noise floor outside of a mathematical perspective and I haven't really had any luck finding a source that discussed quantization noise without using equations to show where quantization error comes from.

That's because it's a purely mathematical effect, there's no physics involved.

It's technically quantization distortion, not noise, but we simplify the analysis by treating it as a noise source.

But it's really not; it's nonlinear distortion, so the spectrum of the quantization error varies with the signal.  Maybe that's where you're confused?

For instance, if you sample a sine wave that's a sub-multiple of the sampling rate and quantize it:

you will get harmonic distortion only at specific frequencies:

If it's not an exact sub-multiple, the infinite number of harmonics produced by the distortion will alias and produce a more noise-like residual:

Since this sort of behavior is more common with arbitrary signals, it's usually treated like a noise source.  In addition, we usually randomize this distortion by adding a small amount of dither noise to the analog signal before quantization, which flattens out the spectrum, reducing the peak value of the distortion and making the noise unrelated to the signal, so we are actually forcing it to be a true noise source rather than a distortion:

from __future__ import division
from numpy import linspace, round, abs, sin, log10, pi
from numpy.fft import rfft
from numpy.random import randn
from scipy.signal import hamming as window
import matplotlib.pyplot as plt

fs = 5000  # Hz
t = linspace(0, 1, fs, endpoint=False)
f = 25  # Hz
x = sin(2*pi*f*t)
plt.figure(1)
plt.plot(x)
plt.figure(2)

x += randn(fs)/20  # Dither
x = round(x*5)/5  # Quantization

plt.plot(20*log10(1/fs * abs(rfft(window(fs)*x))))
plt.ylim(-100, 0)
plt.figure(1)
plt.plot(x)
plt.margins(0, 0.1)

The physics interpretation has to do with the sources for noise. In my experience, if the analog portion of the electronics is "clean", the limiting factor is quantum. Ultimately, you are counting electrons or photons. This gives you a Poisson distribution with standard deviation of the number of particles you are counting. That's the noise floor.

If the electronics is not perfectly clean, it will pick up stray energy. This may be a fluctuation on the zero or full scale reference on the A/D converter. Or, it may have something to do with the detector/transducer. Or it could be longer term variation in the gain of a circuit leading to the A/D, often caused changes in temperature. Or, the detector itself may have a bias that fluctuates.

In my line of work, I often deal with detectors of light. These will often "mistake" a photon of heat for a photon of visible light. They call this dark current. Stray light can also be an issue... catching light that we didn't want to catch.


Building a Signal Model for Samples from a Sensor I have a signal I'm getting from a sensor (shown below in the photo of an excel plot) and I need to process it to get it to look something like the red line I've overdrawn on it.  Sort of a moving average type process that doesn't fall too quickly unless the overall signal is actually falling.

Anyone have any ideas as to what the best method(algo) would be to approach this problem?

Any suggestions on directions I could look towards would be helpful as well, it's been years since I took signal processing in school and I don't remember what methods are even out there that could be tried.

Thanks!

Jonathan

A simple exponential average has the following difference equation:

y(n) = (1-alpha) * x(n) + alpha * y(n-1)

To get a different rise and fall rate, you can dynamically choose alpha depending on whether x(n) is larger than y(n-1). For example, you could build have alpha = 0 when x(n) > y(n-1), and alpha = 0.99 otherwise. This could accomplish the "peak hold" goal, and the second value of alpha would let you control the decay after each peak.

Looks like you're looking for the envelope of that signal. You can actually go back to your original audio signal, load it into Matlab, or something that has the necessary functions, and take the sum of the signal and [imaginary unit] times it's Hilbert transform to get a complex-valued "analytic signal". The magnitude of the analytic signal is the envelope, and should be pretty close to what you are drawing for a red line.

You can read up on it here: http://en.wikipedia.org/wiki/Analytic_signal


Remove brighter regions in an image on frequency domain I'm newbie and maybe this is not the right place to ask this question, or this question has no sense.

I'm very interested in working on filtering images on frequency domain. I have found a few filters on frequency domain and I wondering if it is possible to remove bright regions on an image (with white colour) using filters on frequency domain.

Imagine you have a greyscale image and you want to remove the pixels with white colour.

How can I do that?

It is important to understand what the frequency domain means. The bright regions in an image correspond to pixels with a certain value, whether it be color or gray scale. Frequency domain has the interpretation of looking at how fast the values change, not their values. For example, think about an image of a dark blue car against a white background. The pixels along the edge may change from dark blue to white just by moving a few pixels, this is high frequency because the change in pixel value is high over a just a few pixels. This is different from a picture where you have dark blue which slowly fades to white over many pixels. Filtering in the frequency domain on an image would be affecting how the pixels change in value. For example, using a low pass filter on the dark blue car picture I mentioned, would block the high frequency components. This means the sharp changes in pixel value, like on the edge between the dark blue car and white background, would be smoothed (blurred) out.

Is it possible to remove white pixels from an image with frequency domain filter? No, and if you need to do that it would be better to just scan the image for pixels with a white value or value above some threshold and remove those pixels.


Conjugate reciprocal pairs of zeros and poles in FIR design Assuming the impulse response $h[n]$ of an FIR filter is real for all $n$,

Why are zeros and poles in FIR design found in reciprocal and conjugate pairs?
Is the assumption necessary for this phenomenon to take place?

Look up the complex conjugate root theorem which states that:

  If all the coefficients of a polynomial are real then its roots are either real or if there is a complex root, then its conjugate is also a root.

This theorem can be applied to the denomerator and numerator of a rational transfer function to judge about its poles and zeros.

"conjugate reciprocal roots" is not a necessary consequence of real coefficients. conjugate reciprocal pole-zeros are considered when we need some specific features such as being all-pass or when designing minimum-phase filters. This is because for instance reflection of a zero to its conjugate reciprocal does not influence on the magnitude response.

If the filter has a finite impulse response (and does not come from a truncated IIR filter, as stated in the comments by Robert), then all its poles must be at the origin. So no point stating that its poles are reciprocal and conjugate pairs: they are all at $z=0$.

If $h(n)$ is real, causal and stable, then its zeros are found in reciprocal and conjugate pairs if and only if the filter has (generalized) linear phase.


Is it possible to generate Hamming code with minimum distance $d$ which can correct more than 1 bit error? In literature, all I am seeing is Hamming code with minimum Hamming distance $3$. Theoretically a Hamming code with minimum distance d can detect $d-1$ errors and can correct $(d-1)/2$ error. So minimum distance of $3$ can detect $2$ errors and correct $1$ error.

So here is my question: Is it possible to generate a code with minimum Hamming distance of (let's say) $5$ which can correct $2$-bit errors in each block of coded bits?

If yes, could anyone please give me suggestion about that?

If not, then please give me an explanation.

You cannot detect 2-bit errors (reliably) with a Hamming code distance of 3 if you are correcting 1-bit errors.

At any rate, you can design Hamming codes correcting two errors.  Let's assume 32bit of data and n Hamming code bits.  The states we then need to be able to encode is
0 bit error = 1 state
1 bit error = 32+n states
2 bit error = (32+n)(31+n)/2 states
For n=10, this amounts to 1+32+10+21*41 = 904 states that can be encoded with 10 bits of redundancy.  Depending on how complex you want your correction machinery to get, you might want to use more bits in order to make multiple error correction less complex.


Am i performing Sinesweep correctly? I am bit unsure whether i am doing i sine sweep correctly.
I am feeding my system a sine wave with f hertz within a given time interval and log its output.

And then i do the same for another frequency.

I tried this technique but seem to do it incorrect or so i think.
I tested for 2 frequencies, and fed them into  Ident in matlab and got a fit rate of 98 % which makes no sense that it could build a such a good bode plot based on 2 values..

I am doing something wrong.. should it not have more than one 2 input or is it enough?

The theory behind sweep-sine measurements of LTI systems requires a signal with constantly changing the frequency. You cannot simply playback few tones - the whole frequency range is necessary.

So that if you want to identify your system with the impulse response $h[n]$, you feed the sweep sine signal $s[n]$ into it and record the output. Obviously output will be given by the convolution:

$$y[n]=h[n]\star s[n]$$

In order to obtain the impulse response of the system, it is required to convolve the output $y[n]$ with the inverse filter $f[n]$. Generally, the inverse filter is a sweep signal, inverted in time and with scaled amplitude. So in the end:

$$h[n]=y[n]\star f[n]$$

There are two major types of sweep signals: linear and logarithmic (exponential) frequency change. The latter one is more commonly used, since its spectrum is similar to a pink noise. Personally I always use the logarithmic sweep sine:

$$s(t) = \sin\left[ \dfrac{2\pi f_1 T}{\ln\left( \dfrac{2\pi f_2}{2\pi f_1} \right)}\left(e^{\dfrac{t}{T}\ln\left(\dfrac{2 \pi f_2}{2\pi f_1} \right)} -1\right) \right]$$

where:

$f_1$,$f_2$ - initial and final frequency of a sweep

$T$ - sweep duration in seconds

There is no strict rule on what should be the length of a sweep. For example in acoustics, the rule of thumb is that it should be at least as long as predicted reverberation time.

For more literature on the topic, please refer to the following sources:

  Meng Q., et al. - Impulse Response Measurement With Sine Sweeps and
  Amplitude Modulation
  Schemes

  Farina A. - Advancements in impulse response measurements by sine
  sweeps

  Farina A. - Simultaneous Measurement of Impulse Response and Distortion with a Swept-Sine Technique

A sine sweep usually goes over the entire frequency range of interest. For example with audio you may have a sweep of 20 Hz to 20 kHz (or 16 kHz). All of the frequencies will be swept through in one continuous output. You can do this by generating the samples with sin(2*pi*n*f/fs) with f incremented every arbitrary number of samples n.


Percentage difference between two signals Apologies in advance for what I imagine is a simple question, however I can't work out how to do this.

I have two signals of the same type but with different sensor types, the figure below will give you an idea of some of the data I've got:

I very simply want to calculate the average percentage difference between the two signals. Methods I've attempted so far have calculated outrageous results, in the order of ~200% difference. My attempt is below:

  for i = 1:length(data)
      per(i) = abs((var1(i)-var2(i))/((var1(i)+var2(i))/2)) * 100;
  end

I believe the problem in this method is that I'm only calculating the difference between the two immediate points in question (i) and not on a grand scale of all the data.

Any help would be greatly appreciated.

At the least you should normalize it by the length of the sequence. But what you're trying to do is a pretty standard thing. I'm not sure what percentage difference really means. What you should probably be using is something like the root mean squared error:

$$RMSE=\sqrt{\frac{1}{N}\sum_{n=0}^{N-1}(x(n)-y(n))^2}$$

where $x(n)$ and $y(n)$ are your two sequences.

You can use this command:

plot(T,100*(var1-var2))

Hope it might solve your problem although you don't need a plot command.


Implementing a phase-neutral band-pass filter using SciPy Does the following code properly implement a phase-neutral band-pass filter using SciPy?

freq_fraction = FREQ / (SAMPLE_RATE/2)
sharpness = 0.1
filter_length=128
fir = sig.firwin(filter_length, [freq_fraction*(1-sharpness),freq_fraction*(1+sharpness)], pass_zero=False )
filtered = sig.lfilter( fir, 1, data_to_filter )
filtered = sig.lfilter( fir, 1,filtered[::-1] )[::-1]

The code is based on the answer to:

Is it correct to subtract a low-pass filtered signal from the original signal and use the result as a "high-pass"?

which doesn't give a final code listing for how to maintain phase neutrality.  I'm assuming it's just a matter of running the filter forward and backwards but want to double-check.

Also, if there are better/different ways of doing this (other than the remez and firwin2 methods listed in the other question) I'd love to know.

How about scipy.signal.filtfilt and scipy.signal.butter filter? See the example at the bottom of the documentation.

Here's another example with and without phase neutrality: Applying filter in scipy.signal: Use lfilter or filtfilt?

Calculate the frequency response (https://en.wikipedia.org/wiki/Frequency_response) of the filter and see yourself.


Draw a frequency response graph from a simple equation How is it possible to draw an audio frequency response graph (x-axis : freq from 20 to 20 000hz, y-axis : dB)  from a simple equation like

y[n] = x[n] - x[n-486]

(Example: with this equation I get a nice comb filter that removes frequency 197,5Hz and all its harmonics, for an audio file with sampling rate 96 Khz)

How to turn such an equation into a graphical interpretation like this
http://downloadfreesamples.com/wp-content/uploads/2013/08/EQ-Four-CIS-DSP_3.jpg ?
Then I could see how deep / steep the filtering is.

Thank you

(PS : I continue to read references and online books about filter design that were given to me in previous questions, thanks for these references by the way.)

The tool of choice here is the z-transform, which is applied to your equation. It transforms sequences of complex numbers to functions in the complex plane. Its most useful properties are that it takes delays to powers of the argument z, and for z = exp(i omega) we get the fourier transform of the transformed sequence. Look up the details, I'll just show you how it works here.

y[n] = x[n] - x[n-486]  ---> Y(z) = X(z) - X(z)*z^(-486)

The equation in z-domain factors on the right hand side to give Y(z) = X(z) ( 1 - z^(-486) ).

Since Y(z) is the output and X(z) is the input, for linear systems we define the transfer function to be Y(z)/X(z) = H(z) and the above equation becomes

H(z) = Y(z)/X(z) = 1 - z^(-486)

Now if we substitute z = exp(i omega) we get the frequency response of your system, which is

H_f(omega) = 1 - exp(i omega)^(-486)

and with the properties of exp() it simplifies to

H_f(omega) = 1 - exp(- 486 i omega)

If you need actual frequency f instead angular frequency omega, just substitute omega = 2 pi f/fs, with sampling rate fs.

Finally, you'll probably want to know just the magnitude of the response, so look at abs(H_f).
In this specific case it's very simple to calculate because abs(H_f) = sqrt( H_f * conj(H_f) ) and conj(H_f) = 1 - exp(+ 486 i omega).

As requested, some code. Here's what you can do in matlab/octave to get the magnitude frequency response with a logarithmic frequency axis and magnitude axis in dB.

% Create a logarithmically spaced frequency axis with range 20Hz to 20kHz

fAxis = logspace( log10( 20 ) , log10( 20000) , 10000);

% calculate the complex frequency response at these frequencies

% Set SampleRate before calling this, like SampleRate = 44100;

Hf = 1 - exp(1i*486*2*pi*fAxis/SampleRate);

% plot the magnitude spectrum

semilogx(fAxis,10*log10(Hf.*conj(Hf)));

I finally found a solution (probably not the best, but ...)

1) Generate a 10 seconds white noise (which has flat frequency spectrum on 10 seconds average)

2) Apply the filter

3) Display the spectrum in a Sound editor software (e.g. Sound forge) with averaging on the 10 seconds.


How to understand FFT results of scipy.fftpack? I calculated FFT for a speech wav-file using scipy.fftpack. How do I read (understand) the return of FFT? I have read that it supposed to be like so: y[0] is 0Hz loudness, y[1] is 1Hz loundess, ... y[n] is nHz loudness ... But seems like it is not like that exactly.

Q1: What will I get when I do abs(y)? I know that we get list of complex numbers from FFT and need to square() or abs() them. But what we will have after that? Is this Decibels?

Q2: Why do we need normalize wav-data before doing FFT? What does depend on this? Before and after normalization I get different results from FFT. If I do normalization, then results of FFT are measured by hundreds, if I don't the results are measured by 1.x small values... Is this Decibels also?

# Read wav-data
fs, data = wavfile.read('eric.wav')

# this is a two channel soundtrack, I get the first track
wavdata = data.T[0]

# this is 16-bit track, b is now normalized on [-1,1)
wavdata = wavdata / (2.0**15)

Q3: What is length of the returned list from FFT? Seems like the length of the result depends on length of given sound file... But in Q1 I supposed to get list of frequencies and their loudness independently from a given source of data. For now, if I cut in half wavdata I will get twice shorter resulting list from FFT...

Complete simple code:

import matplotlib.pyplot as plt
from scipy.fftpack import fft
from scipy.io import wavfile

# load the data
fs, data = wavfile.read('eric.wav')

# this is a two channel soundtrack, I get the first track
a = data.T[0]

# calculate fourier transform
y = fft(a)

# show
plt.plot(abs(y), 'g')
plt.show()

Q4: What do I process results form FFT to get it in form Db vs Hz?

Wav-file could be found here: https://aacapps.com/lamp/voices Thanks.

One rarely takes an FFT of a whole audio file, or a half. Usually only a sequence of shorter segments some number of milliseconds in duration.  The length of the "list of frequencies" will be half of the FFT length (for strictly real data) and the frequencies related to each y[i] will then be proportional to the sample rate and inversely proportional to your chosen FFT length.  Thus y[1] will only represent the output of a leaky narrow-band filter centered at 1 Hz if the FFT length is equal to the sample rate.

Decibels are more related to the log of the abs() or log magnitude of the FFT results. But without a calibrated microphone and DAC, the resulting decibels will be relative to some arbitrary level or offset. Normalization merely changes the arbitrary offset.


Sample time in convolution I have a signal sampled at a period of 0.01 and a channel impulse response sampled at 0.1250. I convolved the signal with the  channel impulse response.

It is right to perform the convolution between two signals with diffrent sample times?
If it is,  how can I calculate the sample time of the convolution result?

It is right to perform the convolution between two signals with diffrent sample times?

No, aside from a very few special use cases (if you have to ask, not yours), you need to resample either of them so that both are at the same sample rate.

Also notice that sampling a channel impulse response at 8 Hz only gives you information about these 8 Hz (or 4 Hz, depending on whether this is real direct sampling or a sampled equivalent complex baseband representation of the channel). You're trying to apply it to a 100 Hz (or 200 Hz, depending...) wide channel: that won't give you anything useful, even if you resampled. Your channel impulse response measurement simply isn't sufficient.

Convolution does not change the sample rate, but the sample rates of the two inputs need to match. A simple way to do this in your case is to zeropad the spectrum of your signal by a factor of 2 and the spectrum of the channel impulse response by a factor of 25. That way they will both have a sample period of 0.005.

EDIT: There is nothing wrong with this answer. Zeropadding the spectrum of a signal is a simple and common way to change the sample rate. Below is an example using octave.

% FFT function with swapping.
sfft = @(x) fftshift(fft(ifftshift(x)));
sifft = @(x) fftshift(ifft(ifftshift(x)));

% The frequency sample locations for a given
% sample rate and DFT size.
fvec = @(N,fs) (0 : fs/N : (fs-fs/N)).' - (fs - mod(N,2)*fs/N) / 2;

% Parameters for the input signal.
fs1 = 100;                % Sampling frequency (Hz)
T1 = 100;                 % Signal duration (s)
t1 = 0 : 1/fs1 : T1;      % Sample times (s)
N1 = length( t1 );        % Number of input samples
B1 = 0.8 * fs1;           % Signal bandwidth (Hz).
f1 = fvec(N1,fs1);        % Frequency sample locations (Hz)

% Define the some bandlimited random signal.
x1 = sifft( exp(1j*2*pi*rand(N1,1)) .* ( abs(f1) <= B1/2 ) );

% Parameters for the channel impulse response.
fs2 = 8;
T2 = 20;
t2 = 0 : 1/fs2 : T2;
N2 = length( t2 );
B2 = 0.8 * fs2;
f2 = fvec(N2,fs2);

% Some random numbers to represent the channel response.
x2 = sifft( exp(1j*2*pi*rand(N2,1)) .* ( abs(f2) <= B2/2 ) );

% To convolve the input signal with the channel
% impulse response, they both need to be at the
% same sample rate. We do this by zeropadding the
% spectrums of each such that both are sampled
% at the same rate. We use the least common multiple
% so that the zeropad factor is an integer.
fout = lcm( fs1, fs2 );

% The number of samples in the zeropadded spectrums
% for both the input signal and the channel response.
Nzp1 = N1 * fout / fs1;
Nzp2 = N2 * fout / fs2;

% The indices that correspond to the lower
% frequencies. These are the bins into
% which we copy the specrtums.
zinds1 = fvec(N1,N1) + floor( Nzp1/2 ) + 1;
zinds2 = fvec(N2,N2) + floor( Nzp2/2 ) + 1;

% Insert the spectrum into an array of zeros
% and inverse FFT to get the signal at the
% common sample rate.
tmp = zeros( Nzp1, 1 );
tmp( zinds1 ) = sfft( x1 );
x1_out = sifft( tmp );

tmp = zeros( Nzp2, 1 );
tmp( zinds2 ) = sfft( x2 );
x2_out = sifft( tmp );

% Now that both signals have the same sample rate,
% they may be convolved.
y = conv( x1_out, x2_out );

% The frequency sample locations of the convolution.
fy = fvec( length( y ), fout );

figure();
plot( fy, abs( sfft( y ) ) );
xlabel( 'Frequency (Hz)' );
ylabel( 'Amplitude' );
title( 'Spectrum After Convolution' );
xlim( [ -10 10 ] );


Detecting changes across two images I am trying to calculate the change in pixels across two images.
Visually looking at them, there seem to be too many changes but calculating the distance between pixels gives me around about the same percent change values.

I am trying to calculate the distance here by adding the squared differences for each pixel and its square root which gives me the distance of value at rgb level for each pixel. i.e how much has the pixel value changed (sqrt((img1.r - img2.r)**2+(img1.g - img2.g)**2+(img1.b - img2.b)**2).

Then I am trying to check how many such pixels have changed by a threshold of 10.

I am not able to figure out, why does it show about the same difference in any two pair of images.

Is there a problem in the program or have I been thinking incorrectly

import pylab,skimage
import imageio,math,numpy as np

img1 = imageio.imread("029.jpg")
img2 = imageio.imread("028.jpg")
img3 = imageio.imread("030.jpg")
wh = img1.shape[0]*img1.shape[1]
squared_diff = ((img1 - img2)**2)
squared_diff1 = ((img3 - img2)**2)
squared_diff2 = ((img1 - img3)**2)
print float((np.sqrt(squared_diff[:,:,0] + squared_diff[:,:,1] + squared_diff[:,:,2])>10).sum())/wh*100
print float((np.sqrt(squared_diff1[:,:,0] + squared_diff1[:,:,1] + squared_diff1[:,:,2])>10).sum())/wh*100
print float((np.sqrt(squared_diff2[:,:,0] + squared_diff2[:,:,1] + squared_diff2[:,:,2])>10).sum())/wh*100

Output :

45.2615017361
46.1771918403
46.2858072917

For higher thresholds like 20 it gives 0 . np.max gives max distance as 15.x  everytime.

Is there any other way to find a feature like number of pixel changed with a threshold?

Any help is appreciated.

The value of :

img1 = imageio.imread(".jpg")

are rgb value ? i mean are they ranged between 0-255 or between 0 and 1.
Cause if they where all between 0 and 1 then:

((img1 - img2)**2)

will have always value below 10 indeed the maximum value of this operations could be only 1 and therefore the maximum value obtained here :

np.sqrt(squared_diff[:,:,0] + squared_diff[:,:,1] + squared_diff[:,:,2])

could be np.sqrt(3) and that's always minor than 10 and ">10" operator will return an array with all False values


FFT of any sinusoidal signal -- how many point fft to take? suppose I have to measure the DFT of an 60 hz sinusoidal signal and calculate upto 20 harmonics of this.  Value of the 20th harmonic will be 20*60 = 1200 HZ. Then my sampling frequency should be >= 20*60*2 >= 2400 hz

Now how many point FFT i have to calculate to estimate
information about harmonics ?
Is number of points in FFT(N) equal to the sampling frequency ?
How number of points in FFT(N) related to the sampling frequency
?
Is there any advantage with output off FFT if we sample at
higher rate ?

Is bellow three estimates are correct ?

1> If suppose my sampling frequency is fs =  20*60*2 = 2400 then, What i know is that the frequency range [0,fs] is represented by 2400 point.
   I.e each point represent 2400hz/2400 = 1 HZ information. Then output of FFT will be, really only need to look at the first N/2+1 samples

   X[0] is the constant term
   X[1] is 1 hz
   X[2] is 2 hz
   X[1201] is the 1200 HZ

2> If suppose my sampling frequency is fs =  20*60*2 = 2400 then, What i know is that the frequency range [0,fs] is represented by 2400 point. can  i take 4800 sample points ?
   I.e each point represent 2400hz/4800 = 0.5 HZ information. Then output of FFT will be, really only need to look at the first N/2+1 samples

   X[0] is the constant term
   X[1] is 0.5 hz
   X[2] is 1 hz
   X[2400] is the 1199.5 HZ
   X[2401] is the 1200 HZ

Is there any advantage here to take more number of samples to calculate the FFT ? One thing i note is that spectrum is more dense with each point having 0.5 HZ of information, will my calculation of FFT will be more accurate in this case ?

3> If suppose my sampling frequency is fs =  20*60*4 = 4800  =  then, What i know is that the frequency range [0,fs] is represented by 4800 point.
   I.e each point represent 4800hz/2400 = 1 HZ information. Then output of FFT will be, really only need to look at the first N/2+1 samples

X[0] is the constant term
   X[1] is 1 hz
   X[2] is 2 hz
   X[2401] is the 2400 HZ

Now in this case FFT is generating information upto 2400 hz but i need upto 1200 hz. Is there any advantage with output off FFT if we sample at higher rate ?

Please suggest on this.

Many of your questions are answered elsewhere, so I added additional links below but to summarize in one place relative to your questions I offer the explanations below that together with the links should answer your questions:

Required Sampling Frequency:

If the signal is real, then yes the sampling frequency must be greater than twice the highest frequency for a signal bandwidth that extends down to DC. However if the signal is complex then the sampling frequency should be greater than the highest frequency for same case. This is because a complex signal has a unique spectrum over the Nyquist frequency range of $\pm f_s/2$ where $f_s$ is the sampling rate, while for a real signal the negative frequencies are the complex conjugate of the positive frequencies (identical magnitude, opposite phase).

Note that prior to sampling an analog anti-alias filter is required on the signal, otherwise energy if present in higher Nyquist zones (frequencies above the Nyquist frequency range defined above) will alias into the Nyquist frequency range being observed.  (see this post which explains how aliasing occurs: Aliasing after downsampling). It is not feasible to implement a brick-wall filter, so additional margin to the sampling rate is usually added to allow for the characteristics of the anti-alias filter. For this you need to understand the properties of the signal signal of interest prior to sampling (as once sampled there is no way to distinguish an alias from a signal in the first Nyquist zone!). For example in your case, where you are interested in observing the signal out to 20 Harmonics of 60 Hz: The anti-alias filter should pass the 20th harmonics without distortion and then reject the higher harmonics so that you can minimize the sampling rate needed. You would then choose the sampling rate depending on your filter implementation such that no energy folds in to distort the 20th harmonic.

How Many Point FFT is Needed:

This has been answered several times, so see the links given below for further details and caveats, but the main point is that the frequency resolution of each "point" is 1/T Hz where T is the length of your signal in seconds.  Further the bins in the DFT (the FFT is just an algorithm that computes the DFT) as typically returned extend from DC to one bin less that $f_s$. Due to aliasing the DC is equivalent to $n f_s$ where n is any integer, and likewise for all other frequencies in the DFT, so the DFT that extends frmo DC to $f_s$ equivalently extends from DC to one bin less than $f_s/2$ and then from $-f_s/2$ to one bin less than $f_s$ (see fftshift in Matlab). As detailed in the links given, due to the relationship of the sampling rate and the total length in time T of your signal, each bin will have a frequency resolution of $1/f_s$ when no windowing is used.

Note importantly that each bin of the DFT is the result of a Sinc filter, so if your sampling rate is not an integer multiple of 60 Hz, you will also have the effect of "scalloping loss" that will distort the measurement of each harmonic. (see the "Bank of Filters" explanation at this link: Intuition for sidelobes in FFT). By choosing a sampling rate that is an integer multiple of 60 Hz, each 60 Hz harmonic will be properly centered in each bin.

Suggestion for your specific application:

With relationship to choosing the sampling rate and number of bins for your specific application, if the 60 Hz and its harmonics are the dominant signal and you don't have other signal energy of significance then there is a simple solution that will work quite elegantly in accurately estimating each harmonic with minimum measurement distortion. As described above, choose a sampling rate that is an integer multiple of 60 Hz, high enough to allow for the anti-alias filter used, and then choose the number of points such that each bin passes a harmonic of 60 Hz and rejects all the other harmonics. This means the number of samples in the FFT is $f_s/60$ where $f_s$ is an integer multiple of 60 Hz > 2400 Hz (assuming your signal is real). For example 4800 Hz may be a good choice as you suggest in your options (allowing for a fairly relaxed anti-alias filter!), and in this case $N =4800/60= 80$ points.  So in this case $f_s=4800 Hz$ and with 80 points each point is 60 Hz, so 0, 60, 120, 180 etc centered on each harmonic and perfectly nulling the other Harmonics. The actual magnitude measured for each point is the total area under the Sinc function frequency response for each bin, so if there is other energy present at non-harmonic frequencies, this would distort the magnitude measured in each bin. But if the harmonics are dominant, this would be an excellent and simple approach.  This is somewhat similar to the objective in this link in removing harmonic signals (Different way to separate a particular frequency from a signal) as each bin in the DFT is a rotated moving average filter just as done in this referenced link when only one bin is of interest. Note from that link how every other harmonic is nulled!

Please see these other responses which will further help give insight:

Specific frequency resolution

FFT frequency resolution

Number of FFT points required for a specific frequency resolution for an oversampled signal

FFT and number of samples relations

What happens when N increases in N-point DFT

Suppose 20th harmonic is 20*60 = 1200hz

1> If fs = 20*60*2 = 2400hz
now Taking N=2400  or 4800 points won't change the frequency range you're looking at as long as you have the same  fs...
If N = 4800, then output of FFT will also be N = 4800
and output of the fft will be represent  maximum spectral information upto fs/2  ..
and each index of FFT will represent frequency increment of (fs/2)/(N/2) = 1200/2400 = 0.5
& Spectrum information which can be get from this FFT is upto fs/2 = 2400/2 = 1200

   X[0] is the constant term
   X[1] is 0.5 hz
   X[2] is 1 hz
   X[2400] is the 1199.5 HZ
   X[2401] is the 1200 HZ

2> If fs = 20*60*4 = 4800hz
now Taking N=4800  or 9600 points won't change the frequency range you're looking at as long as you have the same  fs...

If N = 9600, then output of FFT will also be N = 9600
and output of the fft will be represent maximum spectral information upto fs/2  ..
and each index of FFT will represent frequency increment of (fs/2)/(N/2) = 2400/4800 = 0.5
& Spectrum information which can be get from this FFT is upto fs/2 = 4800/2 = 2400 hz

   X[0] is the constant term
   X[1] is 0.5 hz
   X[2] is 1 hz
   X[4800] is the 2399.5 HZ
   X[4801] is the 2400 HZ

If N = 4800 , then output of FFT will also be N = 4800
and output of the fft will be represent maximum spectral information upto fs/2  ..
and each index of FFT will represent frequency increment of (fs/2)/(N/2) = 4800/4800 = 1
& Spectrum information which can be get from this FFT is upto fs/2 = 4800/2 = 2400 hz

   X[0] is the constant term
   X[1] is 1 hz
   X[2] is 2 hz
   X[2401] is the 2400 HZ


A 23 kHz sine wave sampled at 48 kHz, then FFT A 23 kHz sine wave sampled at 48 kHz with 100 samples. Then do a 64 point FFT of the sampled signal. Why is the FFT of the signal not a delta function at 23 kHz. Or why does the  FFT also have values at other frequency?

I think it has something to do with 23kHz / 48 kHz * 64 samples not being an integer. I just couldn’t think the logic through. Please help me! Thank you!!

Edited: sorry guys, i meant 23 kHz

Since you used 2.3k twice, I am assuming you meant that with the original 23k value.

Converting your frequency:

$$ 2.3k \frac{cycles}{second} \cdot \frac{ 64 \frac{samples}{frame} }{ 48k \frac{samples}{second}  }  \approx 3.07 \frac{cycles}{frame}$$

Because this is not an integer value, the peak DFT value will be at bin 3 and there will be values in all the other bins, which is commonly called "leakage".  Being a real valued signal, the upper half of the DFT will be a mirror image with complex congjugate values, so there will also be a peak at bin 61 (64-3).

Your 100 sample figure is irrelevant if you are only using 64 of them for the DFT.

Signals have to be a whole number of cycles per frame in order for the DFT not to have leakage.  The other answers presume you understand that the DFT is a dirac delta, windowed, convolved version of the continuous FT.  This is not an easy way to understand it, but it is how it is usually taught.  It also leads to the misconception that the leakage is described by what is known as the sinc function $ \frac{\sin(d)}{d} $.  That is not true.  The fall off for a complex valued signal is described by what is known as the aliased sinc function $ \frac{\sin(Nd)}{N\sin(d)} $.  Another name is "Dirichlet function".  A real valued signal is actually the sum of two complex signals, so that function doesn't apply either.

You can find the exact formula for the bin values, the peak and "leakage", in my blog article "DFT Bin Value Formulas for Pure Real Tones".

You can also find a much more intuitive explanation of why leakage occurs in another one of my blog articles "DFT Graphical Interpretation: Centroids of Weighted Roots of Unity".

Hope this helps,

Ced

Short windows require a much higher sample rate.  Sample rates only slightly above the Nyquist rate only work for near infinitely long FFTs.  Also, there will be windowing artifacts for any data waveforms that are not exactly integer periodic within the FFT aperture.


Can $\delta(t+\infty)$ be a legitimate signal? Mathematically speaking, when I try to use some signal to disprove a system is invertible, can I use the signal like $\delta(t+\infty)$ ($\delta$ representing the Dirac distribution)? For example, the two inputs $0$ and $\delta(t+\infty)$ are distinct, but their outputs from the system $\int_{-\infty}^{t}e^{\tau}x(\tau)d\tau$ will then both be $0$. In this way, can I say the system is invertible?

It's possible, in mathematics, to complete real numbers with "infinite" values, with sound topological properties; for instance non-standard analysis or the Extended real number line (discussion at Math StackExchange).

However, in this context, for any standard real number $t$, the rule is to set $t\pm \infty = \pm \infty$ (see Arithmetic operations). So, using your notations, for every $t$, you would have [CAUTION ADVISED] $\delta(t+\infty)=\delta(\infty)$, to which I am not able to give a meaning, else than something being a constant evererywhere equal to $\delta(\infty)$. Of course, $\delta$ is a distribution, and should not be treated as a function, but defined as a functional operator, or via limits of functions.

However, so far, I have never seen

a (serious and useful) use of the extended real system in signal processing,
a reference trying to define such a Dirac at infinity.

For the latter, I fear (just intuition) that constructing this with test functions can be troublesome, since two different limits are required (location and amplitude). This is just an analogy, but $\delta(0)\times \delta(0)$ is, as far as I know, not defined.

The Dirac delta generalized function is only defined under an integral

$$
Dirac(t-a) = \int_{-\infty}^\infty f(t)\delta(t-a)dt = f(a)\,
$$

Typically there's also a function under the integral.

Hence

$$
Dirac(t-\infty) = f(-\infty)\,
$$

if one considers $-\infty$ to be a point - not to mention problems evaluating the $Dirac$ at -$\infty$.


Finding the fundamental frequency of a periodic signal Suppose we have the signal $$x(t) = e^{j\omega_1 t} + e^{j\omega_2 t} + e^{j\omega_3 t},$$ where all the frequencies are rationally related (that is, the ratio of any pair of frequencies is a rational number).

How do I prove the fact that the fundamental frequency of the expression above is given by $\omega_0 = \gcd(\omega_1, \omega_2, \omega_3)$?

My attempt is as follows: assume that the signal has a fundamental period $T_0$ which is related to the fundamental frequency by $\omega_0 = 2\pi/T_0$. Now for the signal to be periodic with this period as the fundamental period the following equations must be satisfied $$\omega_1 T_0 = 2 \pi l$$ $$\omega_2 T_0 = 2\pi m$$ $$\omega_3 T_0 = 2\pi n,$$ where $l, m, n \in \mathbb{Z}$ and their values are such that $T_0$ is minimized.

Where do I go from here?
How do I show that the fundamental frequency $\omega_0 = 2\pi/T_0$ is the greatest common divisor of $\omega_1$, $\omega_2$, and $\omega_3$?

1) You are almost there if you have found an $ l,m,n$ and a $T_0$ that fit your equations.  Let $p$ be the largest common factor of $ l,m,n$.  If $p=1$ you are done, otherwise divide $ l,m,n$ and $T_0$ by $p$ and now you are done.

2) From your own equations:

$$ \omega_1 = \omega_0 \cdot l $$
$$ \omega_2 = \omega_0 \cdot m $$
$$ \omega_3 = \omega_0 \cdot n $$

You already know that $l,m,n$ don't have a common factor greater than one.

Hope this helps.

Ced

=================================================

Followup:

Taking it from the top.  You are supposing that you know the frequency ratios:

$$ \frac{\omega_1}{\omega_2} = \frac{a}{b} $$
$$ \frac{\omega_2}{\omega_3} = \frac{c}{d} $$

So $a,b,c,d$ are known.  Therefore:

$$ \frac{l}{m} = \frac{a}{b} $$
$$ \frac{m}{n} = \frac{c}{d} $$

Cross multiply:

$$ lb = ma $$
$$ md = nc $$

Let $ m = bc $ and therefore $ l = ac $ and $ n = bd $

You now have an initial $l,m,n$ set that you can possibly reduce as stated above.  Then the fundamental frequency, $ \omega_0 $, can be solved for using any of the three values.

$$ \omega_0\ = \omega_1 / l = \omega_2 / m = \omega_3 / n $$

Finally, $T_0 = 2\pi / \omega_0 $ gives you the period.

Quite Easily Done

Note, it is possible that none of $l,m,n$ are equal to one.

you show that $$x(t+T_0) = x(t)$$ for all $-\infty<t<+\infty$ and that $T_0$ is the smallest such period that this is the case.


Created Matlab echo signal: Is it an actual "echo?" Or am I just manipulating the signal to make it seem like an echo? This is just an outline of what I did (assume that all my dimensions fit):

    (code that imports audio to read audio as function x)
    X=fft(x);
    (code to define time t and frequency f)

    H=L+M*exp(-1i*2*pi*f*N)+O*exp(-1i*2*pi*f*P);
    % This is where it gets tricky
    % In the frequency domain I have the function above where f is the frequency
    % variable and L M N O P are arbitrary values

    Y=X.*H;
    % again assume dimensions fit

    y=real(ifftshift(ifft(Y)*length(t)));
    %Taking the inverse fft and fftshift was the only way I could get back a signal
    % Once again assume the dimensions fit

Here L M N O P = [.9 .75 1.0 .5 .5]

Is this legit? or am I just pulling my own leg?
Is it an actual "echo?" Or am I just manipulating the signal to make it seem like an echo?

To generate a repetitive echo all you need to do is pass your signal through an appropriate AR filter with a delay corresponding to the echo period. For instance the following code generates (for a damped 1 Hz sinusoid sampled at 100 Hz) an echo which repeats every 1000 samples. The strength of the echo can be adjusted by the AR filter coefficient (0.5 here)

Example Matlab code

% generate 10000 samples of a damped sinusoid at 100 Hz
sig = sin(2*pi*(0:9999)/100).*exp(-(0:9999)/200);

% AR filter to generate a repetitive echo, period 1000 samples (10 seconds)
b = 1;
a= [1 zeros(1,999), 0.5];

echo = filter(b,a,sig);

You may want to experiment with your AR filter to match your echo decay. Also look into 'gapped deconvolution'

Addendum to my earlier post - the AR model for echo generation leads to an infinite number of echos (most of them damped out after a while depending upon the filter coefficient). To generate a finite number of echos use the MA approximation to the AR model, which is:

$$H(z) = \frac{1}{1 - az^{-N}}$$

where $|a| < 1.0$, and $N$ is the reverberant delay in samples.

The MA approximation (series expansion) of the AR transfer function $H(z)$ is

$$H_{MA}(z) = 1 + az^{-N} + a^2z^{-2N} + \ldots$$

So, for instance, to generate just one echo you may want to just filter your signal through an MA filter of the form

$$H_{MA}(z) = 1 + az^{-N}$$

To generate two unequally spaced echoes with differing intensities you could, for example, use an MA filter of the form

$$H_{MA,2}(z) = 1 + az^{-N} + bz^{-M}$$

The Matlab code below generates (a) one echo and (b) two unequally spaced echoes

% generate 10000 samples of a damped sinusoid at 100 Hz
sig = sin(2*pi*(0:9999)/100).*exp(-(0:9999)/200);

% MA filter to generate an echo beginning at 1000 samples (10 seconds)
a = 1;
b= [1 zeros(1,999), 0.5];

% generate echo
reverb1 = filter(b,a,sig);

% MA filter to generate two echos beginning at samples 1000 and 3000
a = 1;
b= [1 zeros(1,999), 0.5 zeros(1,1999) 0.25];

% generate signal with 2 echoes
reverb2 = filter(b,a,sig);

The figure below shows the original signal and the signal with echoes suprimposed from the Matlab output


Camera Calibration using Single Input Image Given a camera, permanently located at the same position, angle and distance from what is being captured, and a single calibration image ( a row of vertical lines [ | | | | ] ), how would one go about using the image to apply correction (barrel distortion) to all subsequent images?

Numerous examples exist online for the camera calibration problem, but most attempt to achieve the correction through guess-work, or highly involved processes to correct real-world images. However, I cannot find anything on calibrating a camera whose position and subject is always at the same position, and where the expected positions of vertical lines are known prior to correction.

The closest I've come to my approach is the plumb-line calibration method described by brown (1971).

Are there any algorithms or pseudo code out there which describes an approach to calibrate using a single input image of vertical lines (where distortion would be clealry visible by the gap between each line)?

Close-Range Camera Calibration - D.C, Brown

Take a look at this paper here: "Straight lines have to be straight" by Faugeras et al.
http://link.springer.com/article/10.1007/PL00013269#page-1

It is straight forward to implement, but essentially removes distortions by observing the straightness of lines.

It should be straight forward to implement if you remove some automation (e.g. you click on the lines manually rather than relying on image processing).

You can take a look at this article http://home.isr.uc.pt/~jpbar/Publication_Source/meloICCV2013.pdf. It calibrated the image from natural image lines. You can recover most of the calibration parameters and then apply them to the subsequent images.


Signification of alternating signs of imaginary and real part of an FFT Let me start by saying that this question is motivated by data analysis that I am doing. It's part of a much larger set of signal processing operations, but in order to keep the question contained I'll ask a general question and hope, someone will be able to help me:

Is there any physical or mathematical characteristic that is known for a signal (containing only the odd Harmonics) whose FFT-imaginary and FFT-real parts have alternating signs ?

As a matter of fact, this alternation happens when I shifted my signal in the time domain. This signal is measured and contains many periods. As I don't know the phase shift between the excitation signal and the measured signal, I am not really sure where the correct signal begins because I don't have the real phase shift. So when I shifted the time signal with (3/4)*Number_of_bins_per_Period, I got a plausible time signal (that suits other operations I am doing on its basis). So I wanted to know if I can find the reason why this shift : (3/4)*Number_of_bins_per_Period, is the correct shift by analyzing its effect on the real and imaginary part of the signal. Before the shifting, the FFT of the signal has these real and imaginary part:
EDIT Please note that I am not asking about the fourier Transformation or why I am getting this behavior. This is clear. As I said, the shift I am doing (3 *Tr/4) I did it intuitively and it gave me afterwards the results I was looking for. I was hoping to find the answer, as WHY exactly this shift is the right one. And I thought maybe the hint would be the alternating signs. So why should be the sogns of the real and imaginary part alternating so that my shift is right . Is this a sort of finger-printing for a specific physical phenomena ?

This is exactly the expected result for a time delay: A delay in time is a linear phase in frequency; so you shifted the amount that would cause 180 degree phase shifts between your frequencies (Note that at your first frequency sample location you have a positive real and positive imaginary, which rotates to positive real and negative imaginary which is a -90° rotation, then follow the same process for your second frequency and you will see it is a +90° rotation; the phase change between the frequencies given is 180°). It is repeating at 180° due to the modulo property of phase, but you are really seeing the evidence of a linear increase in phase that is rotating past 180 degrees at each of your frequency samples where energy exists.

Consider a signal at given frequency going through a fixed length of cable; it will be delayed a certain amount in time, and given the frequency that time will relate to a certain phase shift between the signal at the input and output of that cable. Now consider a frequency that is doubled, the time delay is the same but there will be twice the phase shift between the signal at the input and output in this case.  This is the result of the Fourier Transform of a delay, the magnitude is one but the phase in linearly increasing at a negative rate.

So in your case we are seeing the resulting rotation of 180 degrees for the same distance between your frequencies tones (evidence of a linearly increasing phase).

If you plotted each frequency as a complex vector on an I-Q (real vs imaginary) diagram, as you change the delay of your signal in time, you will see the vector rotate accordingly with no change in magnitude. The higher frequencies will rotate proportionally faster.

Here is another view in case that was not clear describing the same effect and very specific to your case.  The top plot shows sine-waves in time for a first and third harmonic, and the delay at a fixed time offset such that the first harmonic has shifted -90° and the third harmonic as expected has shifted three times as much or -270°.

These last plots demonstrate your case exactly and shows the same thing on a complex plane which is consistent with your complex frequency outputs that have real and imaginary components. Note that I use the common notation I and Q for real and imaginary where "I" is "In Phase" and "Q" is "Quadrature Phase".

Here we see the same thing in that each frequency tone is a vector that rotates around the complex plane at a fixed rate; so the third harmonic is rotating three times as fast as the first harmonic (and the fifth is rotating five times as fast etc). After a given time delay, we see the diagram where the first harmonic has rotated 90° and the third has rotated 270°. To try to keep this clear, I first show the case where the first and third harmonic are in phase, starting at an angle of +45° (Such that I and Q are both positive as in your plots) and we see the exact condition where after a fixed delay the first harmonic has rotated -90° such that I is postive and Q is negative, while the third harmonic in that same time duration has I negative and Q positive. (-270° rotation.)

If we instead have a starting condition identical to your plots, where the first harmonic has both I and Q positive, but the third harmonic is 180° out of phase such that both I and Q are negative, we end up with the diagram below. We see the same thing where the first harmonic shifts -90° and the third shifts -270° after a given time delay. This is exactly the result from your plots: The first harmonic goes from I+, Q+ to I+, Q-, while the third harmonic goes from I-,Q- to I+,Q-.

Note that applying a finite length FFT to a longer length arbitrary frequency sinusoid applies a rectangular window (e.g. cuts off the stuff before and after the finite FFT input vector) to that sinusoid.  A rectangular window in the time domain is the same as convolution with a Sinc function in the frequency domain.  Note that a Sinc function has alternating lobes (above and below zero).  So if the frequency of a sinusoid is shifted off a bin center (e.g. is not integer periodic within the FFT width), the Sinc function's alternating lobes will be sampled by the FFT result, thus resulting in alternating signs in the real and imaginary components of the FFT result.

(For exactly integer periodic in aperture signals, the Sinc is sampled mostly at zero crossings, so becomes invisible in the FFT result.)

Rotating a signal in the time domain sample array (equivalent to shifting a sufficiently zero-padded signal in the time domain) results in twisting the phase of the complex FFT result vector, with the amount of twist in the frequency domain proportional to the amount of shift the time domain.


If humans can only hear up to 20 kHz frequency sound, why is music audio sampled at 44.1 kHz? I read in some places that music is mostly sampled at 44.1 kHz whereas we can only hear up to 20 kHz. Why is it?

The sampling rate of a real signal needs to be greater than twice the signal bandwidth. Audio practically starts at 0 Hz, so the highest frequency present in audio recorded at 44.1 kHz is 22.05 kHz (22.05 kHz bandwidth).
Perfect brickwall filters are mathematically impossible, so we can't just perfectly cut off frequencies above 20 kHz.  The extra 2 kHz is for the roll-off of the filters; it's "wiggle room" in which the audio can alias due to imperfect filters, but we can't hear it.
The specific value of 44.1 kHz was compatible with both PAL and NTSC video frame rates used at the time.

Note that the rationale is published in many places: Wikipedia: Why 44.1 kHz?

In order to faithfully reproduce a signal, the faster the sample rate the better. ~40 kHz was chosen, because it was a low sample rate that most people can't tell the difference for (when reconstructed). When audio sampling was introduced, memory and storage was expensive and higher sample rates were not cheaply possible.

At double the upper limit of human hearing two samples per cycle is very poor reconstruction, even if it meets the Nyquist criteria for sampling signals, a simple chart depicting a sine wave with two samples per cycle will show you how poor two samples per cycle is in reproducing a waveform. You can literally turn a sine wave into square wave; it is a good thing at 20 kHz nobody can tell. I bet a dog could though.


Looking for an algorithm that computes "Similarity" between two signals I'm new in this forum! ... and I'm not an expert in signal processing! ... So I hope to find help here:-)
I was wondering if someone could point me to an algorithm/technique that is used to compare time dependent signals. Ideally, this hypothetical algorithm would take in 2 signals as inputs and return a number that would be the percentage similarity between the signals.
In other terms I’m looking for a technique that allows to establish a degree of “similarity” between the two signals even if there is a delay or a shift between them. Obviously, if this technique could be implemented in a recursive way or with a small number of point ...it would be great! :-)

Thanks

Coherence? As in: http://en.wikipedia.org/wiki/Coherence_(signal_processing). By comparing the signal spectrums, you get real number numbers between 0 and 1, which each number is a measure of the similarity of the two signals at each frequency.  But it might be difficult for a beginner to understand and program.  You might look for simpler explanations and/or examples by doing a Google search for 'signal coherence', which results in about 17 million hits.

something like  sum(x(n)*y(n))/sqrt(sum(x(n)^2)*sum(y(n)^2)).
this would give a 1 when x(n)=y(n)...


Basic question regarding third octave filters I would like to use third octave filters on a time signal.
Do I have to apply all third octave filters to a single frequency or do I only need to apply the one that has its center frequency the closest to the frequency that I am looking at?

I would like to use a simple example to make my quesion clear:
Let's say I have a basic sin signal: signal = sin(2*pi*f*t) with f = 57Hz.
57Hz is almost right between the center frequencies of 50Hz and 63Hz, but the closest would be 63Hz. Do I apply all filters or only the one with a center frequency of 63Hz?

I am pretty sure that the answer is that I will always have to apply all filters in theory. Of course filters with much higher center frequencies will result in insignificant amplitudes...

Someone released my natural inheritences and regulatory systems calculated as it formed my employment, rewards  and inheritences payable to me obly all as one from more than just my g.c., Tricon, Albertsons, Tetra,  and more...as someone or a group if someones have changed my statuses to a protective state and have done so as if court ordered. I gave nevet in my life been so isolated during major life upsets such as missing children and death with zero compliance of the government and of others involved. My personally named for claim and my personal poatal mail and hijacked emails go uncorresponded and still remail displaced like my paymemt and documentations I have  applied for in the last entire seven years and my lifes worth is wromgly dispensed. From my insurance policies to my ach as named, global this state is acting as if I am dead. Even as the shut down goes, these are being processed illegally like my procured custody reformed payments from my base assigned living EXPOSURE money.
A step father and father with uncle's and aunts steal. Like they always have.
I am Lindsey ann maxfields
717 briarwood lane Hurst rd.


Split and Merge image segmentation algorithm in OpenCV I asked google and I've tried to find something in the API of OpenCV, but hadn't any good hints on split and merge image segmentation method.
I want to test it, if it's working better than my proposed method which I posted a while ago here on SP.
Do you know if there exists a free library or OpenCV extension for this segmentation method?
Any good hints for writing it fast by myself are welcome, if there is no available solution.

kind regards

I asked Google again for you, but I did manage to find some hits in the end. There is already a very good question on stackoverflow concerning the exact same thing you are interested in.

There is a very nice explanation of split-and-merge provided in one of the answers, as well as simplified pseudocode.

The other answer provides a link to the implementation (although I did not check it, I'm not sure which language the implementation is in).

If you want to use split-and-merge as a reference segmentation technique to compare your idea with, other segmentation methods offered by OpenCV might also be interesting:

functions useful in image segmentation from the OpenCV manual
a chapter dealing with image segmentation in OpenCV from Gary Bradski; Adrian Kaehler: Learning OpenCV
article about Pyramid segmentation
Mean-shift segmentation: stackoverflow answer, scientific article
Watershed segmentation seems fairly well explained by Wikipedia

In the end, you can choose one of the already offered methods as the reference method after all, or, if you want to present amazing, all-encompassing comparisson, you might also want to use more than one reference method

I stumbled upon this link. It offers an implementation of watershed segmentation using opencv. There is also an example of using grabcut, which can help in segmenting out fo


Bandpass Filter Order vs Sample Numbers of Noisy Signal I'm working on signal filtering in MATLAB. i wrote sample code that generate 50 Hz signal with some noisy frequency and use Kaiser bandpass filter to get only 50Hz frequency. this code work great! but there is a big problem!

Sample Rate is 2000 Hz and Filter Order is about 39000 , and my signal length is 50000 samples! in other hand i have to wait about 25 seconds to gathering my signal! that is too long time!, if i decrease signal length ( to about 6000 samples ) filter result is so noisy because of filter order. if i decrease filter order , filter signal result is too noisy too! the best time that i cloud wait for gathering data is about 3 seconds.

is there any way to get same result with 6000 samples of signal and high order filter? ( 39000 )

thanks.

Fs = 2000;                    % Sampling frequency
T = 1/Fs;                     % Sample time
L = 50000;                     % Length of signal
t = (0:L-1)*T;                % Time vector

Hd =  KAISERFilter;

for i=1:50

a = 5;
b = 500;
r1 = a + (b-a).*rand(1,1);
r2 = a + (b-a).*rand(1,1);
r3 = a + (b-a).*rand(1,1);
r4 = a + (b-a).*rand(1,1);

n = r1 *sin(2*pi*49.2*t) + r2 *sin(2*pi*50.5*t) + 0.1 * randn(size(t)) +  r3 *sin(2*pi*15.5*t) +  r4 *sin(2*pi*110*t) ;
x = 30 *sin(2*pi*50*t); % main frequency
y = x+n;

subplot(2,1,1);
plot(y)
title('Signal')
xlabel('time (milliseconds)');

yf = filter(Hd.Numerator,1,y);

and Kaiser filter code is :

function Hd = KieserFilter

% All frequency values are in Hz.
Fs = 2000;  % Sampling Frequency

Fstop1 = 49.5;              % First Stopband Frequency
Fpass1 = 49.9;            % First Passband Frequency
Fpass2 = 50.1;            % Second Passband Frequency
Fstop2 = 50.5;              % Second Stopband Frequency
Dstop1 = 1e-06;           % First Stopband Attenuation
Dpass  = 0.057501127785;  % Passband Ripple
Dstop2 = 1e-06;           % Second Stopband Attenuation
flag   = 'scale';         % Sampling Flag

% Calculate the order from the parameters using KAISERORD.
[N,Wn,BETA,TYPE] = kaiserord([Fstop1 Fpass1 Fpass2 Fstop2]/(Fs/2), [0 ...
                             1 0], [Dstop1 Dpass Dstop2]);

% Calculate the coefficients using the FIR1 function.
b  = fir1(N, Wn, TYPE, kaiser(N+1, BETA), flag);
Hd = dfilt.dffir(b);

You are dealing with a time frequency uncertainty in estimation problem.  It takes longer than 3 seconds for there to be enough information in a signal to tell the difference between 3 clearly separated sinewaves between 49 and 51 Hz, and only 1 fat (slightly modulated) one.  And you can't filter out a signal unless there is enough information to tell whether it's even there or not.

Here is some tutorials which you may need to refer.  Zero phase filtering adjust the time delay produced by digital filters.

zero phase filtering to remove noises from signal

50 Hz and 60 Hz power hum removal from signals.

https://youtu.be/iwR1QIspjM8

high frequency noise removal tutorials in matlab


Orthogonal signals in frequency and time domains I have 3 orthogonal signals obtained by Fourier transformation. But when I do the inverse Fourier transform in Matlab, my signals are no longer orthogonal in the time domain.

Here is my question:

How to keep the orthogonality of signals in the time domain?

Preamble: a scalar product is a positive semi-definite bilinear form, and we will suppose here we want the standard definitions.

Using their Fourier domain is (one of) the simplest way to construct orthogonal signals. It suffices to choose signals with disjoint frequency supports. Note that is not a necessary condition: some signals can be orthogonal with the same frequency supports, like $\cos$ and $\sin$, and more generally Hilbert transform pairs.

This is due to the Parseval-Plancherel theorem (inner product and norm preservation under Fourier transform): with disjoint frequency support, the scalar product is zero, and hence the vectors are orthogonal.

So in the following results, we take a random spectrum for a real signal. We then segment it into three subsets of points, mutually exclusive, and whose union give the full frequency axis. And then do an inverse FFT to obtain three subsignals whose inner products are:

S1.S2 = -4.7705e-18
S2.S3 = 1.7564e-17
S3.S1 = -1.2645e-17

The corresponding signals and spectra are displayed below, and the code follows (in Matlab):

function seDsp51117(nSample)

if nargin < 1
    nSample = 512;
end
nSampleHalfSpectrum = nSample/2+1;
halfSpectrum = randn(nSampleHalfSpectrum,1)+1i*randn(nSampleHalfSpectrum,1);
halfSpectrum(1) = real(halfSpectrum(1));
halfSpectrum(end) = real(halfSpectrum(end));
nSampleSpectrum1 = ceil(nSampleHalfSpectrum*rand(1));
nSampleSpectrum2 = ceil((nSampleHalfSpectrum-nSampleSpectrum1)*rand(1));
nSampleSpectrum3 = nSampleHalfSpectrum - nSampleSpectrum1 - nSampleSpectrum2;

idxSpectrumRand   = randperm(nSampleHalfSpectrum);
bndSpectrumRand = round(linspace(1,nSampleHalfSpectrum+1,4));

idxSpectrum1 = idxSpectrumRand(bndSpectrumRand(1):bndSpectrumRand(2)-1);
idxSpectrum2 = idxSpectrumRand(bndSpectrumRand(2):bndSpectrumRand(3)-1);
idxSpectrum3 = idxSpectrumRand(bndSpectrumRand(3):bndSpectrumRand(4)-1);

[signal,halfSpectrumSelect] = IFFTReal(halfSpectrum);
[signal1,halfSpectrumSelect1] = IFFTReal(halfSpectrum,idxSpectrum1);
[signal2,halfSpectrumSelect2] = IFFTReal(halfSpectrum,idxSpectrum2);
[signal3,halfSpectrumSelect3] = IFFTReal(halfSpectrum,idxSpectrum3);
figure(1);clf
subplot(4,2,1)
plot((signal));axis tight;grid on
xlabel('Original signal')
subplot(4,2,3)
plot((signal1));axis tight;grid on
xlabel('Subsignal S1')
subplot(4,2,5)
plot((signal2));axis tight;grid on
xlabel('Subsignal S2')
subplot(4,2,7)
plot((signal3));axis tight;grid on
xlabel('Subsignal S3')

subplot(4,2,2);hold on
plot(abs([halfSpectrumSelect]),'o');
plot(abs([halfSpectrumSelect]),'x');
xlabel('Original spectrum')
axis tight;grid on
subplot(4,2,4);hold on
plot(abs([halfSpectrumSelect]),'o');
plot(abs([halfSpectrumSelect1]),'x');
axis tight;grid on
xlabel('Subsignal spectrum for S1')

subplot(4,2,6);hold on
plot(abs([halfSpectrumSelect]),'o');
plot(abs([halfSpectrumSelect2]),'x');
axis tight;grid on
xlabel('Subsignal spectrum for S2')
subplot(4,2,8);hold on
plot(abs([halfSpectrumSelect]),'o');
plot(abs([halfSpectrumSelect3]),'x');
axis tight;grid on
xlabel('Subsignal spectrum for S3')

sp12 = dot(signal1,signal2);
sp23 = dot(signal2,signal3);
sp31 = dot(signal3,signal1);

disp(['S1.S2 = ',num2str(sp12)]);
disp(['S2.S3 = ',num2str(sp23)]);
disp(['S3.S1 = ',num2str(sp31)]);

function [signal,halfSpectrumSelect] = IFFTReal(halfSpectrum,idxSpectrum)
if nargin < 2
    idxSpectrum = 1:length(halfSpectrum);
end
halfSpectrumSelect = zeros(size(halfSpectrum));
halfSpectrumSelect(idxSpectrum) = halfSpectrum(idxSpectrum);

%% Hermitian glue
wholeSpectrum = [halfSpectrumSelect;conj(flipud(halfSpectrumSelect(2:end-1)))];
ratioImagReal = 10^(-6);
tmp = ifft(wholeSpectrum);
normRatio =  norm(imag(tmp),1)/norm(real(tmp),1);
    signal = real(tmp);
if normRatio > ratioImagReal
warning('Check the imaginery value');
    end

There is nowhere it will occur if the signal is orthogonal in freqecey domain then it is orthogonal in time domain.the Fourier transform is just a tool to convert one to another it will not change signal properties(maybe even and odd), the orthogonality in time domain is defined between two signals f(x) & g(x) as$ \int( f(x).g(x)dx) =0$ for one period if it is periodic or the entire signal, hope it is helpful


Combined signal folding and shifting Let $g[n] = f[-n]$(n is the time index). Then what will be expression for $g[n+k]$.

Is $g[n+k] = f[-n + k]$ or $g[n+k] = f[-n - k]$. Which one is true and Why?

The statement $g[n]=f[-n]$ is meaningless unless you assign some numerical value to $n$, e.g. $g[3]=f[-3]$, and if you do say that $n=3$, then the equality
$g[3] = f[-3]$ tells you nothing about the relationship between $g[n+k]=g[3+k]$
to the value taken on $f[\cdot]$ for any choice of argument.

What you really meant to say is that

  $$\text{For all integers } n, g[n]=f[-n]$$

which asserts that the equality holds for all choices of $n$, and in
this case, it follows by definition that $g[n+k] = f[-n-k]$;
you just need to negate the argument $n+k$ of $g$ (to get $-n-k$)
and use it (it meaning $-n-k$) as the argument of $f$ exactly as
the definition says you should.  The
complicated arguments in another answer that claim that $g[n+k]=f[-n+k]$
are wrong, as discussed in the comments following that answer.

The answer is f[-n+k].

I'll try to explain this with example. Let the signals be discrete & k=2. The bold numbers indicate sample at time 0(I am not well versed with editing here)

Let f[n] = {1,1,1,1}
Thus f[-n] = {1,1,1,1} = g[n]
Now g[n+2] = {1,1,1,1,0,0}

Now how can you obtain the same answer from f[n]
f[n-2] = {0,0,1,1,1,1}
And if we invert this signal ie f[-(n-2)], we get {1,1,1,1,0,0} which is g[n+2].
Thus g[n+2] = f[-(n-2)] = f[-n+2]
Hence proved.
The theoretical answer are the rules followed in Shift & Scale
y[n] = x[an+b]
If you shift first & then scale, then
1) y[n] = x[n+k]
2) Then scale y[n] by a

If you scale first & then shift, then
1) y[n] = x[a(n+(b/a))] ie scale x[n] by a first
2) Shift y[n] by (b/a)


Intuitive explanation of cross-correlation in frequency domain According to the cross-correlation theorem : the cross-correlation between two signals is equal to the product of fourier transform of one signal multiplied by complex conjugate of fourier transform of another signal. After doing this, when we take the ifft of the product signal, we get a peak which indicates the shift between two signals.

I am not able to understand how this works? Why would i get a peak which indicates the shift between two signals. I got the math from : http://mathworld.wolfram.com/Cross-CorrelationTheorem.html
but i am not able to understand as to what this means intuitatively. Can somebody please provide some explanation or point me to the right documents?

Thanks!

The concept is based on the convolution theorem, which states that for two signals $x(t)$ and $y(t)$, the product of their Fourier transforms $X(f)$ and $Y(f)$ is equal to the Fourier transform of the convolution of the two signals. That is:

$$
\mathcal{F}\{x(t) * y(t)\} = \mathcal{F}\{x(t)\}\mathcal{F}\{y(t)\}
$$

You can read more on the derivation of this theorem at the above Wikipedia link. Now, convolution is a very important operation for linear systems in itself, so the theory on its properties is well-developed.

However, what you're looking for is the cross-correlation between $x(t)$ and $y(t)$. Here's the key: you the cross-correlation integral is equivalent to the convolution integral if one of the input signals is conjugated and time-reversed. This allows you to utilize theory developed for evaluating convolutions (like frequency-domain techniques for calculating them quickly) and apply them to correlations.

In your example, you're calculating the following:

$$
\mathcal{F}\{x(t)\}\left(\mathcal{F}\{y(t)\}\right)^*
$$

Recall that in the Fourier domain, complex conjugation is equivalent to time reversal in the time domain (this follows directly from the definition of the Fourier transform). Therefore, using the first equation given above, we can state that:

$$
\mathcal{F}\{x(t) * y^*(-t)\} = \mathcal{F}\{x(t)\}\left(\mathcal{F}\{y(t)\}\right)^*
$$

If you then take the inverse Fourier transform of this equation, the signal you're left with is the cross-correlation between $x(t)$ and $y(t)$.

% Matlab function for frequency domain cross correlation
function [Lag,C]=xcorrf(X,Y,L)
% X, Y ---> Input vectors
% L --->  maximum lag (must be less than minimum of (length of X, Y)
% C ---> correlation vector
% Lag ---> lag times
X=X(:);
Y=Y(:);
s1=size(X);
s2=size(Y);
D=min(s1(1,1),s2(1,1));
for i=1:L
    X1=ifft(fft(X(1:D-i,:)).*conj(fft(Y(i+1:D,1))));
    C(i,1)=X1(1,1);
end

C=flipud(C);
X1=ifft(fft(X(1:D,:)).*conj(fft(Y(1:D,1))));
C(L+1,1)=X1(1,1);
for i=1:L
    X1=ifft(fft(Y(1:D-i,:)).*conj(fft(X(i+1:D,1))));
    C(i+L+1,1)=X1(1,1);
end
Lag=-L:1:L;
end


The impact on the frequency of adding zero samples and non zero samples Can someone give me a brief explanation of what is going on in here? I struggle to understand the full differences between the impact on the frequency of added zero samples and adding samples which are not zeroes. (In between the samples).
I learnt this course 100% from the lectures so if I have some misconception regarding the basics, I'd be glad if you could kindly help me understand.
The key points I know so far:

Adding zero(ed) samples in between the existing ones in the time domain, squizes the frequency domain and adds height to it.
The squized frequency domain is squized while it's normalised (once we turn it to be periodic with period of 2π).
Adding non zero(ed) samples to our time domain, we get an added zeroes vector to the end of the frequency domain's period. (Padding zeroes in frequency). Which means that when it's being normalized into 2π period's length, then it's also being squized.

From that we conclude, that adding zero samples and non zero samples do the same exact thing.
[Question:] where is my mistake in the above logic sequence? Where is the difference between the impact of adding zero samples in the time domain on the frequency vs adding non zero samples in the time domain on the frequency?

As soon as you sample in the nuiqist rate, sampling faster won’t add information(!). Sampling zeroes can be looked at as approximation mistake(!) It gives the same frequency as the original signal mathematically because it’s multiplying by zero and IF by contradiction, there would have been any other value instead of the original value, which WOULD [again by contradiction] add some information aside from what we have with the niquist samples, THEN we would need a higher niquist rate and that would be a contradiction to our given Max. Niquist frequency.
Yes, indeed sampled zeroes gives the same signal and they are considered as a measure mistake given the Niquist frequency!


The reasons for filter explosion Assume signal $x$, sampled at $f_s = 44100\; \mathrm{Hz}$. I tried to filter it using the Butterworth bandpass filter ( $30\; \mathrm{Hz} - 70\; \mathrm{Hz}$) of order $8$. However, as a result I get a vector with most elements being NaN (and some of them extremely small, approx.  $-2.5 \cdot 10^{306}$`).

If I try the same filter of order $6$, I get results as expected. What could be possible reason for order $8$ filter to 'explode'?

Here is the MATLAB code, just in case I made an error which I don't see:

[b, a] = butter(4, [60 / fs, 140/fs]);
x_filtered = filter(b, a, x);

Basically you never want to use the Transfer Function representation (with b and a) and rather use the Zeros-Poles-Gain (z,p,k). This will allow you to avoid the numerical errors. In your case you might design your filter in following way:

fs = 44100;           % Sampling frequency
Wp = [30, 70]/(fs/2); % Pass band frequencies (as normalized frequency)
Ws = [20, 90]/(fs/2); % Stop band frequencies
Rp = 3;               % Ripple at pass band
Rs = 50;              % Ripple at stop band

[n, Wn] = buttord(Wp, Ws, Rp, Rs);     % Get order and omega vector
[z, p, k] = butter(n, Wn, 'bandpass'); % Design filter accordingly
[sos, g] = zp2sos(z, p, k);            % Convert to state matrix
Hd = dfilt.df2sos(sos, g);             % Create the filter object

Which for some dummy random signal:

x = rand(1,100000);
y = filter(Hd, x);

Will produce a stable output:

And here is the filter frequency response (everything looks as requested):

Read the limitations section  of the documentation of the butter function. http://www.mathworks.com/help/signal/ref/butter.html#bt5pm8a


Is there a difference in bicubic interpolation of images vs patches? I need to extract low resolution patches from an image. For example, 32x32 patches from a 500x500 image.

Method One

Down res 500x500 image (by bicubic interpolation, both down res and up res).
Extract 32x32 patches.

Method Two

Extract 32x32 patches.
Down res (by bicubic interpolation, both down res and up res) 32x32 patches.

Will the patches extracted be identical? If they aren't, how different will they be?

The short answer is no. The order of operations doesn't matter and the final result would be similar within reasonable approximation errors.

The long answer is:

A given image $I$ has a pixel resolution of $M \times N$ pixels and a pixel size of $P \times Q$mm. Therefore, the image physical size is $M \cdot P \times N \cdot Q$mm.

If you downsample that (with any interpolation method) by some factor $d \in \left(0 \ldots 1\right)$, then its new physical size will be $d \cdot M \cdot P \times d \cdot N \cdot Q$. Obviously, here $d \cdot M \le 32$ and $d \cdot N \le 32$.

By the way, since at this point we are talking about reduction in physical size, this means that some of the original pixels, some of the original actual information content, had to be thrown away

If you then select the $32 \times 32$ patch, you are selecting some "space" out of the "shrunk" image, some of which information content has now been thrown away. ($32 \cdot d \cdot M \cdot P$)

If you were to select a $32 \times 32$ patch out of the original $M \cdot P \times N \cdot Q$ image, then you select a patch at full resolution.

Once you have selected the patch, the rest of the image is "thrown away". Any downsampling of that patch will produce a new patch whose dimensions (and actual information content) will be lower than $32 \times 32$. If you are thinking at this point "Oh, its alright, I will cut out a $32 \times 32$ patch, downsample it to (for instance) $12 \times 12$ and then upsample that patch back to $32 \times 32$" then, in terms of actual information content you are not gaining anything and you are back to the first case.

In fact, practically, because of the double interpolation (one to downsample and then one to upsample), this second case will have more opportunity for numerical errors but in terms of actual information content, actual spatial resolution, you don't gain anything.

The "dangerous" bit here is the downsampling. Once that takes place, that is where actual information is thrown away. This is why, when you are trying to resample a signal by a non-integer factor $\frac{u}{v}$, you first upsample by $u$ (where no new information content is added but the same information is distributed over larger space) and then downsample by $v$ (where actual information get's limited).

In your case, you have two different routes that both include downsampling by the same factor at some point. Their end result will be similar within reasonable approximations afforded by the interpolation steps.

Hope this helps.


Filter Implementation I am rather new to the world of signal processing, and am struggling to understand a fundamental concept: How are filters actually implemented?

I have read a significant portion of this online book, and scrounged the internet, finding snippets of useful information here and there, but I cannot quite tie it all together.

In brief, I have a vibration profile (~200,000 samples) in the time domain, which I would like to analyze using a particular frequency weighting. This involves a four-step filter, using a (Butterworth) high-pass and low-pass, followed by an a-v transition and an upward step. The analog equations are shown below.

Yet I am still struggling on how to implement them. I have seen mentions of the bilinear transform, used in programs in MATLAB/Python, but such implementations seem to omit the 's' or 'p' variable from the filters shown, typically creating the numerator 'B' and the denominator 'A' from the coefficients of the analog filters.

Other sources suggest taking the Fourier of the filter equations and the signal, multiplying, then applying the inverse Fourier to bring it back into the time domain. Yet I cannot discern how to take the Fourier of an equation.

My same problem occurs with using convolution in the time domain, as I have not understood how to convolve an equation yet.

I am sure I have come near the answer, and I know this is a basic principle, but I cannot seem to comprehend it, and my Googling is off the mark. Any help on this would be appreciated. General advice is also welcome, as I have a lot to learn.

** If important, I am currently working in Python, using the scipy.signal package, with access to MATLAB for testing and evaluation purposes.

Currently I am looking at substituting j2πf for p, then simply running the data through the equations and taking the inverse Fourier to bring it back into the time domain.

Alternatively,  my eyes have just noticed what was staring me in the face: the second equation given for each filter is in terms of 'f', which would suggest I am able to simply run the frequency-domain data through this equation before taking the inverse Fourier to recover the time-domain. If this is a correct realization, this question can be closed (or I can 'answer' my own, now self-evident, question).

There are a few different ways to proceed. They involve doing some algebra but it is more an issue of tedium. The notation of your analog filters is a little non standard, essentially your $p$ is $s$ in most books.

Traditionally, phase wasn't considered important in hearing so standards tended to be specified by magnitude, so there is some ambiguity with respect to the digital implementation's phase response. You have 2 basic simple choices, a digital IIR filter or FIR filter. FIR is likely to be easier but easy and better are typically not the same. FIR can be direct or FFT based.  IIR has a number of choices such as biquads or direct. Floating point math has its own easy/better issues.

One path way is to focus on the terms inside the magnitude brackets, i.e. the complex analog transfer functions. The other path is to focus on the analog magnitude.

So, given the analog transfer functions, you can either use a table of Laplace transforms and deduce the impulse response of each component you listed, or use the bilinear transform to express each component that is in $s$ in terms of $z$.

If you derive the impulse responses, you can sample points of the impulse response and provided a generally exponential decay, use those samples as coefficients of a FIR filter of sufficient length. The FIR filters are then applied sequentially. One nice property of ideal LTI filters, the order that you apply the filters are ideally interchangeable.

If you choose the bilinear transform and do the algebra to derive the discrete time numerators and denominators, you need to be aware that the bilinear transform squeezes the response as it approaches the Nyquist frequency. One needs to plot the discrete time response to asses how the corresponding analog response is distorted. If some modest tweaking results in an acceptable response then, the numerators and denominators need to be factored so that they are compatible with the IIR implementation, such as a cascade of biquads.

If you go with the magnitudes, the task becomes a direct sampling from the analog frequency domain to the discrete time periodic frequency domain. One can choose a linear phase response, and appropriately impose symmetry and define the frequency domain FIR filter. This might not sound like an analog implementation.

The essential challenge is mapping a frequency response in $-\infty < \Omega < \infty$ to $-\pi \le \omega \le \pi$

There are other approaches like state variables, so this is not the entire story.

you already have the basic knowledge about analog filter, then get a book of "digital signal processing", check the chapter about "z transform","transfer function",  "digital filter design"


Efficient frequency determination Is there any algorithm that tells you the most dominant frequency component, i.e. does argmax(fft(x)) or similar, that's more computationally efficient than argmax(fft(x))?

I need to continuously find which fft bin contains a tone, of (up to) 4096 bins in the full signal bandwidth. This is part of an RF detection system implemented in FPGA. As the input is taken from the air, the input signal is very noisy. Calculating argmax(fft(x)) with a length 4096 fft does what I need, but I was wondering if there's a way to do it without calculating all the 'unnecessary' fft outputs, to save fpga resources.

For example could there be some fft-type algorithm that bisects the spectrum to find the peak without calculating many of the outputs? Though I suspect this wouldn't really work without high SNR...

Is there any algorithm that tells you the most dominant frequency component, i.e. does argmax(fft(x)) or similar, that's more computationally efficient than argmax(fft(x))?

Mostly no.

The "trouble" is not the argmax, but fft. One alternative to this is to use autocorrelation and then determine whether that function is periodic and where does it crosses zero. The other alternative is to not evaluate the FFT in all 2048 possible points via the use of the Goertzel algorithm.

This covers why the answer is "mostly" no.

But here is a case where the answer might be yes.

Instead of trying to estimate the frequency via the typical "sinusoidal" concept, you can pass your signal through a sign function which will make it look like a square wave and then count its zero-crossings. (For complex $x$, you would have to do it twice). This is an inexpensive way of getting some estimate of "dominant" frequency.

Hope this helps.

imin = max(0, int(np.floor(not_to_fftbin(NOTE_MIN-1))))
imax = min(SAMPLES_PER_FFT, int(np.ceil(note_to_fftbin(NOTE_MAX+1))))

(np.abs(fft[imin:imax]).argmax() + imin)*FREQ_STEP


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));