Can you detect circles only using the hough line transform? I'm not talking about the Hough Circle Transform (where you replace the representing eqns, hence the parameters with a 3D set) but a Hough Line Transform.

My solution to this:

A hough line transform maps lines to a single point in (p,theta) space. This point contains the number of points lying on the line in $(x,y)$ (cartesian) space. Conversely, a point $\mathbf{A}=(x,y) \in \mathbb{R}^2$ maps to a curve in polar $(r,\theta)$ space. This is because a point belongs to an infinite number of lines. Now, imagine a circle:

where the tangent $\mathbf{AB}$ is perpendicular to the line connecting the center and the point of contact ($\mathbf{NA}$). Let's call this line the radial line. Since the circle is a symmetrical object, this applies to every point on it.

This means essentially, that a circle is constructed of such infinite tangent lines which are perpendicular to the radial line. The tangents map to the bin in polar space: $\mathbf{A} \rightarrow (r, \theta)$. So essentially, the circle in $(x,y)$ space becomes a line in $(p,\theta)$ space. Am I correct in saying this?

Then could we apply a Hough transform on the Hough matrix itself to detect circles in the original $(X,Y)$ space?

(I know academics can be ruthless but newbee here so please be kind when tearing me a new one.)

No this won't work. Simply because we do not know the center of the circle in advance (no oracle telling us) and therefore we cannot find a voting space parameterization where $r$ is fixed and only $\theta$ changes. For each point both quantities will change and this will still generate a curve in the Hough space. Having said that, I think below there is something interesting which uses tangent lines to detect circles:

Although what we will be doing is not exactly a Hough line transform, your way of thinking is on the right track. There is a strong connection there but because the tangent is not well represented in the Hough-circle transform, it probably won't work. Yet, what you have described brought me to the following idea:

It is possible to create an accumulator in the cartesian space ($X$-$Y$), where each pixel casts multiple votes as a partial line segment, along its gradient direction (or opposite gradient direction). The center of the circle will lie at the intersection of those segments. The width of the line segment can be related to the radius of the circle we like to detect. I think this approach is tightly coupled to operating in the $\textit{dual space}$ instead of the original and has demonstrated advantages, especially for conic fitting (well, a circle is basically a conic). A dual space is nothing but the space composed of the tangents. Each tangent is a point, just like in Hough.

Out of interest and because I think this is not a bad idea, I implemented this operation in MATLAB. Of course, instead of explicit edge-detection, I would simply use a gradient magnitude thresholding and use the image gradient as a cue of the direction of $\mathbf{NA}$ vector.

Here is the original test image:

And here is the accumulator space I get when I search for circles with $r=30$ pixels:

With peaks detected by the houghpeaks function:

The same peaks visualized back on the image:

It is possible to see that thanks to the voting along the line, the approach is pretty robust to the noise in the radius, i.e. we do not need to specify the radius explicitly.

For interested reader, I also provide the MATLAB code I used to generate those images:

function [ x,y,V ] = hough_circle_dual( I, maxRadius, threshold, numMaxCircles )

[M, N] = size(I);
[Gx, Gy] = imgradientxy(I);
[Gmag, Gdir] = imgradient(Gx, Gy);

V = zeros(M, N);

for i=1:M
    for j=1:N
        if (Gmag(i,j)>threshold)
            g = [Gx(i, j) Gy(i, j)]./Gmag(i,j);
            p2 = maxRadius.*g + [j, i];
            [x, y] = bresenham(j,i,p2(1), p2(2));
            ind = (x>0 & y>0 & x<=N & y<=M);
            x = x(ind); y = y(ind);
            ind = sub2ind([M,N], y, x);
            V(ind) = V(ind)+1;
        end
    end
end
peaks = houghpeaks(V, numMaxCircles);
x = peaks(:,1);
y = peaks(:,2);
end

% generates a bresenham line (compactly written)
function [x, y]=bresenham(x1,y1,x2,y2)
x1=round(x1); x2=round(x2);
y1=round(y1); y2=round(y2);
dx=abs(x2-x1); dy=abs(y2-y1);
steep=abs(dy)>abs(dx);
if steep t=dx;dx=dy;dy=t; end

if dy==0, q=zeros(dx+1,1);
else, q=[0;diff(mod([floor(dx/2):-dy:-dy*dx+floor(dx/2)]',dx))>=0];
end

if steep
    if y1<=y2 y=[y1:y2]'; else y=[y1:-1:y2]'; end
    if x1<=x2 x=x1+cumsum(q);else x=x1-cumsum(q); end
else
    if x1<=x2 x=[x1:x2]'; else x=[x1:-1:x2]'; end
    if y1<=y2 y=y1+cumsum(q);else y=y1-cumsum(q); end
end
end

Here is how I would call this function in MATLAB 2016:

I=imread('coins.png');
[x, y, V] = hough_circle_dual(I, 30, 10, 10);
figure,imshow(V,[]);
hold on, plot(y,x,'r+');

The Bresenham line generation used here is taken from the FileExchange post.


why upsampling contracts the spectrum? I'm trying to understand why Matlab function upsample causes contraction in the spectrum, like the following figure shows.

Also showed in mathwork online document.

enter link description here

Filtering after upsampling removes the imaging artifacts, but contraction still stays, as shown below. Could anyone help explain in details why?

enter link description here

I prefer to think of it as expanding the spectrum due to increasing the sample rate.  When you expand the frequency spectrum you reveal the aliases that were already there in the other Nyquist zones.  See this answer for more detail about this.


How to multiply two FFTs in MATLAB? I want to find the output signal(y) from input signal(x) and impulse response(h). To do that first I take the fft of both x and h, then multiply them, and finally taking the ifft of the result. The input signal is a 5 second sound. I tried to equalize the lengths of the FFTs but it didnt work. How can I overcome this?
[x,Fs]=audioread('hello.wav');

T=1/Fs;

L=220500;%length of the original sound

d=8820;%delay in terms of samples

a=0.8;%attenuation coefficient

h = [1 zeros(1,d-1) a];%impulse response h(n)

X=fft(x);

h1 =[h; zeros(length(x)-length(h),1)];

H=fft(h1);

Y=H.*X;

y=ifft(Y);

You have two signals of length $N = 220500$ and $ L =  8821$ samples long, and you want to obtain their convolution $y[n]$ of length $K = N + L-1 = 229320$ samples long, by using a frequency-domain DFT//FFT method...
Then you have to compute $K = N + L -1 = 229320$ sample long FFTs of both signals $x[n]$ and $h[n]$, and then multiply them, and then invert the resulting spectrum back into time-domain as given by the following MATLAB line.
% assuming x[n] and h[n] exist....
N = 220500;  % length of x[n]
L = 8821;    % length of impulse response h[n]
K = N+ K -1; % length of output y[n] = x[n] * h[n]

y = real( ifft( fft(x,K).*fft(h,K) , K ); % y[n]:


Taking the area under a square pulse Hey Stack Overflow community,

I've been given a series of datasets that resembles a square pulse, and ideally looks like the following:

This particular plot is composed of many sets of data.

I'm wondering if there is a way to reduce the noise found in the data, and retrieve just the area under the square pulse with the assumption that "spikes" may appear before or on the actual pulse itself.

The method I've tried is to cut everything before and after the pulse, and look for spikes that may occur on the actual pulse by looking at how much each value deviates from the mode peak found. However, this seems less than adequate as my implementation requires going through each point multiple times generating quite a large time complexity, and the file sizes are pretty huge.

Thanks for helping out!

A good strategy in the presence of white noise would be to do a moving average over the number of samples that equals the length of the pulse. The peak of this moving average will be your area with the noise optimally minimized.

So if you use MATLAB or Octave it would be as simple as this, assuming the pulse is 200 samples wide:

out = filter(ones(200,1),1,in);
area = max(out);

Or otherwise convolve your data with a sequence of 200 ones.

Using the Python scipy.signal package this is equivalently:

out = sig.lfilter(np.ones(200,1),1,in)
area = np.max(out)


Correlation between 2 signals of uneven dimensions As a part of my work, I am trying to correlate the audio signal in a video with the pixels of each frame. The steps I follow are:

Audio sampling rate and frame rate of the video are known. So, division of the audio sampling rate with the frame rate gives us the number of audio samples associated with each frame of the video.
That being done, I want to correlate the audio signal (the collection of audio signal amplitudes with 2 channels) with the 2-D pixel intensity of each frame, to basically decipher whether which part of the frame has maximum correspondence with the audio.
But, obviously the 2 signals ( audio signal with 2 channels and the 2-D frame pixel intensity values ) are having different dimensions. In this scenario, how do I compute the correlation between them?

I did encounter the idea of performing canonical correlation, but then again, the unequal dimensions of the 2 signals makes it a problem.

The intuitive way to go about this would be to consider every single pixel in your video (assuming that's only intensity, not e.g. color) a 1D signal over time.

Then you'd get width×height number of crosscorrelation functions between pixel intensities and your audio.

The different sampling rates of audio and video just mean that you'd need to reduce your audio to a signal of a rate that is identical to the video frame rate. (You could also interpolate the video to the audio rate, but then you'd be mostly making statements about the correlation properties of your interpolation function, not of your video).

The question is how you can reduce your audio to something that only takes say 25 values in a second. My wild guess is that this is something where you either use your inherent ideas of what would be a clever measure (e.g. "loudness in 40ms steps") or you just do the same for all 12.5 Hz wide bandwidths of your audio, and see which correlates best.


What is meant by "arbitrary" in the context of digital filter design? I try to understand the 'arbitrary' what does it meant?
I had read many references ,such a one it is 'begin the process with a transfer function of your choice' ,in other reference :related by the starting point,

In other hand i read about Chebyshev approximation with arbitrary magnitude,
how can i find the better design with arbitrary magnitude ?
any method of design have a criteria ,the 'arbitrary' related by the method of design or related by the transfer function?

The term "arbitrary" in this context refers to the desired frequency response, and what it means is that the respective design method accepts any desired response. So unlike required by many standard design routines, the magnitude response does not need to be piecewise constant or linear, and, more importantly,  the phase response need not be linear or minimum-phase. The corresponding design problem is a complex approximation problem because the desired frequency response which must be approximated is a complex-valued function.


How to detect underdamped harmonic signal? I would like to ask what could be the easiest way to detect the presence and beginning (no more, no less) of under-damped harmonic signal in a collect of ultrasound data. I have some acquired data which is presented in 2-D image form, where vertical axis are related with time. Under-damped harmonic signal, which I need to detect, looks like here, but of course it is noised.

I would like to detect presence of this kind of signal automatically. Currently I've developed algorithm, which calculates short-time Fourier transform for every vertical line of image, collect values of five the most often frequencies for every short-time window, and chooses the vertical line with the most harmonically nature - with the most repeated frequencies. Unfortunately in this solution every time one of the line is chosen. I can use some kind of threshold for sure, but i don't think so it's good solution. However, it isn't the most important problem - my algorithm calculate which line has under-damped harmonic signal, no where it starts. Second problem is that this solution is too slow for my purpose - 20 times per second.

I have tried to make my own research but solutions like Kalman filter, or estimation by least squares method, in my opinion, it's like attempting to kill ant by tank - I need smaller caliber method. Currently I don't need to know features of signal like values of angular frequencies and amplitudes.

Any ideas signal-processing community? Thanks in advance.

First, de-noise the waveform using a technique suitable for non-stationary signals. Admittedly, the spectrum of the signal doesn't change too much over time. There are efficient iterative de-noising algorithms for this purpose, such as basis pursuit de-noising [1]. But even a simple low-pass filter may be sufficient depending on how bad the non-stationarity is.
Since your record is mostly oscillatory, the underdamped signal effectively begins when the successive peaks of your signal are lower in amplitude. With your de-noised smoothed signal, attempt to pick up local extrema. This can be done in real-time by looking at the first derivative or discrete difference. Once you have extracted the extrema, using the values of the extrema, you can construct the signal envelope. Empirically, you can decide that if the envelope is decaying, then it signifies the starting point for under-damped behavior.

Fourier and short-time Fourier transforms may not be as effective in this problem since you are looking for the exponential decay feature, which is a temporal signature.

[1] Chen, Scott Shaobing, David L. Donoho, and Michael A. Saunders. "Atomic decomposition by basis pursuit." SIAM journal on scientific computing 20.1 (1998): 33-61.


DFT independent variable as fraction of sampling rate When the frequency domain's independent variable is labeled as a fraction of the sampling rate, the values along the horizontal axis always run between $0$ and $0.5$, since discrete data can only contain frequencies between DC ($0$) and half the sampling rate.

But why is the highest frequency pegged at exactly $0.5$ of the sampling rate? Isn't it possible to have a scenario where the highest frequency is less than half the sampling rate?

E.g., consider a speech signal that has been filtered to remove all frequencies above $3.3$ kHz and has been sampled at $10$ kHz. The highest frequency ($3.3$ kHz) is less than $0.5$ of the sampling rate ($5$ kHz). How would the frequency domain of this signal (labeled as a fraction of the sampling rate) look like?

The DFT calculates spectral components up to $f_\mathrm{s}/2$, no matter what the input signal is. If the input signal to the DFT contains (absolutely) no spectral components between $0.33 f_\mathrm{s}$ and $0.5 f_\mathrm{s}$ the corresponding DFT outputs are zero (assuming that there is no spectral leakage).


Function for response over image scales I want to write code in matlab to find a characteristic scale of an image patch.
I have an image patch X, I need to maximize function F(X,sigma) so that sigma is the characteristic sigma.

For my experiments, F can be the laplacian of the image, sigma is the "power" of a gaussian filter on the image. F(x,sigma) should return for each x (image patch) its response over scales of increasing sigmas (i.e. a discrete value for each X and sigma).
My questions are:

The laplacian is obtained by convolving laplacian kernel with the image, where
does the sigma takes place ? (is it by smoothing the image before (the laplace) with G(sigma)?)
How can I define F(X,sigma) to return discrete value (X is nXm), ? i.e. what does this function do exactly.

Dima and Jean, Thank you for answers, the reason I'm asking this question is because I need to do something similar on 3D data (probably/maybe with some other normalized derivatives), and I first want to see how it works with 2D data (i.e. image).
I understand now what I need to do for constructing the scale space, there is still one thing I'm not sure of and it is finding the local maxima over F. For each image (x,y) and σ - F is actually the derivative of the image. What does is mean local maxima of F over x,y and σ ? look for local maxima(s) at each scale ? there are several local maxima at each scale, are these local maxima correspond (in position) over all scales? if so, I understand that for each such position I need to find the σ that gives best maxima over all scales..
Is this is what you meant guys ?
I'm sorry if I'm going to much into details, I just want to make sure I completely understand,
THANK YOU !

Conceptually this works as follows: first you construct a Gaussian scale space by convolving the given image with the Gaussian filters of increasing $\sigma$. This is where the sigma comes from.  Thus you have a 3D scale volume, where each plane is an increasingly blurred version of the original image.  Then you evaluate your function $F$, which in your case would be a Laplacian, over each scale plane, i. e. at each $\sigma$.  Then you find local maxima of $F$ over $x, y$, and $\sigma$.

Note that in practice you can use the fact that convolution is a linear operation in various ways.  For example, you can generate each scale plane from the previous plane, rather than from the original image.  This allows you to use the same filter, rather using several filters with increasing $\sigma$.  Or you can generate $F(\sigma)$ directly from the original image by using the Laplacian-of-Gaussian filter of the corresponding $\sigma$


Convolution : Anechoic Speech signal and Impulse response I am trying to perform convolution in MATLAB of an anechoic speech signal (2.1 s) and a impulse response (20 s) using the following code:
[y_anech,fsSig]=audioread('MA03_03.wav'); %speech signal fsSig = 48000

[y_IR,fsIR]=audioread('IR.wav'); %reads the Impulse response, fsR = 44100

%%% Resampling the speech signal
[P,Q] = rat(fsIR/fsSig);
y_anech_resampled = resample(y_anech, P,Q);

sig_conv_meas = cconv(y_anech_resampled,y_IR)/fsIR;
sig_conv_meas = sig_conv_meas/max(abs(sig_conv_meas));

player1 = audioplayer(sig_conv_meas,fsIR);
play(player1);

Running the above code, it gives me the required convolved signal, however at 20s, I also hear the anechoic speech signal. I have the following questions:

What the reason behind that I hear the original signal after the convolved signal? How should I remove it?

See the attached plots. In the last plot, around 20 s, there appears the anechoic signal:

ETC plot of IR:

When performing convolution, why do I need to divide the convolved signal with sampling frequency as I do in my code (sig_conv_meas = cconv(y_anech_resampled,y_IR)/fsIR)?

Thank you!

I have no access to your audio files so I've downloaded:

IR from here (mono/r1_omni.wav) - it's a really long one
Anechoic recording from here (operatic-voice/mono/singing.wav)

Resampled voice signals:

Final convolved signal:

As for your questions:
1.
As you did the plot of IR in logarithmic scale it's clearly visible that towards its end there is another peak at 20s. Once convolving with speech it will yield another repetition of the speech signal. If you look at the first plot, your audio starts around 1s mark. As you are convolving the IR with the input signal, it will repeat it after 20s, and this is what you see in your convolved output around 21s.
Think about how an IR would look like if you want to produce an echo after 20s - one strong peak at 0s and another, smaller at 20s. I am guessing that this "peak" at the end is, in fact, the build-up before the main peak and should actually precede the main one (you can tell it by digital zeros and low energy samples in-between), or potentially a result of harmonic distortion that was moved to the end. Definitely, someone must've extracted/saved this IR in such a way.
Nonetheless, you can easily ignore anything past the 10s mark as it doesn't carry any extra information during convolution. TBH, I would ignore the part after 5s, where it roughly meets the noise floor. Otherwise, it's a waste of computational resources to convolve with full length.
Possibly you can automate the cropping. It looks like the IR file is stored in 16-bit audio, and when the noise-floor ends (around the 10s mark), you get sample values equal to 1. If you look when they start in the vector, you can find your cut-off point automatically, should you have more IR's.
2.
First of all your division by fsIR has no effect as it's normalized by the maximum value in the line that follows. Anyway, most likely your IR is reverberant and generally increases the energy. Think of it as a filter that has a lot of frequency boost. Here is a comparison for resampled anechoic recording and unscaled auralised signal.


Bode plot from from step response data For quality control, I want to measure the frequency response of hardware I built. What I can easily generate on the hardware is a step response. Thus, I want to derive the magnitude and phase from the step response. This works well in general, however I am looking for ways to improve and understand basic relationships.
Remark: Quality control is probably done best by comparing the measured step response directly with a "golden" step response. However, I want to provide also numbers and visuals of magnitude and phase as the original system specification is given in these terms.
Find below the python code of what I have, comparing my approach with a direct Bode calculation from the transfer function.
Questions:

Is there a formula or rule of thumb of how the duration/length in time of the recorded step response and the sample rate used affect the quality/uncertainty of magnitude and phase? One obvious thing is, that Bode data for low frequencies with longer period than the recorded duration will have lower quality.
I always get a phase jump at the nyquist frequency, regardless of the sample rate and signal duration used. I am struggling to identify the root cause of this. Does it have an explanation from DSP theory or is it my implementation?
What other things should I pay attention to to get high quality Bode data?

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from scipy import signal # to generate example data

def stepResp2Bode(stepResponse, fS, minFreqRes=None):
    if minFreqRes is not None:
        if fS/stepResponse.size > minFreqRes:
            print('Minimum frequency resolution not achieved!')
            return # not enough samples for required frequency resolution
    impResp = np.diff(stepResponse, axis=0, prepend=0) # get impulse response by differentiating step response
    fftResult = np.fft.rfft(impResp, axis=0) # do fft
    bodeF = np.fft.rfftfreq(impResp.size, d=1/fS) # create frequency vector for fft-Results
    bodeAmplDB = 20 * np.log10(np.abs(fftResult)) # extract amplitude in dB
    bodePhase = np.angle(fftResult) # extract phase in radians
    results = pd.DataFrame({'frequency':bodeF,
                        'amplitude in dB':bodeAmplDB,
                        'phase in radians':bodePhase,
                        })
    return results

if __name__ == "__main__":
    # generate example data - in the application this would be real world measured data
    fS = 50 # sampling frequency
    wc = (2*np.pi*1) # corner frequency of example low pass
    t = np.arange(0, 5, 1/fS) # vector of time values
    sys = signal.lti([1.0], [1/wc, 1.0]) # generate lowpass filter
    stepResponse = signal.step(sys, T=t)[1] # calculate step response of low pass filter at given time values
    # Calculation
    bodeStepResults = stepResp2Bode(stepResponse, fS, minFreqRes=1) # from step response to magnitude and phase
    w, mag, phase = signal.bode(sys) # bode plot directly from system description
    bodeScipyResults = pd.DataFrame({'frequency':(w/(2*np.pi)),
                    'amplitude in dB':mag,
                    'phase in radians':(phase/180*np.pi),
                    })
    # plot results
    fig = plt.figure(1)
    fig.clear()
    ax1 = fig.add_subplot(211)
    ax1.semilogx(bodeStepResults['frequency'], bodeStepResults['amplitude in dB'], '-', label='from step response')
    ax1.semilogx(bodeScipyResults['frequency'], bodeScipyResults['amplitude in dB'], '-', label='from signal.bode')
    ax1.set_xlabel('Frequency  in Hz')
    ax1.set_ylabel('Amplitude  in dB')
    ax1.grid(which='both', axis='both')
    ax1.legend(loc='upper center')
    ax2 = fig.add_subplot(212)
    ax2.semilogx(bodeStepResults['frequency'], bodeStepResults['phase in radians'], '-', label='from step response')
    ax2.semilogx(bodeScipyResults['frequency'], bodeScipyResults['phase in radians'], '-', label='from signal.bode')
    ax2.set_xlabel('Frequency  in Hz')
    ax2.set_ylabel('Phase in radians')
    ax2.grid(which='both', axis='both')
    ax2.legend(loc='upper center')
```

First of all - I'm not familiar with Python and didn't check your code.
1a) For FFTs it is the usual that the frequency resolution equals the inverse length of your data record.
1b) As far as I understand your code (I may be wrong here), you take the FFT of the step response. The Bode diagram usually shows the FFT of the impulse response. Take care - the result may not be what you expect.
1c) In practice (due to transients, noise,...) the usual approach to measure the FRF is

excite the system with a broadband signal, estimate the spectral density of the input and the output signal, divide the output spectral density by the input spectral density. This approach is well covered in the literature (and called "spectral method" or something similar

exite the system with a sinuosoid, wait until the transient has died out, and read the magnitude and phase from your data. Repeat for every frequency where you want to know the FRF.

There are other approaches to directly measure the FRF, like the Local Polynomial Method (or Local Rational Method), TRIMM, and more - they are covered in (scientific) publications, together with a lot of beautiful maths. The approaches I mentioned are called non-parametric methods (there are also parametric methods, but that may lead to far away for your use case)

The phase jump at Nyquist frequency is to be expected. The FFT spectrum of a rational signal will be symmetric to the DC-frequency - the magnite is axisymmetric, the phase point symmetric. For some FFTs the "negative" and "positive" frequencies are wrapped (this may depend on your FFT implementation). The result is a "jump" the phase.

The input signal should excite the entire frequency range of interest - check the spectral density function of the input signal (i.e. the step).


Which filters are the most commonly used and easier to implement for ECG noise removal? ECG signal is affected by interferences such as the Powerline Interference, Baseline Wandering and Muscle Noise.

There are many filters available which can reasonably cancel these types of noise.

I've found that Butterworth, Chebyshev, Wiener and Kaiser filters, along with wavelets could be used for this purpose.

Now, my main questions are:

Which filters are the most commonly used and easier to implement
(both from the above list and from other filters you have in mind)?
Can you provide me some pseudo-code samples on how to implement some simple filters, maybe Butterworth or Wiener?
In which references should I look at concerning ECG digital
filtering, I mean which provide an overview of the various filters instead of foccusing on a specific one?

Thanks in advance

Different filters are used for different purposes. Removing powerline noise, which usually is concentrated at a single frequency, is usually performed with a notch filter. Removing baseline wandering is usually done with a DC blocking filter, which is a very narrow high pass filter, and is essentially another notch filter with 0 Hz center frequency. Removing EMG is usually more difficult, but you can try doing it with a low pass filter, since EMG has a lot of content in high frequencies.

In order to process EEG signals for simple things like heartbeat rate detection, look at the Pan-Tompkinks algorithm.

Questions 2 and 3 require some effort on your part, since they're very broad. I've used this book in a graduate class on the topic, and I found it to be very insightful.


Create struct from column vectors I do not know how to create sturct variable in Maltab from column vectors. I have three column vectors, Voltage1 (32768x1 double), Voltage2 (32768x1 double) and Voltage3 (32768x1 double). I would like to create 1x1 struct variable with 3 fields (Voltage1-32768x1 double, Voltage2-32768x1 double and Voltage3-32768x1 double). Please help.


How do I show if this system is time-varying? $\displaystyle y(t) = x(t), t \ge 0$ and $-x(t), t \lt 0$

I can anticipate this system is time-varying.

However, I cannot be sure it and prove it.

please help me.

In order to show that a system is time-varying, it is sufficient to find one counter-example, where the response to a shifted input signal does not equal a shifted version of the original output signal.

In your case you might want to consider the signal

$$x(t)=\begin{cases}&1,& 0<t<1\\
&-1,& -1<t<0\\
&0,&\text{otherwise}\end{cases}$$

Now try to figure out the corresponding output signal $y(t)$. Next, take a shifted version of $x(t)$, e.g. $x(t-1)$, as the input signal and see what output signal you get. Is it equal to $y(t-1)$? If not, you've just shown that the system is indeed time-varying.


Why FFT of a rectangular pulse has imaginary part while the sinc function is real? I have a rectangular pulse:

w = 20; pln = 10;
v =[zeros(1, pln), ones(1, w),zeros(1, pln)];
figure(); plot(v,'-o');

by perform fft and fftshift the pulse center to zero to get its frequency domain representation:

yf = fftshift(fft(v));
xf = [-0.5000   -0.4750   -0.4500   -0.4250   -0.4000   -0.3750   -0.3500   -0.3250...
-0.3000   -0.2750   -0.2500   -0.2250   -0.2000   -0.1750   -0.1500   -0.1250...
-0.1000   -0.0750   -0.0500   -0.0250         0    0.0250    0.0500    0.0750...
0.1000    0.1250    0.1500    0.1750    0.2000    0.2250    0.2500    0.2750...
0.3000    0.3250    0.3500    0.3750    0.4000    0.4250    0.4500    0.4750];
yf = yf.*exp(-2j*pi*xf*(pln+w/2)); % this shifts the time domain signal to zero

yf is close to the sinc function - the Fourier transform of a rectangular function - but not exactly the same, it has imaginary part while sinc is is purely real:

yf_sinc = w*sinc(w*xf);

Why is the difference and why ifft of samples taken from a sinc function will not produce an exact rectangular pulse but with ripples close to the edge while ifft(fft(v)) will exactly reproduce the original rectangular pulse:

% has ripples around the edges
figure(); plot(fftshift(ifft(fftshift(yf_sinc))));
% reproduces the exact rectangular pulse
figure(); plot(fftshift(ifft(fftshift(yf))));

herer yf, if we assume it was directly sampled form a underlying function, that is ingoring the fact that it is produced by fft of a rectangular pulse (it's just some numbers anyway), what function would it be? Is it possible by sampling a function in the frequency domain (not need to be the sinc function) and take the ifft will produce an exact rectangular pulse?

The even part of a real signal transforms to the real part of the DFT.  The odd part of a real signal transforms to the imaginary part of the DFT.  For example, the DFT of cosine, an even function, is strictly real, while the DFT of a sine, an odd function, is strictly imaginary.


DAC and sample and hold effect I am having some trouble understanding how a digital signal gets converted back to an analog signal.

Why is the "sample and hold" effect a problem?

so the sampling theorem says basically:

let   $ \ x[n] \triangleq x(nT) \ $  where $ \ f_s = 1/T \ $   is the sampling frequency.

then $ x(t) $ can be reconstructed from the samples $x[n]$  using

$$x(t) \ = \ \sum_{n=-\infty}^{\infty} x[n] \cdot \mathrm{sinc} \left(\frac{t-nT}{T} \right) \ $$

which is what you get when you pass this ideally sampled signal

$$  T \sum_{n=-\infty}^{\infty} x[n] \cdot \delta(t-nT) \ $$

through an ideal brickwall low-pass filter with cutoff at $f_s/2$ and a passband gain of 1.

but what comes out of a DAC is not that, but is this:

$$  \sum_{n=-\infty}^{\infty} x[n] \cdot \mathrm{rect} \left(\frac{t-T/2 -nT}{T} \right) \ $$

so it's like the DAC has a hypothetical filter with impulse response of

$$ h(t) \ = \ \frac{1}{T} \mathrm{rect} \left(\frac{t}{T}-\frac{1}{2} \right)
 = \begin{cases}
\frac{1}{T} & \mbox{if } 0 \le t < T  \\
0           & \mbox{otherwise}
\end{cases} \ $$

so the "sample-and-hold effect" really should be called the "zero-order hold effect" and what the effect is, is to pass your signal through a filter with the above impulse response and this transfer function:

$$ H(s) \ = \ \frac{1 - e^{-sT}}{sT} \ $$

and this frequency response:

$$ H(j 2 \pi f) \ = \ e^{-j \pi f T} \mathrm{sinc}(f T) \ $$

now, these equations came out of the wikipedia references, which is why i suggested to look there first.


How to simulate a signal of a particular frequency according to the number of points in a period of the signal? I have one period data ($n$ points in total) of a signal. I want to generate a periodic signal with a total of $m$ points and corresponding signal frequency of $F$ according to this data. How do I do that? Besides, I want to design a low-pass filter according to the simulated periodic signal, $F_s$ and $F_{stop}$ need to be set in fdatool. How to calculate $F_s$? Thanks!

I suppose that the $n$ points (hereafter, the chunk) precisely divide the base period (which you don't know yet)  evenly. You are lucky. You can verify that by concatenating  the $n$  discrete samples  (or possibly $n-1$ consecutive ones)
by an integer number $K$ of (periodic) repetitions (copy-paste, several times, hereafter the chain), with $Kn\ge m$ (or $K(n-1)\ge m$. It should be smooth  visually. The FFT of the chunk and the chain should be superposable (with appropriate scaling in time and amplitude).
Then, crop the chain of $Kn$ samples to $m$.
Then, using a non-parametric (discrete Fourier transform, FFT), you can estimate the spectrum of the (underlying) continuous data, and find the minimum Nyquist sampling frequency, which may differ from the original one (the initial sampling).


Making objects stand out? Lets say I have two objects very close to each other. There is a colour difference between them, but the difference is small. Therefore when finding contours, both the objects are treated like on object. Is there a way to make the difference in colour more significant? As in exponentially increase the small difference in colour? Or if not possible how would I make these objects stand out from each other?

I am using OpenCV please suggest solutions possible in that library.

Thank You

Consider converting the image into L*a*b* color space — colors with subtly-different component values in RGB color space will often have significantly-different a*b* component values in L*a*b* color space.


FIR filter compensator when using a CIC decimation filter When decimating a narrow-band signal with a cascaded integrator-comb (CIC) filter, which  FIR  filter is more suitable to compensate the CIC frecuency response?

There was a similar question, https://dsp.stackexchange.com/a/1551/306, and the following is a subset of the answer from the other post.

Generally, to compensate a CIC filter the inverse of the CIC filters response
can be used to generate the compensation filter.  The CIC has
a response of 2

$$
  H(\omega) = \left| \frac{sin(\omega D/2)}{sin(\omega M/2)} \right|^N
$$

Where D is the differentiate delay, M is the decimation rate, and N is the
filter order (number of cascaded filters).  The inverse can be specified as

$$
  H(\omega) = \left| \frac{sin(\omega M/2)}{sin(\omega D/2)} \right|^N
$$

Once we have the
frequency response of the compensation filter, we can simply choose the length
of FIR filter that we desire.  The length of the FIR is application
specific.  Obviously the longer the FIR filter the better compensation.

The following are plots of this straight forward compensation.

The following is the Python code to create the frequency responses and plots.

    import numpy as np
    from numpy import sin, abs, pi
    import pylab

D = 1; M = 7; N = 3

Hfunc = lambda w : abs( (sin((w*M)/2)) / (sin((w*D)/2.)) )**N
HfuncC = lambda w : abs( (sin((w*D)/2.)) / (sin((w*M)/2.)) )**N

w = np.arange(1024) * pi/1024

G = (M*D)**N
H = np.array(map(Hfunc, w))
Hc = np.array(map(HfuncC, w))
# only use the inverse (compensation) roughly to the first null.
Hc[int(1024*pi/M/2):] = 1e-8
plot(w, 20*log10(H/G))
plot(w, 20*log10(Hc*G))
grid('on')

See, Altera, "Understanding CIC compensation filters", for other approaches and $sinc^{-1}$ approximation 1.

1 Altera, "Understanding CIC compensation filters"

2 R. Lyons, "Understanding Digital Signal Processing",
    2nd ed., Prentice Hall, Upper Saddle River, New Jersey, 2004


Data measuring interface I have an hardware related question.

I want to create a measurement set-up so that I can record data for DSP analysis. When searching around I see I have a few different options. The ones I have found are:

Data acquisition box / interface
Logic analyzer
Scope
Soundcard (have specific input type (XLR/ TRS-Jack) )

From DSP perspective what are the differences between these types of measuring devices? Specially the difference between a DAQ and a logic analyser?

As Questioned in the comments: what are the parameters i want to measure. Here is the list for now(witch can change later on):

Torque
Acceleration(vibration)
(velocity)
Force
if possible Ammeter
different switches for safety sensors(just low and high over time)
some pulse generators

I'm not sure what the bandwidth is of every signal. But i assume that the maximum bandwidth is from the acceleration sensors. I assume that the bandwidth is f_max = 10kHz so f_sampling = 25 or 30 kHz.

Eventhough it's not fully clear what sort of an environment you are referring to by the measurement set-up, in the most broadest sense of the term you need either a data acquisition interface or a sound card for the cheapest semi-alternative.

A data acquisition interface to a PC is what you mainly need. There are a bunch of alternatives. You must set your requirements (data rate, #channels, bit depths, ADC linearity, SNR etc.) to get the optimal match for the price. Note that the software interface to the acquisition hardware is another concern that makes big functional differences in the end.

The cheapest alternative to a professional data acquisition hardware can be a PC audio soundcard. They will record AC coupled and distorted signals. So no DC is captured unles you modify (if possible) by shorting the soundcard's ADC chip's input-DC blocking caps. There are a number of soundcard capturing software which also simulate an oscilloscope as well. Most typical PC (laptop) sound cards do provide 16 bits, 2 ch, 96 kHz sampling these days.

A logic analyser is almost useless from a DSP point fo view. Logic analysers are used for digital electronics such as logic gates, microprocessors, interface protocols etc. to analyse their bits, high and low values, ones and zeros.

An oscilloscope is again useful mainly for analog electronics diagnosis purposes. It's an indispensible hardware tool for any hardware project, but not exclusively required for data capturing, unless you want to observe in real-time what you are capturing.


How virtual ads work in stadiums? I was looking at a technology which enables virtual advertisement on the perimeter ads in soccer. It's relatively an old technology for a few years back, but I'm surprised to see the quality of it. They overlay the physical perimeter ads with new ads on different camera feeds. Basically from this link. And here is a video on Youtube.
LED perimeter system with integrated virtual replacement technology - makes it possible to change the content seen by different viewing audiences around the world
How does this overlaying mechanism work? Is there any hardware involved (like the LED boards itself), or it is purely based on Computer Vision and Augmented Reality? How can it be such high-quality without glitches purely based on computer vision?!

Since "green grass" and "viewers/whatevers behind the LED boards" do nicely contrast from the boards themselves, and because the shape of these is pretty precise known beforehand, I'd probably go with "pure computer vision". You could even make the image processing aware of what's displayed on the LEDs, so that recognition gets even easier.
Of course you could do something like blinking the LEDs at a specific rate and synchronize cameras to that blinking, so that it gets even easier on the processing side, but honestly, I don't think that's necessary. Might still be done to make the whole thing failsafer.


No-lag causal bandpass filter kernel formula Can anyone point me towards a formula for a good (for any values of "good") no-lag causal bandpass filter?

I am doing sound processing and I need to identify starting points of certain sounds (target precision <10 ms.) I've been using a Lanczos filter followed by convolution with a sine wave to do spectral analysis, but that gives me skewed results because amplitudes at any given time are being affected by stuff that's happening at later times (the filter is non-causal).

Simply cutting off the right half of the Lanczos filter resolves that particular problem, but the resulting spectrum is badly messed up (I get heavy sideband spillover). All other simple filters I tried are even worse.

Half an hour of googling turned up several mentions of Butterworth filters, lots of transfer functions (which I could not care less about), but no formulas of actual kernels.

EDIT: Let me try to be more specific. I start with a narrow-band filter with impulse response:

$$ h(t) = \begin{cases}
 C \, e^{i 2 \pi t/T} \operatorname{sinc}^2\left(\tfrac{t}{NT}\right) \qquad & |t| \le NT  \\
 0 &  |t| > NT \\
\end{cases} $$

where $\operatorname{sinc}(u) \triangleq \frac{\sin(\pi u)}{\pi u}.$

With $T=$2 ms and $N=$10, that gives me 0 dB at 500 Hz and <-40 dB below 440 or above 560 Hz. The support is $2NT$=40 ms. If I use it on an input that goes from zero to a pure 500 Hz sine wave at time $t_0$, response goes from exactly zero to just over zero at $t_0$, to half the maximum level at $t_0+NT$, and stabilizes when the entire area under the filter is a sine wave. The "effective" lag (time to 50% level) is 20 ms.

What I'd like to do is to bring this "effective" lag as close to zero as possible, while maintaining the overall discrimination ability (being able to distinguish between a 440 Hz sine wave and a 500 Hz sine wave.)

What you need is not a "no-lag" filter, as such a thing doesn't exist if the filter is still causal (how should it filter, without knowledge of the future and not being aloud to wait).

What you need is a filter with a constant group delay. Any linear phase filter has that. Any time-symmetric filter has linear phase. And by the simple method of shifting a time-symmetric filter by half its length, you make it causal. And have half its length in delay, at any frequency.

So, any filter kernel that is symmetric works. There's far too many bandpass filters that have linear phase, and which one you'd use completely depends on your application. But it's usually easier to find a tool to design a linear-phase filter than to find one that gives a different kind of filter, so you should be fine.


Reliability of normalized cross-correlation with high lags Say you have two signals with 100 data points in each signal. Then there are 20 time lags where the computation of the normalized cross-correlation between those signals only considers 10 or less pairs of data points. It seems like the normalized cross-correlation is less meaningful at these lags because it is computed from less data. My question is, is there a general way of accounting for the uncertainty introduced when we look at the correlation with higher lags? Or does that vary too much based on the context?

Describing "correlation" as just a sample by sample multiplication and then a sum of the products of your two data sets [note that I am calling them data sets, not signals to avoid later confusion], correlated components will increase in magnitude at rate N (20Log N in dB), while independent components will increase in magnitude at the square root of N (10Log N  in dB). If we consider the correlated components as the "Signal" of interest and the independent components as noise, we can see that we get a 20Log N - 10Log N = 10LogN increase in SNR due to the process of correlating.

With that in mind, consider your example if the first data set is a signal of interest "buried" in noise. Let's say it is a pseudo-random sequence, which is nice for this explanation in that we will get a sharp correlation at one delay and minimum correlation elsewhere. The second data set is the pseudo-random sequence itself. We correlate to find the presence of the sequence in the noisy data set, along with the delay (the lag where we get maximum correlation). Now here you should see the answer to your question in light of my description above:
Assume the sequences are completely aligned in the two 100 sample datasets (lag=0), in this case we get 100 samples where the signal accumulates in the correlation, to 100x= 40 dB, while the noise component accumulates to sqrt(100x) = 20dB, and therefore you got a 20 dB net processing gain in SNR due to your 100 samples (important! Assuming here the noise from sample to sample in the first data set is independent, meaning white in frequency).

Now consider if the sequences were not aligned until you are at a lag of 90 samples, such that only 10 samples really are multiplied and then summed. In this case the signal accumulates by 10x or 20dB and the noise accumulates by sqrt(10) = 10dB for a net processing gain in SNR of 10 dB.


How to verify if my cross-correlation algorithm works? So I have a signal generator which can produce two separate signals.

I use a spectrum analyzer to read these signals, and can interface with Python.

I wrote a Python code which takes in these two signals, perform the FFT of both, take the complex conjugate of one of them, and multiply together which gives me the cross-correlated frequency spectrum?

My problem is, I want some validity checks to see if this works or not. The main objective is noise floor reduction, but I am not sure how much the noise should drop by, etc. This signal generator is connected to the spectrum analyzer using BNC coaxial cables and that is it. I am sending in simple sinusoidal waves of frequency 2kHz.

There are two issues / concerns with this approach in that you may not be getting what you want. The primary one is the equivalent of post detection averaging by using the post-processed results from a spectrum analyzer, and the second is that the result of your complex conjugate multiplied FFT is actually in the time domain since you started in the frequency domain (so the FFT processing is actually converting to the the time domain where you perform the multiplication.-- I suspect but am not certain that your goal is to observe the frequency domain spectrum after cross-correlation of your signals?).

Most significantly, spectrum analyzers "detect" the signal in that the result is dB magnitude, and regardless of the dB conversion which could easily be undone, all phase information is lost such that any processing gain through coherent averaging can no longer be done. This is equivalent to the difference in adjusting the resolution bandwidth (RBW) versus the video bandwidth (VBW) of the spectrum analyzer . Since you are already using a spectrum analyzer, playing with those knobs will actually provide a great demonstration of pre-detection versus post-detection averaging, which will give you much further insight into what you are trying to accomplish. For that reason, I will first go into that in more detail of that, starting with a simple functionally equivalent block diagram of a spectrum analyzer:

In this simplified but functionally equivalent view of a spectrum analyzer an input signal (in this diagram a 70 MHz tone) is frequency translated by the voltage controlled oscillator (VCO) and mixer to be swept through the Bandpass Filter with resolution bandwidth RBW. It may be helpful to envision the sweep rate as being extremely slow in comparison to the inverse of the filter bandwidth, such that the operation for explanation here is similar to stepping the VCO through each frequency one at a time and statically processing each result (pixel on the display) before stepping to the next frequency. This will avoid us getting into constraints on the sweep rate which exist. Thus at any given moment, the power detector is converting the entire power that is within the bandwidth RBW to a power level. This power level passes through a low pass filter with bandwidth VBW, in converted to dB and then used to control the vertical position of the pixel on the screen representing the power level of our signal under test. Similarly the ramp rate controls the horizontal position representing the frequency at that moment in time.

By adjusting the bandwidth RBW we reduce the total power that would be presented to the power detector. If the signal is centered on our tone (in this case 70 MHz), the power of the tone would dominate and we would see no change in the power detector output regardless of the bandwidth of RBW, as long as the bandpass filter remained centered on our signal under test, and assuming the signal under test was a pure tone. In this case the apparent bandwidth of that signal, as displayed in the diagram above is actually the bandwidth of RBW shown as the single tone is swept through our filter during the horizontal trace. This would be very apparent by adjusting RBW and observing the result. More notably you would see the noise floor everywhere else go down at $10Log(N)$ where N is the decrease in RBW. (If you half RBW, the noise floor goes down by 3 dB). I believe this is a direct demonstration of the noise floor reduction you seek to that the extent correlation has mathematical similarity to averaging, as in this case you would observe that the signal level would remain unchanged. In this case, this is the result of "Pre-detection" averaging which is exactly what a bandpass filter is doing (with the averaging operation mathematically translated to the center frequency of the filter). In the plot shown in the figure above, the noise floor is approximately -90 dBm. Both the noise floor and the width of the tone suggest that the RBW is about 1 MHz, so I will assume that was the case. If we adjusted RBW to 500 KHz, we would see the noise floor drop to -93 dBm while the peak of the tone would remain unchanged.

Similarly if you adjust the bandwidth of the lowpass filter (VBW), you are simply smoothing (averaging) the noise power already measured but not decreasing it! In the plot above we would see the "noise on the noise" get smoother, but it will remain at the -90 dBm level. Thus we are simply averaging the result of the noise without reducing it.

To accomplish what you want, I would consider cross-correlating the signals directly in the time domain. The complex conjugate multiplication of the FFT of these time domain signals would then represent the spectrum of the cross correlation, if that is what was ultimately desired.


Cross-correlation peak How get cross-correlation peak and based on it calculate correlation score for similarity of two audio samples. SO far I've

FFT two samples
complex conjugate second
multiply results
IFFT
cross-correlate with itself(autocorrelate)

Thanks for any advice

As Matt stated, you should use the correlation coefficient!.

Points 1 to 4 calculate the crosscorrelation. From that you have to find the highest peak (or lowest, if it has a higher absolute value). This value is the value of the nominator.

The denominator consists of the two autocorrelation values. Those are obtained by using the same algorithm where both signals are the equal. Here the peak should be in the middle (t=0) as already stated by welcomedungeon. Taking the square root of both autocorrelation values and multiplying them, gives the denominator.

Edit:
Maybe this description is more clear:

$\frac{max(abs(ifft(fft(x_1)*fft(x_2)')))}{sqrt(max(abs(ifft(fft(x_1)*fft(x_1)'))))*sqrt(max(abs(ifft(fft(x_2)*fft(x_2)'))))}$

The apostroph means conjugate complex.

Edit: Two examples with Matlab code:

Using the same signal:

x = rand(1000,1)-0.5;
max(abs(ifft(fft(x).*fft(x)')))/(sqrt(max(abs(ifft(fft(x).*fft(x)')))).*sqrt(max(abs(ifft(fft(x).*fft(x)')))));

gives 1;

Using a sine and a cosine, should also give 1 because they are delayed versions of each other:

x = sin([0:pi/100:10*pi]);
y = cos([0:pi/100:10*pi]);
max(abs(ifft(fft(x).*fft(y)')))/(sqrt(max(abs(ifft(fft(x).*fft(x)')))).*sqrt(max(abs(ifft(fft(y).*fft(y)')))))

gives approximately 1

Using windowing before transformation to frequency domain should improve results.


How can I realize a complex communication system in a real time application? So I am working on a project where I have to model a communication system in MATLAB, then implement it on an FBGA.

My transmitter goes something like that:

M        = 4;   N = 1000;
Symb     = randi([0,M-1],N,1);
Mod_Symb = qammod(Symb,M);

Now, as you can see the output is complex consisting of four possibilities [+1+1i, +1-1i, -1+1i, -1-1i].  The question is (which I don't know after 4 years mobile com. school!) how do I go from those complex representations to a real time system? I mean, when I want to feed my antenna amplitude data, how do I do that?

The output of qammod in your program is a sequence of amplitudes. If $a_k$ is the real part of the sequence, and $b_k$ the imaginary part, then the QAM baseband signal is:
\begin{align*}
s_{BB}(t)&=\sum_k a_k p(t-kT_p) + j\sum b_k p(t-kT_p) \\
&= A(t) + jB(t).
\end{align*} In this equation, $T_p$ is the pulse (baud) rate, and $p(t)$ is the pulse shape.

In order to convert this to a real, physical passband signal with carrier frequency $f_c$, you upconvert the baseband signal to frequency $f_c$: $$s_{PB}(t)=\text{Re}(s_{BB}(t)e^{j2\pi f_ct}).$$ While the previous equation is mathematically nice, in practice this is implemented as: $$s_{PB}(t)=A(t)\cos(2\pi f_ct)+B(t)\sin(2\pi f_ct)$$ with possible scale factors to keep the signal energy constant.

Once you have $s_{PB}(t)$, all that remains is amplify it, possibly bandpass filter it to eliminate out-of-band emissions, and connect an antenna.

Note that Matlab works in the discrete-time domain and the signals I've described are continuous time. At some point in your system, you will need to do a digital to analog conversion. Exactly where that will happen depends on many factors. I recommend that you take a look at some common RF architectures, such as heterodyne, homodyne, zero-IF and so on. They can give you ideas on how to implement the analog part of your system. One good reference is "RF System Design of Transceivers for Wireless Communications" by Q. Gu. Another very good book is Steven Tretter's "Communication System Design Using DSP Algorithms: With Laboratory Experiments".


adjusting swept signal equation When I do a spectrogram in matlab / octave I can create a swept signal that looks like the RED plot line below.  But how can I create a swept signal like the BLUE line in the 1st plot using the equation below.

thanks to Daniel and David for getting me this far with the code is below

startfreq=200;
fs=44100;
endfreq=20;
dursec= 10;%duration of signal in seconds
t=(0:dursec*fs)/fs; %Time vector
alpha=log(startfreq/endfreq)/dursec;
sig = exp(-j*2*pi*startfreq/alpha*exp(-alpha*t));
sig=(sig/max(abs(sig))*.8); %normalize signal
wavwrite([sig'] ,fs,32,strcat('/tmp/del.wav')); %export file
specgram(sig,150,400);

1st plot

2nd plot

How can I fix the the equation in the variable sig to get it to look like the BLUE line in the 1st plot?

3rd plot

The general formula for the exponential rise in frequency is given by

$f=A+B(1-\exp(-\alpha t))$

We need to solve for the values of $A$, $B$ and $\alpha$. Let's assume you want an initial frequency of $20$ at time $t=0$, and an asymptotic frequency of $200$ Hz at time $t=\infty$.

Substituting $t=0$ and $f=20$ into the above equation we find $A=20$. If we substitute $t=\infty$ and $f=200$ we find that $B=180$. For $\alpha$ you really have a free choice - it controls how fast the frequency rises. If you have another control point say at time $t=2$ we want the frequency to be at $95%$ of the final frequency i.e. $f=.95 \times 200=190$ then we can substitute these into the above equation (Since we know what $A$, and $B$ are).

$190=20+180(1-\exp(-2\alpha))$

Solving gives $\alpha= -0.5\ln(10/180)=-1.445$. If you want a faster rise then you'll need to choose a larger value of $\alpha$.

To find the expression for phase we simply integrate with respect to time to give

$\phi =2\pi(At+Bt+\frac{B}{\alpha}\exp(-\alpha t))$


How to calculate Signal-To-Noise Ratio I am doing some analysis in python working with a photodiode. The signal is more or less periodic with some phase shifting in certain areas that give rise to faster or slower periods but overall still very periodic.

I am trying to see how "good" our signal is by computing the Signal-to-Noise ratio (SNR), but so far I have only read from a book on DSP that SNR is defined only as the mean divided by the standard deviation of the signal. Also the book defines Coefficient of variation (CV) as the standard deviation divided by the mean, multiplied by 100. They say a good signal has high SNR and low CV.

Is this simply the definition of SNR or is it something that should be adjusted depending on the signal? How would i go about calculating SNR for my signal?

If you can get input signal when it doesn't contain any useful signal (I mean only noise is presented), you can estimate average noise power at first. Simply find a power of such a signal:

$P_n=1/N \cdot \displaystyle\sum_{k=0}^{N-1}|s(k)|^2$.

Choose $N$ in $2^{12}\ldots 2^{15}$ range for example.

Then you can measure a power of $signal+noise$ mixture, $P_s$, by the same way, when input signal is presented at the input. Caution: the measurement of $P_s$ is only valid if the signal exists all the measurement time. SNR could be calculated by:

$SNR=10 \cdot \log_{10} \frac {P_s - P_n}{P_n}$.

If your noise is stationary and white it's possible to measure it once in the setup and use $P_n$ as a constant.

If you can't measure noise power at any reason, you can try to estimate SNR in the frequency domain via FFT. But you have to know your signal's band. The point is the same but you should use FFT frequency indexes instead of sample numbers $k$ while estimating integral power. $P_s$ is now power of $signal+noise$ in the band of interest, while $P_n$ is the noise power in total band ($0...fs$). In this case power spectral density of the noise in your system should be evaluated at first. This method could be far too complex for your task.

Both methods show you pointless result in the case of low SNR in the signal band (not total band in general). Frequency domain method will be preferable if your signal's band is quite less than $f_s$, otherwise time domain method will be good enough. To estimate SNR of the order of 5 dB or lower it's better to know signal and noise power (or power spectral density) a priory. So the task has robust analitical solution.

Hope this helps.


Real time slope and peak detection and calculation I have a signal that I sample at 500khz. I am trying to detect a rise, fall and the peak in the incoming data. The base of the peak could be  for 250 usec or 2.5msec, amplitude could be 6db or 15db above the noise floor. I don't have good snr unfortunately.  The dc level of the signal is not constant but move much more slower  than the ac component.

At the decision point, I need to know the slope of the rise and fall. This is a hard realtime system and I really need to make a decision in the 100usec after the downward slope reach to dc level.

I am looking for suggestions how can I efficiently implement an algorithm that is decent.

Currently I do a running average (past 25 data points added together) and try to detect the trend. Once I detect the trend up I start seeking trend down and once I do that I collect perhaps another 50 samples and start calculating.

Noise now easily screws this algorithm, hence the question.

UPdate

For the benefit of others, I have end up implementing a Moving Average followed by integrator. Moving average of past 64 data smoothed enough but lost rise to a degree, integrating last 8 values gained back the rise and I simply seek for rise and fall, later I did a linear regression for the slope. Works ok, not great but ok.

You should start with a bandlimited differentiator (equivalent to a differentiator followed by a low-pass filter). The differentiator will remove the low frequency trend and will respond sharply to your peaks and slopes. The low-pass component will remove noise beyond the cutoff frequency.

You should design your cutoff frequency so that you get clean pulses for your slopes.

Positive slopes will slow as positive pulses; negative slopes as negative pulses, and the peak will correspond to the zero-crossing between positive and negative.

This type of filter is normally implemented as a FIR filter. The number of samples for your filter will then depend on your real-time constraints, the sharpness at the cutoff frequency, and the cutoff frequency itself.


Projecting a 3D point into new camera coordinates Given I have:

a camera with estimated rigid motion 3x4 matrix $P = [R|t]$ that starts from the origin of the reference coordinate system
an intrinsics 3x4 projection matrix $K$
a 3D point $X$

How would I find its projection in the second camera/image? Is the following correct?

$x = K*[R|t]*X$

Your formula is correct. Here are some precisions:

the matrix $P$ takes a scene point to an image reference coordinate frame;
the matrix $K$ takes the coordinates of the point in the reference image frame and makes the necessary frame changes to fit with the actual camera (origin of the frame, pixel size...);
$X$ is a 3D point in projective coordinates, i.e., $X = (x, y, z, 1)^T$, otherwise the dimensionality of the point is not correct with the projection matrix;
$x$ will be actually a column of 3 rows: $x = (u,v,w)^T$, and the final image coordinates are given by diving $u$ and $v$ by $w$.

A few answers to your comments:

the equation assumes an ideal pinhole camera model. Thus, the coordinates are undistorted. You can actually find all these equations by yourself using basic triangular geometry and Thales's theorem;
you don't add the projective coordinate ("the 1") twice! The product $P \times (X,1)^T$ will create as output a projected (2D) point that has 3 coordinates, i.e., a 2D point in projective coordinates. What you need only is to "remove" the projective coordinate by dividing by $w$ at the end;
to get the order always right, the best way is to think of these equations as I wrote then in plain English. You take a 3D world and project it onto the world plane given by the camera position ($P$, extrinsic parameters). Then, you consider what happens for this particular camera model/sensor (intrinsic parameters). The extrinsic parameter matrix describes the camera with respect to the world, while the intrinsic parameters are always teh same wherever you put your camera.


Hilbert transform in the frequency domain As I know, the Hilbert transform
$$ H(x(t))=\frac{1}{\pi t}\star x(t) $$
in time domain is equal to
$$ -j \operatorname{sgn}(f) \cdot X(f) $$ in frequency domain. so I tried simple example using MATLAB as below,
x=[1,2,7,3];
y1=imag(hilbert(x));
f=[0,1,2,-1];
y2=ifft(-1i*sign(f).*fft(x));

but the result of y1 and y2 are different as below
y1 =
    0.5000   -3.0000   -0.5000    3.0000
y2 =
   0.5000 - 0.7500i  -3.0000 + 0.7500i  -0.5000 - 0.7500i   3.0000 + 0.7500i

just only the real part of y2 is same with y1.
Any one who knows why please explain.

The discrete-time Fourier transform (DTFT) is always periodic. This is also the case for the frequency response of the discrete-time Hilbert transformer. For this reason, the ideal frequency response is not only zero at DC but also at Nyquist, which corresponds to index $2$ for a signal of length $4$. Consequently, the correct way to do what you're trying to do is:

x = [1,2,7,3];
X = fft(x);
hil = [0,-1i,0,1i];
Y = hil .* X;
y = ifft(Y);
y =

  0.5  -3  -0.5  3


Shift a signal by fraction of a sample I have a signal with a sample time of 0.5 micro seconds and I would like to shift this signal by a fraction of the sample time, say by 3 nano seconds.

I have read a few online resources about fractional delay filtering and about using the FFT and IFFT to perform such delay. Can someone point me to some theory on this or give me an idea on how to implement it.

For the regular shifting of the signal for integer samples, I have implemented this by shifting the signal by the required number of samples and adding zeros at the beginning. Is this approach correct?

There's a good overview article which appeared in 1996 in the IEEE Signal Processing Magazine: Splitting the unit delay: tools for fractional delay filter design. The nice thing about it is that there's also a set of related Matlab files available. These routines will allow you to design such a system.

As for shifting by an integer number of samples, what you did is obviously correct.


How are time resolution and signal bandwidth  related? I am confused by the dual concepts of time-resolution and bandwidth. Often I will hear that a pulse-compressed radar application 'doesnt have enough BW' for some specific time resolution that is sought.

Isnt the maximum time resolution simply the reciprocal of your sampling rate?

How are those concepts related?

Dilip's points in his answer are correct. Speaking more to the context that you referenced of pulse-compression radar, I think you're getting confused by differing meanings of the oft-used word "resolution." In a broad signal processing sense, your time resolution is defined to some extent by your sample rate. But in the specific problem domain of building a radar receiver, you are concerned with being able to identify multiple echoes from distant objects and precisely observe their arrival time. "Resolution" in this context refers to resolving and separating multiple received echoes so that they can be processed independently.

A typical radar receiver uses a sliding cross-correlator to locate echoes from objects that have reflected the transmitted radar signal. The receiver knows the format of the transmitted pulse, so cross-correlation between the RF receiver's output and the transmitted pulse waveform is the optimum scheme for detecting the presence of reflected pulses in AWGN. The correlator output will contain copies of the transmitted pulse waveform's autocorrelation function (which typically has a sinc-like shape) for each received echo, shifted in time based on the range to the target that caused the reflection. In order to discriminate between the targets, their corresponding lobes at the correlator output must be sufficiently separated in time.

A "high-resolution" radar is able to finely discriminate between multiple targets in the range dimension. If your radar has multiple targets at approximately the same range, then their echoes will reach the receiver at nearly the same time. Therefore, their autocorrelation lobes will appear at the output of the correlator at nearly the same time. The ability of the radar to discriminate between the echoes depends upon the time duration of the waveform's autocorrelation lobes; a narrower autocorrelation function (ideally one that looks like an impulse) is better.

This long-winded introduction brings us to the idea of pulse compression. Pulse-compression radar waveforms are typically implemented using linear frequency modulation (also known as "chirping"); instead of transmitting a constant-frequency pulse, the transmitted frequency is swept linearly over the course of the pulse. In practice, the sweep can be done over tens or even hundreds of MHz of spectrum. What's the benefit? An autocorrelation function with nice properties:

$$
<s_{c'}, s_{c'}>(t) = T \Lambda \left(\frac{t}{T} \right) \mathrm{sinc} \left[ \pi \Delta f t \Lambda \left( \frac{t}{T}\right) \right] e^{2 i \pi f_0 t}
$$

The equation above is borrowed from the Wikipedia article; I'll defer the full explanation to that source. What's important here is the $\Delta f$ term; it refers to the amount of frequency covered by the linear frequency chirp. Since $\Delta f$ is a factor in the argument of the $\mathrm{sinc}$ function, it can easily be seen that by chirping over a larger bandwidth, the main lobe of the pulse's autocorrelation function will be narrower. Narrower lobes are more easily discriminated by the radar receiver, giving such a radar high "resolution" in terms of differentiating between similarly-ranged targets.

Just to wrap up a bit, this sort of finding should make sense. Recall that a wide-sense stationary signal's power spectral density can be defined as the Fourier transform of its autocorrelation function. The ideal autocorrelation function for a radar pulse would be an impulse; separating a bunch of echoes with "zero" width is easier than separating a bunch of more broad lobes. The Fourier transform of an impulse has infinite frequency extent. Qualitatively, it follows that autocorrelation functions with very short time extent would be comparatively wideband in the frequency domain. This is the basis for the oft-used rule of thumb in detection and estimation theory that you need a high-bandwidth signal in order to make high-resolution time-of-arrival measurements.


Self studying, getting a quality spectrogram So I've been self studying signal processing, and been trying to use Python to get quality spectrograms, quality spectrogram like the one that would be produced with audacity etc. The steps I have done are:

Me saying 'computer', was recorded on my mac in .m4a, then I used ffmpeg converting it into .wav

$ ffmpeg -i computer.m4a -acodec pcm_s16le -ac 1 -ar 16000 computer.wav
Imported into Python via:

sample_rate, samples = wavfile.read('/path/to/computer.wav')
Now I just do a quick plot of the samples:

plt.plot(samples)

and I get this kinda graph:

Why are my peaks capped?
Anyway, I carry on and as I understand it, I perform a Fourier transform that will take my signal, decompose it into the sine waves of different amplitudes and frequencies, and then get the real part of the complex numbers that are generated:samples_fft = np.fft.rfft(samples)
frequencies = np.abs(samples_fft)
I now plot a graph of frequencies:

This is what I guessed would be like, as vocal frequencies are typically around the range of 100 to a thousand, right?
A spectrogram is when I divide my samples into a certain amount of windows, and then take the Fourier transform of each of the windows, resulting in a time frequency graph.
last_wind = 0
windows = np.linspace(0,len(samples),1000)
spectrogram = []
for window in windows:
    if window != 0:
        samples_in_range = samples[last_wind:int(window)]
        samples_in_range_fft = np.fft.rfft(samples_in_range)
        frequencies_in_range = np.abs(samples_in_range_fft)
        spectrogram.append(list(frequencies_in_range))
        last_wind = int(window)
spectrogram2 = np.asarray(spectrogram).T
Plot the spectrogram:audio_length = len(samples)/sample_rate
audio_time_intervals = np.linspace(0,audio_length,999)
tempfreq = np.linspace(0,1,19)
plt.pcolormesh(audio_time_intervals, tempfreq, spectrogram2)

So now I have a spectrogram, but of very poor quality. but I've been watching youtube videos such as this one, how do I get better resolution and clarity so as to read stuff such as formants in vowels, i tried changing the size of my partitions, but still get very similar results, more windows get less data of frequencies in each window, and less windows get poor time resolution.

I am very new to signal processing, so any help would be greatly appreciated, such as a point in the right direction or pointers as to what I'm doing wrong, in order to get clear spectrograms so as to study.

Note, I have tried python packages such signal.spectrogram, but the results are very poor aswell. Plus I'd rather get to know the nitty gritties so I really understand what is going on.

Why are my peaks capped?

Your amplification gain is set to too high, or you are too close to the microphone. The amplifier is driven to its limits and it clips the output. Keep this recording and make another one where you are a little bit further away from the microphone to later compare the differences in the spectrogram. It will be interesting to see how this clipping manifests itself in the frequency spectrum.

  ...how do I get better resolution and clarity so as to read stuff such as formants in vowels, i tried changing the size of my partitions, but still get very similar results, more windows get less data of frequencies in each window, and less windows get poor time resolution.

On line 7 of your point #6, you make use of numpy.fft.rfft.

From the documentation ...the length of the transformed axis of the output is therefore n//2 + 1..

Your signal is about 2 seconds long at a sampling frequency (Fs) of 16kHz and you are dividing it into windows of 1000 samples that overlap by 965 samples. That's about 97% overlap.

This 97% overlap is a very high temporal resolution.  What you are lacking is in frequency resolution and the only way to increase this is to increase your window size. This will give the Fast Fourier Transform (FFT) the opportunity to resolve more frequencies.

These are Audacity's parameters for deriving spectrograms:

There are two things to notice here:

All of Audacity's window sizes are powers of 2. This is done to exploit the structure of the calculations in computing the FT and do it F(ast).
The window size is inversely proportional to the resolution of the FFT. A size 8 window is marked as "most wideband" and a size 32768 is marked as "most narrowband". Of course, your window size is not independent of your Sampling Frequency and there are limits to how many frequency components it can resolve given the data that has already been captured at sampling. Therefore, higher true frequency resolution means capturing more data which means increasing your sampling frequency. With 16kHz you can resolve (theoretically) up to 8kHz.

So with the exception of a few details and optimisations here and there, you are down the right track with your code.

You are missing:

The windowing functionality (nothing to do with the length of the sliding window you apply in your code (1000)). This will become important later on if you want to resolve things like formants.
"Proper" plotting which sets the x and y axes of your plot to the time and frequency scales required to interpret the spectrogram. This is not something incredibly difficult to do, but you will find yourself needing it, otherwise you will be doing "back of the envelope" calculations every time you want to read frequencies and offsets from your plots.

Ultimately, you will add those to your code too and you will have replicated, to a large extent, the functionality of scipy.signal.spectrogram.

Therefore, I would recommend to you to focus your studying efforts on the specifics of the spectrogram in order to understand how it works and what the impact of each one of its parameters is to the quality of its output.

A good introduction to the spectrogram itself is this one. And of course, it is not the only one out there.

scipy.signal.spectrogram includes all of the key parameters and it will allow you to explore spectrograms comfortably.

Finally, if you are interested in audio analysis, there are some really good libraries in the Python space, that implement a lot of functionality you will need. Some examples are pyAudioAnalysis and librosa (and here for something closer to what you are asking). Many more pointers available here.

Hope this helps.


From lowpass to bandpass I have to do two frequency shifts:

the first: to transform a lowpass filter into a highpass filter
I've implemented it in this way:

h_hpf=h_lpf.*exp(j*pi*(0:N))

the second: to transform a lowpass filter into a bandpass filter centered at $\pi/2$.
I've implemented it in this way:

h_bpf=h_lpf.*exp(j*pi/2*(0:N))
I used the fvtool() command to show frequency responses, and I've observed that the frequency response of h_hpf is also replicated in the negative frequency axis, while it doesn't happen for the frequency response of h_bpf.
Indeed, the frequency response of h_bpf appears only in the positive frequency axis, and the replica in the negative axis doesn't appear.
Moreover, its impulse response h_hbf has complex values, so I think that the correct shift is:
h_bpf=h_lpf.*(exp(j*pi/2*(0:N))+exp(-j*pi/2*(0:N)))
Which is the correct shift?

You just shifted the low pass filter to the right, so you generated a complex-valued filter, as you've observed. Multiplying a real-valued impulse response with a complex exponential naturally results in a complex-valued impulse response. What you actually have to do is shift the spectrum to the right and to the left:

$$h_{BP}[n]=h_{LP}[n]e^{jn\omega_0}+h_{LP}[n]e^{-jn\omega_0}\tag{1}$$

where $\omega_0$ is the desired center frequency. Eq. $(1)$ is equivalent to

$$h_{BP}[n]=2h_{LP}[n]\cos(n\omega_0)\tag{2}$$

from which it is obvious that the resulting band pass impulse response is real-valued.

The reason why it works with the high pass filter is that $e^{jn\pi}=\cos(n\pi)=(-1)^n$, so for this special case ($\omega_0=\pi$) it makes no difference.


What exactly mean X bit depth? I am not sure if I understand meaning of “bit depth” in audio.

As I see in many audio applications the audio sample values (when not distorted) are between minValue = -1.0f; and maxValue = 1.0f;.

I need to know the smallest possible value next to 0.0f of one sample.

As I suppose “x bit depth” means that one sample can be expressed by 2^16 various values. Am I right?

So does it mean that the smallest possible sample value change is (maxValue - minValue) / (pow(2.0f, x)) ??? So for 16 bit depth is it 2.0f / pow(2.0f, 16.0f) ?? Which is 4.6566128742e-10 ?

I have mess in mind, couse I have even problem with simple float variable. As I know it has 32 bits. So the biggest value should be 2^31 (and one bit for sign) which is 2147483647. And for int it works, but how is it possible that float max value can be 3.40282e+38 ? Which is much more than 2147483647.

So let’s imagine I have audio with 32 bit depth. So then the smallest possible sample value change is 2.0f / 2147483647.0f, or maybe 2.0f / 3.40282e+38 ???

And other issue is: should I realy devide (maxValue - minValue) / (pow(2.0f, x)) - which means bit depth concern to amplitude. Or it concern to gain - which means I should calculete it by 1.0f / (pow(2.0f, x)) ?

Please help me to understand it because my math calculations (a lot of logarithms and powers) gives me a lot of “NaN” and “inf”. So I need to narrow the possible inputs and outputs to avoid such errors. But first I need to understand how to do that.

For any help great thanks in advance.
Best Regards

Bit depth is typically applied to a fixed point representation. 16-bit means the number can go from -2^15 to 2^15-1 so from -32768 to +32767. For 8 bits that would be -128 to +127.

If you normalize to a maximum of 1, your least significant bit (LSB) value is simply the inverse, 1/32678 = 3.0518e-05 for 16 bit or 1/128 = 0.0078125 for 8 bit

Floating point works completely differently. A 32-bit floating point number has a 24 bit mantissa and an 8 bit exponent. So it's basically a 24-bit integer number combined with an 8 bit scale factor that gets "optimized" and recalculated on every floating point operation. See https://en.wikipedia.org/wiki/Floating-point_arithmetic


Normalization factor in the convolution theorem Maybe it's a trivial question, but I couldn't find any explanation for this.
According to the convolution theorem, in the continues case we add normalization factor, i.e.
$$
\mathcal F\left\{h\star g\right\} = \sqrt{2\pi}\cdot \mathcal F\left\{h\right\}\cdot \mathcal F\left\{g\right\}
$$
Do we still need to add a normalization factor also in the discrete case? If we do, what should its value be?

The choice of the normalization factor is just a matter of convention. Note that the specific correspondence between convolution in the time domain and multiplication in the frequency domain with a scaling of $\sqrt{2\pi}$, as shown in your question, applies only to the unitary definition of the Fourier transform with angular frequency as the independent variable in the frequency domain:
$$F(\omega)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}f(t)e^{-j\omega t}d\omega\tag{1}$$
If we use the non-unitary definition (as is common in the field of signal processing)
$$F(\omega)=\int_{-\infty}^{\infty}f(t)e^{-j\omega t}d\omega\tag{2}$$
then the correspondence is simply
$$\mathcal{F}\{(f\star g)(t)\}=F(\omega)G(\omega)\tag{3}$$
The same is true if we use frequency instead of angular frequency.
For the discrete-time Fourier transform (DTFT), we have exactly the same choices.
For the discrete Fourier transform (DFT), where we deal with samples in both domains, the same holds for the normalization factor $N$, where $N$ is the DFT length. In the field of signal processing, the DFT is commonly defined as a non-unitary transform:
$$X[k]=\sum_{n=0}^{N-1}x[n]e^{-j2\pi nk/N}\tag{4}$$
The inverse transform is then given by
$$x[n]=\frac{1}{N}\sum_{k=0}^{N-1}X[k]e^{j2\pi nk/N}\tag{5}$$
The unitary DFT definition is
$$X[k]=\frac{1}{\sqrt{N}}\sum_{n=0}^{N-1}x[n]e^{-j2\pi nk/N}\tag{6}$$
and the corresponding inverse transform is
$$x[n]=\frac{1}{\sqrt{N}}\sum_{k=0}^{N-1}X[k]e^{j2\pi nk/N}\tag{7}$$


How to perform sound synthesis? I am generating a simple musical instrument using a ST board. basically I have a sensor which detects motion or spatial angle and generates sounds of differnt pitch and volume depending on the angle. I already have the driver for the audio codec so all I need to do is to generate sound samples and feed the samples to it. Now I am able to get the angle readings from the sensor but the challenging part is how to generate sound. From google search so far, I think I need to generate a sin function with frequency and amplitude(volume) set according to the reading from the sensor. But my codec assumes a sampling rate 0f 48KHz so how would I go about generating sinusoids with different frequencies for a fixed sampling rate?

So far I have done this:

samplingRate = 48000;
n = 0;//reset once there is a change in frequency

//this function is called 48000 times a second
int generateSineWave(float frequency,float volume)
{
  int temp = volume*(sin(2*pi*frequency*n) ;
  n = n + 1;
  if(n == samplingRate )
    n = 0;
  if(abs(temp) > MAXVAL)
    return ERROR_CODE;
  return temp;
}

This seems to be working (I am hearing something),but I am not sure if it's generating the right frequency sinusoid. also the sound I am hearing is not very pleasant, how would I go about generating complex tones(like the ones in piano for example)? I guess also my control variables(spatial angles) needs to be low pass filtered. But apart from that any idea on how I could generate more audibly pleasant waves?

I am a newbie in sound synthesis so any suggestion/help would be greatly appreciated.

First of all, there is something wrong in your computation of the sine wave. It should be something like:

volume * sin(2 * pi * frequency * n / SAMPLE_RATE);

At the moment, your code is generating very high frequencies and what you are hearing are their mirror images.

This is the very first thing to fix. Other problems in your code include:

Resetting the variable n every time the frequency is changed - this creates discontinuities ("clicks") in the waveform.
Resetting the variable n every time it reaches the sample rate. You are not sure the waveform has completed a series of full cycles because the frequency might not be an integer.
Directly computing the phase by multiplication rather than by integration of the frequency. This causes problems when frequency is modulated.
Computing the sin function directly, which is an expensive operation on embedded hardware, that will have to be replaced by a table lookup + interpolation.

Also, make sure that your samples are scaled appropriately with respect to the scale of your DAC (for example, if this is sent to a 12-bit DAC with unsigned codes, volume must be lesser than 2047, and you should add a 2048 offset).

If you get it working, you will have a pure tone - I wouldn't call this tone "unpleasant", but it is static and not reminiscent of any music instrument.

There are many sound synthesis techniques to try, an overview of which wouldn't fit in this post. In short:

Subtractive synthesis works by generating harmonically rich tones with simple waveforms (sawtooth, square), and then shaping them through filters. It is not really good at producing "realistic" sounds, but it is relatively easy to implement (the trickiest step is bandlimited waveform generation, read about minBLEP), the sounds it produce are expressive, and it is easy to understand how each synthesis parameter relates to each aspect of sound. Used by synthesizers in the 65-85, with a revival since 95 due to the popularity of the sounds they produce in electronic/dance music. A way of "cheating" to get subtractive synthesis sounds more realistic is to use single cycle waveforms of real instruments in place of oscillators, a technique known as wavetable synthesis.
Additive synthesis directly builds up a complex tone by summing sinusoidal partials at various frequencies/amplitudes. Given the number of partials involved and the complexity of their variations over time, it is a rather unpractical method.
FM synthesis works by interconnecting simple sine wave oscillators which modulate each other's frequency. It is surprisingly good at imitating natural sounds (esp. brass, metallophones), but it is very difficult to understand how its parameters (frequency ratios and modulation indices) relate to meaningful musical parameters (sound "brightness" for example). Had a golden era in the eighties and faded away.
Sampling uses recordings of real instruments, and adjusts the playback speed of the recording (aka sample) to emulate the different notes. With this method it is trivial to get something that is immediately recognizable as a piano or guitar or coughing Ferris Bueller - though if you wish to capture all the playing nuances of an instrument you will need an increasingly large number of recordings. Currently the dominant form of sound synthesis in many applications, driven by the lowest cost of data storage.
Physical modeling attempts to simulate the physical processes happening in musical instruments - using techniques such as finite element models or digital waveguides. This process can be computationally expensive and/or require a complex modeling and data analysis effort. The results can be very expressive in reproducing all playing nuances of an instrument. Currently quite rare, except for some very narrow application (Piano synthesis).

It is really a case of "realism/expressivity/simplicity, pick any two".

I am not sure about your goals - how many instrument sounds should your system provide? do you also want to modulate other parameters of sound besides pitch and volume with the accelerometer data? I am not sure either about your target platform (the ST MCUs go from 24 MHz/32kb flash/no FPU to 168 MHz with FPU and DSP instructions and 1M of onboard flash). Another factor to take into account is polyphony - do you want your system to be able to play several notes at the same time?

I would recommend:

Sampling if you have at least 128kb of flash and need only one sound/instrument, or if you have very limited CPU.
Subtractive synthesis if you have less flash and/or need more sounds/instruments ; or if you need to use external sensor data to easily modulate sound parameters besides pitch/volume.


Understanding the DCT even symmetry I wanted to prove to myself why the DCT is better than the DFT and heres a brief of what I understand so far.

The DTFT works with a finite sampled input signal, but the frequency response obtained is continuous and periodic. The DFT samples the DTFT since it is still not practical for real world applications. The DFT inherits this periodicity and X[k] can be shown to repeat from the definition itself. Taking the IDFT of the periodic spectrum, results in a periodic waveform, and this is why the DFT assumes x[n] is periodic. Now the reason the DCT is better because its periodicity is even and that eliminates any discontinuity.

Why does the DCT assume that x[n] is mirrored? Is this established in the derivation of the DCT beforehand? Furthermore we say that x[n] doesn't have any discontinuity because X[k] is mirrored for X[k+N], and taking the IDFT of that would result in a periodic x[n] mirrored with no discontinuity. Any answer is appreciated whether its descriptive or mathematically.

Just as the DFT assumes a periodic signal by construction, the DCT assumes an even signal by construction. We can think of taking a DFT of a non-periodic signal by extending it periodically, but restricting our "viewing window" to the original signal itself. Now, what happens when you take the DFT of an even signal? The DFT is then real (which means it can be constructed with only cosines as basis functions). So, in a similar spirit to the periodic extension, we can take an arbitrary non-periodic, non-even signal, mirror it, then extend it periodically. This results in a real, cosine-based DFT. Again, we restrict our "viewing window" to the original signal itself. Note that although we've doubled the length in the frequency domain, we only have real components. Effectively this means we store the same amount of information compared with the traditional DFT. The reason the DCT is useful in compression is that compared with the traditional DFT, there is more energy in the lower frequency components than in the higher frequency components (energy compaction), so we can compress those higher frequency components. This is true because, as you said, this extended mirrored periodic signal lacks the sharp discontinuities of the traditional DFT input and as a result is smoother.


Hashing a video file Bit of a long post, but this is a complex topic:

I want to hash of a series of video files so that they can be quickly checked for duplication. I know that any amount of data can be reduced to a hash, but I'd like to do this in a way that is:

Format independent, i.e. not a binary hash, but one based on what you see
Resolution independent - aspect ratio notwithstanding
Offset independent - this is the tricky part

So far this is what I've got:

I use handbrake to extract still images from a video stream - this makes an FLV come out the same as an MP4
I wrote some code that will blur (nearest neighbor) the image into a smaller one, e.g. 640x480 => 64x48. It then reduces the colors to a smaller palate. This eliminates minor differences in aspect ratio and small pixel differences due to different codecs/bitrates.
Some more code that will hash that image into a 128-bit byte array. These can be stored in a database and compared very efficiently

Where I'm really stuck is the offset. If video A and video B are the "same" video, but B is missing the first 2 seconds, or has additional minute/hour/etc. of footage, the whole thing is off; hashes either match or they don't. In fact, I might even have issues with my sampling rate in handbrake as I don't totally understand key-frames - I'm just taking X frames per second.

So what I'd like is some way to sample sections of the video, and store many hashes per video file, and if any hashes from A match any hashes from B, do further analysis.

I figure this must be possible because of the way SoundHound and Shazaam work on music. Does anyone have insight into this?

Extracting one single hash per video and expecting it to work for subsegment matching is unrealistic.

The approach - which is the one used in audio identification systems but could be extended to video - is to extract many hashes, one for each small segment of content (say a few seconds of video). The matching process consists in matching your "query" sequence of hashes to the sequence of hashes extracted from your video. This is quite similar to a text indexing/retrieval problem, and traditional data structure for such problems (like inverted indices) would work. You can probably build your implementation on top of an existing text indexing system like Lucene/Solr. See for example, the echoprint server code, built on Solr.

This addresses the sequence aspect. Remains the problem of generating robust hashes for images. What you are doing is probably going to still be sensitive to minor differences in compression/resolution - You'll need a very rough quantization of a small number of DCT / wavelet coefficients for it to be robust.


Correctly scaling FFT of different lengths I have been looking at using two FFTs of different lengths and displaying the output magnitudes on the same graph (a shorter length window for the higher frequencies, and a longer one for lower frequencies).

Currently, I find my impulse response by performing an FFT of size 32k on my entire measurement signal, and my reference signal.  I divide each by the scaling factor 1/sqrt(M) which I found through this post
FFT averaging with different block sizes - how to scale the amplitudes

I then divide the two signals, and perform an inverse FFT, and then scale by dividing by sqrt(M).

This gives me the impulse response of the signal. I then multiply the impulse response by my two windows (one short window, and one longer. window).

I then once again perform an FFT on each window, and again scale with the 1/sqrt(M). And then get my magnitudes using:
     i = sqrt (re/im)
    mag = 20 * log10(i/M)

I would like to use the crossover frequency set by the short window length, so that above this, values from the shorter window are used, and below this values from the longer windows are used.

I can find the crossover frequency using:
    crossover = 1 / ((1/sampleRate) * Length * 2)

Currently if I use this, there is a jump in decibels where the crossover frequency is, as something is not scaled correctly.  I am pretty certain that the issue is to do with having not used any scaling factors for my windows used.

Can anyone suggest where I am going wrong at the moment, or point me towards some applicable reading material.

The aim is that I am measuring an acoustic signal from a loudspeaker, and would like to window out the environment, to try and minimise reflections in my measurements at higher frequencies.  I need a longer window in order to have enough data to represent the lower spectrum, and there is no way to avoid the reflections in this part of the measurement, but I can get the higher part of the spectrum as clean as possible.  Giving the user the tools to filter the measurement how they desire.

Below shows what I am. experiencing.  The top graph is magnitude, middle is phase, bottom is impulse, with a half hann window used. The shorter window is drawn here, the longer window is the same start time as the short but a longer length.

Currently, I find my impulse response by performing an FFT of size 32k
  on my entire measurement signal, and my reference signal. I divide
  each by the scaling factor 1/sqrt(M)

Scaling here is unnecessary since you are going to divide the two spectra anyway. Any scale factor will cancel out.

  I then divide the two signals, and perform an inverse FFT, and then scale by dividing by sqrt(M).

Scaling here is wrong. Dividing the spectra creates a transfer function and the impulse response is simply the inverse FFT of the transfer function, regardless of FFT size (unless it creates time domain aliasing)

  I then once again perform an FFT on each window, and again scale with the 1/sqrt(M).

Wrong again. Don't scale FFTs for Transfer/Function<->Impulse response relationships

  And then get my magnitudes using: i = sqrt (re/im) mag = 20 * log10(i/M)

This makes no sense to me. You seem to divide the real part by the imaginary part. That could potentially be part of a phase calculation but has nothing to do with the magnitude. If you have $H(\omega) = a+ j \cdot b$ then your level in dB is simply $L(\omega) = 10 \cdot log_{10}(a^2+b^2)$

  I am pretty certain that the issue is to do with having not used any scaling factors for my windows used.

This has nothing to do with the windows. Stop scaling anything at all, and this should just work.


Digital Image Processing Textbook with Specific Topics I am looking for a Digital Image Processing and analysis textbook, but there are specific topics that I would like the book to cover. Among other topics, I would like it to include a discussion of:

(Orthonormal) basic functions
Mahalanobis distance
Subspaces
Image registration
Error accumulation
Background estimation

What textbook(s) would you recommend?  I imagine this lack of a direct question is the reason for the two unhelpful downvotes.

A fundamental book on image processing for electrical engineers is

Two-Dimensional Signal and Image Proccesing_Jae S. Lim

A highly recommended one, again, for electrical engineers is

Fundamentals of Digital Image Processing_Anil Jain

A hands-on book on basic practical image procesing is

Principles of Digital Image Processing_Wilhelm Burger

If specifically looking for machine vision, then have a look at

Bernd Jahne - Digital Image Processing.