Signal values we will 'miss' between sampling instances during sampling of band limited signals According to the Nyquist–Shannon sampling theorem, any continuous time signal with a bandwidth $B$ smaller than Nyquist frequency $f_N=f_s/2$ (with $f_s$ the sampling frequency), which is sampled at sampling frequency $f_s$ can be perfectly reconstructed by sinc interpolation (i.e. the Whittaker–Shannon interpolation formula).

Assume we sample an unknown, bounded in magnitude, continuous time signal with constant sampling time $T=1/f_s$ at sample instances $kT$, ($k\in\mathbb{Z}$), without sampling jitter or quantization. We add the constraint that $B=\alpha f_N$, with $0\leq\alpha\leq 1$.

What I would like to figure out is the following:
At sample instant $k$, I want to determine for each $\alpha$ a worst case fractional 'overshoot' of any continuous time signal between samples $k-1$ and $k$, that I could have had. I.e. how much the continuous time signal was higher than the highest (absolute) sampled values at sample instants $k$ and $k-1$. The continuous signal, or reconstruction (since the sinc interpolation is perfect!!), that we have 'missed' by sampling.

Example:
We set $\alpha=1$ and assume discrete time signal [1,0,1,0,1,1,0,1,0,1] (notice the double 1 near the middle, and does this signal even have $\alpha=1$?). Its sinc reconstruction (blue line) from the samples (black impulses) looks as follows (I have plotted the sincs belonging to each sample in gray):

The 'overshoot' between samples $k=0$ and $k=1$, is $\approx 0.7$ or $70\%$. So we missed a peak of value 1.7 in our original band-limited continuous time, or the 'perfectly band limited reconstructed', signal. If I would have put 3 or more consecutive 1's the overshoot would have been less (the Gibbs phenomenon is in the end much smaller). Therefore, 2 consecutive continuous samples like this is 'worst case'.

Extending the signal in both directions will make the overshoot grow:

Which shows relative overshoot of $\approx 1.1$ to a value of almost 2.1.

For any sequence length $2m$, this 'overshoot' $o(m)$ will grow indefinitely, $o(m)\propto\ln{(m)}$, which goes to $\infty$ when $m\to\infty$. This is because each sample the sincs will create constructive 'interference', and the sum of $1/\pi n$ (the contributions of all the envelopes of unit sinc's) for $n\to\infty$ does not converge.

This (I think) simlar to the following: if constantly sampling a value 0, I could also reconstruct a continuous time signal with infinite amplitude that is only sampled in the nodes at values of 0, e.g. $\sin{\pi f_s t}$. This tells me the same thing: that if I allow a signal to be at the Nyquist frequency the worst overshoot I could 'miss' is infinite.

We can now state that $o(m\to\infty)|_{\alpha=1}=\infty$. And we can reason that that $o(m\to\infty)|_{\alpha=0}=0$ (sampling a constant signal of which you know that it is band limited has a unique constant reconstruction).

What if $\alpha<1$?

If we now assume we do this same sinc interpolation, but know for sure $\alpha<1$, like $\alpha=0.5$. Then, (my gut feeling says) this effect should go down and should even remain finite (when $m\to \infty$)!. Since for any signal brick-wall limited to bandwidth $\alpha f_s/2$, we get a filter impulse response of $h(t)\propto \text{sinc}{\left(\frac{t-kT}{\alpha T}\right)}$ (right?). Therefore, signal transitions cannot be as fast as for the changing impulse train example above, and therefore contributions of each sinc function during reconstruction cannot create infinite constructive interference.

My problem:
I don't know how to proceed from here; how to form a 'proof' of the worst case overshoot I could ever have find between 2 consecutive samples, knowing that $\alpha<1$, for $any$ signal (not necessarily these unit impulse-train like examples). A given value for $\alpha$ gives me a slope $\frac{dh(t)}{dt}$ of the band-limiting convolution kernel $h(t)$, which should tell me something about how much consecutive samples need to be different, but I do not see the steps to take from there to reach any generic conclusion.

I don't have a real answer but I have the feeling that this result will help you out: Bernstein's inequality says that, if the signal $x(t)$ is bandlimited to $|f|\leq B$, then $$\left| \frac{\textrm{d}x(t)}{\textrm{d}t}\right|\leq 4\pi B \,\textrm{sup}_{\tau\in\mathbb{R}}|x(\tau)| ,\,\,t\in\mathbb{R}$$ where $\textrm{sup}$ stands for "least upper bound".

I found about this inequality in Amos Lapidoth excellent (and free in PDF format) book "A Foundation in Digital Communication". A proof can be found in M. A. Pinsky, "Introduction to Fourier Analysis and Wavelets".


What is the first derivative of Dirac delta function? Could you please help me in a simple way, what is the first derivative of a Dirac delta function?
I found this answer:

The informal answer is a positive Delta function immediately followed by a negative-going Delta function.

Could you please explain this?

If you imagine a Dirac delta impulse as the limit of a very narrow very high rectangular impulse with unit area centered at $t=0$, then it's clear that its derivative must be a positive impulse at $0^-$ (because that's where the original impulse goes from zero to a very large value), and a negative impulse at $0^+$ (where the impulse goes from a very large value back to zero).
Actually, the generalized derivative of a Dirac delta impulse $\delta(t)$, denoted by $\delta'(t)$, is a generalized function (distribution) with the following properties:
$$\begin{align}\int_{-\infty}^{\infty}\delta'(t)f(t)dt&=-\int_{-\infty}^{\infty}\delta(t)f'(t)dt=-f'(0)\tag{1}\\\delta'(t)f(t)&=f(0)\delta'(t)-f'(0)\delta(t)\tag{2}\\\int_{-\infty}^{\infty}\delta'(\tau)f(t-\tau)d\tau&=f'(t)\qquad\qquad\textrm{(convolution)}\tag{3}\end{align}$$
Property $(1)$ is basically the definition of the derivative of a distribution. The right-most equality is of course a consequence of considering the special distribution $\delta(t)$.
Property $(3)$ means that convolution with the derivative of a Dirac impulse results in the derivative of the convolved function. I.e., the distribution $\delta'(t)$ is the impulse response of an ideal differentiator.
From property $(1)$ (with $f(t)=1$) it follows that
$$\int_{-\infty}^{\infty}\delta'(t)dt=0\tag{4}$$
Furthermore, from $(1)$ it also follows that $\delta'(t)$ is odd, because for every even function $f(t)$ whose derivative exists at $t=0$ we have
$$\int_{-\infty}^{\infty}\delta'(t)f(t)dt=0$$
since for even (and differentiable) $f(t)$, $f'(0)=0$ holds.


ADC quantization +thermal noise ADC noise can affect the overall noise performance of a receiver. The noise is partially quantization noise and partially thermal noise and as a general rule it extends up to the Nyquist frequency = Fs/2 . A simplification or assumption in receiver overall noise performance is to consider the noise spectrum to be flat or white. Are there any general rules or situations in terms of sampling rate (or anything else) that make this assumption grossly wrong?

Treating ADC quantization noise as a uniform white noise source is generally valid.  Please refer to this post where I further detail practical limitations on using white noise assumptions for quantization noise as well as oversampling considerations: What are advantages of having higher sampling rate of a signal?
It is typically good design practice to place the ADC quantization noise to be at least 10 dB below the amplified thermal noise of the input signal to minimize the ADC quantization impact on overall noise figure (As the waveform we wish to measure/capture typically is the input waveform including its signal and noise components, not the ADC quantization noise!).
When choosing 10 dB for example as I noted above, the ADC contribution to increasing the noise floor (and hence the noise figure given the same signal level) using the white noise assumption for quantization noise is 0.4 dB given as:
$$10 \log_{10}(1+10^{-\text{dB}/10})$$
where dB is 10 in this case. (This is simply from adding the relative powers and converting back to dB).
The trade off here is dynamic range vs noise figure degradation. Certainly you can place the quantization noise to be much lower below the input amplified noise floor to minimize overall increase in noise (by providing more gain prior to the ADC) but that would be at the expense of dynamic range. For example, if instead of placing the quantization noise 10 dB below the input amplified noise floor as we did above, we placed it 30 dB below, we would lose 20 dB of dynamic range in the process.  Typically I use numbers in the range of 10 dB to 15 dB as starting parameters for input noise floor vs ADC quantization noise, but ultimately is part of an overall noise figure budget.
This is identical to a cascaded noise figure computation for those familiar with that in a receiver design. In this case the ADC has an equivalent noise figure, thus constraining the overall gain requirements prior to the ADC in order to meet an overall noise figure objective while balancing instantaneous dynamic range constraints.
Below are my slides where I walk through a typical ADC receiver design showing how this was used in terms of a "Noise Figure" for the ADC. If there is interest in the rest of this (which shows the same considerations for the maximum input signal), please ask as an additional question specific to that and I will post that there as well (since I already have the slides):

Note: The "Noise Figure" shown in the graphic below refers to the "front-end" noise figure which is the entire receiver up to the input of the ADC. This predicts the amplified noise floor level since that would be the receiver gain up to this point plus the receiver noise figure. This is not the "ADC Noise Figure" that I reference in the text and prior graphics.


When to use cross spectral density versus cross correlation? All my life, I've used cross-correlation to analyze the strength of the relationship between an input signal and the output signal.  I use this to determine which input signal most closely influences the output signal:

Yesterday, a friend told me that he uses cross spectral density (same as coherence, he say).

To my interpretation, they seem to quantify the same thing.

When would someone use one over the other?

The cross-spectral density is in the frequency domain while the cross-correlation function is in the time domain. The two are Fourier Transform pairs, the FT of the cross correlation function is the cross-spectral density.  The the two provide the same information, just that one is in the time domain and the other is in the frequency domain.
This is just as the relationship between the auto-correlation function and power spectral density (power spectrum) for a signal which are also related as Fourier Transform pairs (as given by the Wiener- Khinchin theorem).
One common application of the cross-spectral density is finding the transfer function in a noisy system. We can cross correlate the input and output and assuming the signal is correlated but the noise is not (between input and output) we can optimally estimate the magnitude and phase response of the system.
One common application for cross-corrrelation is receiving GPS signals which at the antenna are up to 20 dB below the system noise floor! We cross correlate the extremely noisy signal received to the known pseudo-random noise sequences transmitted by the satellite and once correlated can extract the data messages and determine range to the satellites (and ultimately our position and time by triangulating on many satellites).
What you may also see is that we can do either approach in time or frequency- in the GPS example you can do the cross correlation I describe with FFT's given the relationship:
$$R_{xy}(\tau) = \text{ifft}(\text{fft}(x)\text{fft}^*(y))$$
Where:
$R_{xy}(\tau)$: Circular cross-correlation function.
$\text{fft}(x)$: The DFT of x.
$\text{fft}^*(y)$: The complex conjugate of the DFT of y.
$\text{ifft()}$: The inverse DFT
For an example with regards to GPS acquisition, see this MATLAB code I have posted on Mathwork's file exchange site.


fmcw radar sampling requirements I have been reading about "all-digital" or "software-defined radars". Lately, I was trying to simulate an FMCW short range radar in Matlab. One thing that struck me during the simulation is that to achieve a range resolution of around a 100 cm or less, the signal bandwidth is quite large -- in multiple GHz. This means, if one wishes to realise a digital implementation, one needs a really high speed ADC - multiple GigaSamples/second. Such ADCs are either not available and if you find something towards the top end of the high speed ADCs out there, the processing requirements are going to be ridiculously high. And of course, if you are clocking your baseband processor in GigaHertz, the power requirements of the device (typically an FPGA) are again going to be impractical.

My questions:

How, if at all, high-resolution, high-bandwidth digital radars are realised in practice?
what would make sense in case of FMCW radar is to mix the transmitted signal and the target return signal and digitize only the beat frequency for Inverse FFT and further analysis. In that case, could you still call it a digital or software defined radar?

A quick scribble on paper shows that FMCW radar has a resolution of

$$
\Delta d = \frac{c}{2b}
$$ with $c$ being the speed of wave propagation (aka speed of light, essentially), $b$ being the bandwidth, and $2$ coming into play because the wave has to trave both ways. For your $\Delta d\overset != 1\mathrm m$,

$$
1\mathrm{m} = \frac{3 \cdot 10^8 \frac{\mathrm{m}}{\mathrm{s}}}{2b}
$$

which means you'd only need $b\ge 150\, \mathrm{MHz}$ in my calculation.

That's something that software defined radio devices can do nowadays.

Hence, to answer your questions:
My questions:

  How, if at all, high-resolution, high-bandwidth digital radars are realised in practice?

Simply be mixing the signal down to baseband (or a low IF), and then sampling it with an ADC.

  what would make sense in case of FMCW radar is to mix the transmitted signal and the target return signal and digitize only the beat frequency for Inverse FFT and further analysis.

Yes, that's the way it's done in most "small" FMCW radars, for example automotive radars.

  In that case, could you still call it a digital or software defined radar?

Puh. Digital: yes.
Software defined: well, if everyone asked me, no, but paper doesn't blush, so there's probably someone out there who calls this kind of radar "software defined". Please don't do that.


Can Kalman filter remove 50 Hz noise? I'm dealing with a brain stimulation.

My EEG contains large spikes (>3 mV) when stimulation occurs. I want to remove the 50 Hz noise, but I can't use notch filter because it creates large artifacts.

Do you think that Kalman Filter could remove the 50 Hz noise?

Thanks for your comments. This is my registration

If i apply a notch filter i will have ringing artifacts due to the large spikes.If so I can't use the signal anymore. So my point is : i have to remove the 50hz component without using a notch filter.( I cannot delete the spikes with interpolation) My idea is to estimate the signal rappresenting the 50 hz component  with the kalman filter and then subtracting from the original signal

With your edit it becomes clear that you've modeled your problem incorrectly:

While the offending signal appears shortly with a frequency of 50 Hz, that is by no means the frequency content of the interference!

(also, your filter isn't well-designed, probably too short, judging from the impulse response it displays, to filter out 50 Hz)

You'll find that interference to look more or less like dirac impulses – and these have a very high bandwidth! So, in spectral domain (plot a PSD!) you'll notice that there's not only a peak at 50 Hz, but wideband interference.

Technically, you seem to have a comb of such, and that would actually lead to a line spectrum again, which you could "kill" with a multi-notch filter.
But then again, these aren't really diracs, and you'd need a very long filter to even come close to being able to notch these out.

I think a time-domain solution is much more promising than filtering:

Simply build a threshold that detects when a burst interferer starts and ends, and zero out all samples in between.

If you want to be really elegant, you could try to use for example linear prediction to fill these gaps instead of zeroing, but my guess is that you'd only win a tiny degree of signal quality.


How do I generate sound using MATLAB? I want to generate and play sound for tones with frequency 100 Hz and 200 Hz for 10 sec and 5 sec respectively using MATLAB with sampling frequency equal to 20500 Hz.

Also I want to plot their frequency spectrum. Is it possible to do this using MATLAB?

Something like that should do. Please keep in mind that this code can be further improved by windowing your signal, using decibel scale for spectrum, taking only first half of a spectrum, etc. This is all cosmetic stuff and you can do it on your own (or maybe I will edit this answer at some point).

%% Time domain parameters
fs = 20500;    % Sampling frequency
dt = 1/fs;     % Time resolution
T = 10;        % Signal duration
t = 0:dt:T-dt; % Total duration
N = length(t); % Number of time samples

%% Signal generation
f0_1 = 100; % fundamental frequency of first sinusoid
f0_2 = 200; % fundamental frequency of second sinusoid
x1 = sin(2*pi*f0_1*t); % first sinusoid
x2 = sin(2*pi*f0_2*t); % second sinusoid

% We want 200 Hz signal to last for half of a time, therefore zero-out
% second half - use the logical indexing of time variable
x2(t>5)=0;

% Now add two signals
x = x1+x2;

% Calculate spectrum
X = abs(fft(x))/N;
% Prepare frequency axis
freqs = (0:N-1).*(fs/N);

%% Plotting time and frequency domain
% Time domain signal plot
subplot(211)
plot(t, x)
grid on
xlabel('Time [s]')
ylabel('Amplitude')
title('Time domain signal')

% Spectrum plot
subplot(212)
plot(freqs, X)
grid on
xlim([0 max(freqs)])
xlabel('Frequency [Hz]')
ylabel('Amplitude')
title('Spectrum')

%% Playing back signal
% Normalize the audio:
x = 0.99*x/max(abs(x));

% For MATLAB R2014a use audioplayer
player = audioplayer(x, fs);
play(player)

% For older versions use wavplay
% wavplay(x, fs);

And apart from playing back it will also produce:


CWT Disapointing Frequency Separation I'm attempting to perform multi-resolution analysis via Continuous-Wavelet Transform (CWT) using Pywavelets. I've heard that CWT is supposed to be superior to STFT due to varying frequency content as a function of the time-window.

My test signal is two sinusoids of 1Hz and 5Hz, each lasting 10 seconds (see picture): f=np.sin(2.*np.pi*t)*((t>=10)&(t<=20))+np.sin(2*np.pi*5*t)*((t>=30)&(t<=40)). The sampling period is 20Hz.

Using Pywavelets, I perform the CWT as follows with the resulting spectrogram:

scales = np.arange(0.6,65,step=0.2)
coef, freqs=pywt.cwt(f,scales,'cgau1', sampling_period=dT)

As you can see the frequency resolution is quite lousy, and the peak (complex magnitude) doesn't even seem to line up at 5Hz for the second segment.

In contrast, a STFT using the Gaussian window with a standard deviation of 5 results in much sharper frequency resolution (at the expense of time sharpness):

Am I doing something here? I'm willing to sacrifice time resolution but I do need to sharpen up the frequency.

The CWT is not, per se, superior to STFT. Due to its variations in scale, it can be better, for instance, when:

you address natural signals where transients are shorter than more stationary parts,
you do not know the appropriate scale of observation.

which is neither the case from your signal: precise frequency on the same support.

On the one hand, you can work with a continuous wavelet with more oscillations, and better refine the scales and voices.

On the other hand, there are so-called STFT with windows whose width vary with frequency.

Add-ons like hybrid time-scale/time-frequency transformations, reassignement, etc. may be used, but choosing the appropriate tool requires to refine your goal.


Load cell noise filteration I'm currently working in a project as apart of a student team.
Setup
A platform is attached (placed over and screwed) two 1 Kg load cells in order to measure the weight of the system. A vibration motor is also a part of the above mentioned system.
Problem
The load cells, coupled with the entire setup causes noise while the vibration motor is switched on, preventing me to take real time data from the load cells, which is something I require.
What I think
The mechanical stresses get amplified by the vibration motor, giving me erroneous readings. (because the load cells itself don't cause too much problems when tested separately)
Are there any filtering techniques I could use in order to be able to get better real time data from the load cells?

The mechanical stresses get amplified by the vibration motor, giving me erroneous readings.

Nothing to do with amplification. Basic high school-level physics!

"Weight" is just the force that an object exerts due to gravity acting on it. Your load cells hence simply measure the force pressing down on them, and, that's the whole idea of vibration motors, that is the sum of the gravity pull and the acceleration applied to the vibrating mass (see: Newton's laws: when your vibration motor moves some mass upwards, it's getting pushed down by a force identical to the lifting force, and vice versa).

  Are there any filtering techniques I could use in order to be able to get better real time data from the load cells?

Any kind of sufficiently low-cutoff low pass filter will do, because no matter how much you vibrate, the mass of the system doesn't change, and what goes up needs to come down.

You'd typically either simply do this with a mechanically dampening system, or with a simple RC low pass operating on your analog electrical signal, which you need for anti-aliasing anyway, because otherwise you're exposing your digital signal to aliases.


How to choose sampling rate for non-band limited signal given as an expression? I have a signal given as an expression $x(t)$.  In frequency domain most of the information is around 0 and higher frequencies progressively get smaller and smaller compared to DC component.

I was wondering how to determine the sampling rate I should use if I wanted to compute discrete samples and pass them through some filter?

I was considering that if I pick some $F_s$ sampling rate, aliasing will cause frequencies above $|F_s/2|$ to be folded back into the region $F_s/2$.  Using Parseval's theorem the sum of square errors because of this appears to be $\int_{|\Omega > F_s/2|}|X(j\Omega)|^2d\Omega$ .  So I could try to figure out the acceptable error and pick $F_s$ based on that.  But that integral goes to infinity and what if the integral doesn't converge, or it's not something I can compute?

I've also seen discussions on anti-aliasing filters that mention the stop band attenuation should be higher than the SNR (eg anti-aliasing filter).  In my case I don't really have an anti-aliasing filter because I'm just computing x(t).  But it seems if I want 16 bits resolution, SNR=98dB so if I compute the frequency at which $20log(\frac{|X(j\Omega)|}{|X(0)|})<-98$ it will be as if I used anti-aliasing filter.  However, I don't understand why the stop-band attenuation is chosen to be equal to SNR, and what that means in terms of signal error.  Because while -98dB is really small compared to the rest of the signal, there are infinitely many aliased frequencies being folded and added up.

Thanks for any clarification.

Edit

I have frequency vs relative magnitude (relative to magnitude at 0 frequency) plot calculated by sampling at 100kHz and 1MHz and taking DFT of the samples.

What reasonable conclusion about required sampling rate can I make based on something like this?  How can I justify it?

For example at sampling rate 100kHz, the smallest frequency component appears to be 158 dB smaller than the largest (or 10^(7.9) times smaller).  Is that small enough to be unimportant to overall signal?  How should I go about determining what's small enough?

It depends on the spectral distribution of the signal under concern. If, for example, as you said the signal has most of its energy around zero frequency and gets progressively smaller as frequency increases then it may be possible to sample it with an acceptable error, provided that the progress is rapid enough. So it depends on the progress: Provided that it decays fast enough, you may find an acceptable sample rate which would provide a bounded error. After all many practical signals will have nonzero (but very small) energy above their half the sampling rate even after aa filtering, but nevertheless they are sampled with sufficient accuracy.

And your bit resolution vs stop band attenuation concern. That is the noise floor of the N-bit representation. Any error less than noise floor, and being uncorrelated with the signal, may be considered ok.

Finaly please note that, when the aliased frequencies fold up (say at 0 hz) for every next shifted spectra to be added, the folded tail will be progressiveley smaller, as it comes from progressiveley farther away from its center. Hence the sum of such infinitely many decaying terms may converge. Provided the folded terms are getting smaller and smaller.


How can I obtain the transfer function from these two transfer functions in this case? Using Scipy signal in Python I want to obtain a transfer function for an input which passes first from sys1 and then sys2. So explicitly
$$\text{input}{\longrightarrow}\boxed{\tt sys1}{\longrightarrow}\boxed{\tt sys2}{\longrightarrow}\text{output}$$
sys1 = signal.lti([a], [b, c])

sys1 = signal.lti([x], [y, z])

How can I cascade sys1 and sys2 programatically to obtain the overall transfer function (call it sys_total) between the input and output?

The numerator and denominator of the total system are respectively the convolution of the numerators and the denominators of the filters to be cascaded.
$$
b(n)=b_1(n)*b_2(n)
$$
$$
a(n)=a_1(n)*a_2(n)
$$
You can validate it by the following code in MATLAB
b = conv(b1, b2);
a = conv(a1, a2);

sys1 = dfilt.df2t(b1, a1);
sys2 = dfilt.df2t(b2, a2);
sys = cascade(sys1, sys2);
[b, a] = tf(sys); % convert system object to coefficients


Estimate the wavelength of an oscillatory signal I need to estimate the dominant wavelength of a signal that is looking like this:
My first attempt was using scipy.signal.argrelmin in Python in order to detect the minima of the signal, but it detects also very small local minima. Is there any useful method for this that I should look into? (The wavelength that I would like to measure is the one between the pronounced minima in the above figure)
The data for this solution can be downloaded from here.

Consider the code snippet below to find the three peaks in the graph. Note, that peak search usually only searches for positive peaks. Hence, I negated the signal to search for the down-peaks. Also, I adapt the SNR threshold in the code to detect only the very strong peaks.

import urllib2
import scipy.signal as signal

url = 'https://raw.githubusercontent.com/Omer80/wavelength/master/oscillations.dat'
response = urllib2.urlopen(url)
cr = np.loadtxt(response)
X = cr[0,:]
Y = cr[1,:]

plt.subplot(2,2,1)
plt.plot(Y)

plt.subplot(2,2,2)
plt.plot(X, Y)

peakind = signal.find_peaks_cwt(-Y, np.array([100]), min_snr=5)
peakind = np.array(peakind)
print peakind

plt.stem(X[peakind], 1+0*peakind);

period = int(np.mean(np.diff(peakind)))
print "Period between peaks:", X[period]-X[0]

Output:

[1919 2500 3080]
Period between peaks: 0.145036259065


How to recover data/information from real complex baseband signal? I have recorded real complex baseband signal. I want to recover data from the same. The complex baseband signal is basically a composite CDMA signal. How can i recover that information ?
Original Tx Composite signal

Real Rx Baseband signal

I know fundamentally I just have to multiply the composite CDMA signal with the relevant PN code in-order to get the original signal back. code length is 1024.
p_recovered = np.array([], dtype = float)
for i in range(p_len):
    p_recovered = np.append(p_recovered, 1.0*sum(composite[i*codelength:i*codelength+codelength])/codelength)
p_recovered = np.repeat(p_recovered, codelength)

r_recovered = despread(composite, r_code, codelength)
q_recovered = despread(composite, q_code, codelength)

I did that in simulation(below), and recover the original messages back.

My question is how can i recover information from actual modem received signal ?

In an actual received signal you will also have to address the frequency, phase and time offsets between the transmitter and receiver since they run off of independent clocks that aren’t otherwise synchronized, as well as frequency offsets introduced through Doppler if the transmitter and receiver are in motion relative to each other. Further information on the details of this can be found under the topics of Carrier Recovery and Timing Recovery.
For CDMA this can be done with a course frequency search by doing the OP’s correlation repeatably at different frequency offsets and using an acquisition threshold on the correlator output to determine the rough offset. Once within the acquisition bandwidth of the correlator (approximately 1/T where T is the duration of the PN sequence in time), the complex I Q outputs of the correlator can be used to measure phase change per correlation interval (which as  $f=d \phi /dt$  is the fine frequency measurement (carrier offset).  This approach can be part of a frequency lock loop or the individual phase measurements can be part of a phase lock loop; either being the carrier recovery implementation.
With CDMA one approach to timing and carrier recovery (common in GPS receivers) is to use three correlators each using the matching PN sequence separated by half a chip such that the correlator outputs form “Early”, “Prompt”, and “Late”.  The prompt is used for carrier acquisition and tracking as described above by both sweeping frequency for carrier offset and the PN code for delay offset (at each frequency test). Once acquired, the early and late can maintain tracking (time offset) by subtracting the two (zero time error when early minus late equals zero).
I have other posts detailing this further as well as timing recovery that I will add links here. The above approach is sequential so slow but minimizes resource requirements.  One link I show also provides further details on a direct (one-shot) frequency and delay (PN offset) acquisition that would be relatively fast (can be completed on a single capture) when the processing resources allow.
See:
Early - late gate spread spectrum
High modulation index PSK - carrier recovery
QPSK constellation from baseband signal


Removing Noise from time across signal data I am new to the field of signal processing. I have sensor data (Impedance across time), that output is like this.

On the right side, I have a noisy signal portion in my data, marked by the red block. I want to remove this noise, its a hand movement artifact noise.

How can I remove this noise?

I have tried to find the frequency of the noise by applying FFT on the data file attached below but I got this result

Data file

Is the width of your 5000-amplitude pulse always the same? If so, then just discard your signal samples that occur after seven seconds. If the width of your 5000-amplitude pulse varies then try lowpass filtering (experimenting with different lowpass filter bandwidths) your signal to see if the filtered signal contains the information you desire. If the width of your 5000-amplitude pulse varies then you might try differentiating your signal and measuring the time instant $t_o$ of the large negative amplitude transition. Then just discard your original signal samples that occur after $t_o$ seconds. Good Luck.


Why does this reconstruction produce a phase shift from the original signal? After an FFT of a signal is done, it is plotted as in the image below, with the original signal is on the first subplot.

Using the magnitude and phase data of each frequency, it's reconstructed, and the produced signal is as the image below.

The reconstructed is signal is good except it's shifted 4 milliseconds earlier from the the original signal. I also has done this with another signal, but this 4 milliseconds shift is there too. Why does this happen?

Additional question: I used this for harmonic analysis of 50 Hz fundamental frequency (electric power system frequency). I had a slightly hard time determining which frequency of each harmonic (50*n Hz) frequency to use because, from the FFT result, the phase angle is θ at (50*n∓0.001) Hz but suddenly becomes θ±180° at (50*n±0.001) Hz. For example, at 149.999Hz the phase angle is 223.3° but at 150.001Hz the phase angle is suddenly 403.3°.
Is this hard time determining which angle of each frequency to use for reconstruction, normal?

By the way, this is the Matlab code I use to do the FFT and the plotting.

Fs = 1/(8e-6);                % Sampling frequency
T = 1/Fs;                     % Sample time
L = 2^21;                     % Length of signal
t = (0:L-1)*T;                % Time vector

SampleNoBallast;              % The original signal

NFFT = 2^nextpow2(L);         % Next power of 2 from length of y
Y = fft(y,NFFT)/L;            % The Fast Fourier Transform producing FFT complex
f = Fs/2*linspace(0,1,NFFT/2+1);     % Frequencies to plot
P = rad2deg(unwrap(angle(Y)));       % Phase degrees from FFT complex

% Plotting the original signal
subplot(3,1,1); plot(t(1:12500),y(1:12500))
title('Arus Masukan LED T8 Opple 18 W tanpa Ballast')
ylabel('arus (mA)')
xlabel('waktu (s)')
grid on
set(gca,'ButtonDownFcn','selectmoveresize');
% Plot single-sided amplitude spectrum.

% Plotting the magnitude of each frequency
subplot(3,1,2); plot(f(1:17000),2*abs(Y(1:17000)))
title('Hasil FFT')
xlabel('Frekuensi (Hz)')
ylabel('Amplitudo (mA)')
set(gca,'ButtonDownFcn','selectmoveresize');

% Plotting the phase of each frequency
subplot(3,1,3); plot(f(1:17000),P(1:17000))
title('Sudut Komponen Harmonik')
xlabel('Frekuensi (Hz)')
ylabel('Sudut (derajat)')
set(gca,'ButtonDownFcn','selectmoveresize');

The shift is due to using an FFT with a different length than the length of the data, and likely using a non-symmetric arrangement of zero-padding to increase that original length to the zero-padded length.

Non-symmetric zero-padding rotates the phase results of an FFT, spiraling across result bins.


What's the reasoning behind changing pilot tones in OFDM? Many OFDM wireless protocols (e.g. 802.11, DVB-T2, LTE, etc.) utilize pilot tones to provide channel state information for the receiver. Within the transmission, the pilot tones generally assume a BPSK or QPSK modulation depending on the protocol.

Why do specifications have the pilot tones' phase follow a known psuedorandom sequence across frequency and time as opposed to utilizing the same phase or a straightforward rotation for simplicity?

Advantage of pilot pseudo-random sequences in standards (with examples are taken from LTE) :

Interference from other sources with different identities (UE ID, Cell ID, ...) are mitigated if detection by correlation operation is used (cell search for example).
Allow including identity in sequence by changing the value of sequence initialization. For example, pilot sequence in LTE is 31-length Gold code which is reinitialized at the start of each OFDM symbol with a value that depends on cell id.


Coefficients of the Fourier Series expansion of physical signals having finite power Question 2.7. Show that for all periodic physical signals that have finite power, the coefficients of the Fourier series expansion $
x_n $ tends to 0 as $ n \to \infty $

I have computed $ |x_n|^2 $ by multiplying $ x_n $ with $ x_n^*$. However, as per my calculations, $ |x_n|^2 = P_x $ i.e. power of $
x(t) $, which is a constant and independent of $n$. I know that this is fundamentally incorrect because $ |x_n| $ should vary with $n$. If anyone can guide me about where I went wrong that would be very helpful. Thanks and regards.

  MY SOLUTION:

Please refer to the solution from the half of the first page.

A finite power periodic signal will have the following property:

$$ P_x = \frac{1}{T_0} \int_{0}^{T_0} |x(t)|^2 dt  < \infty \tag{1} $$

where $P_x$ is the power averaged over a period of the signal. Then from Parseval's theorem, the total power can also be shown to be (for a real signal I'm assuming):

$$P_x = |a_0|^2+2 \sum_{n=1}^{\infty} |a_n|^2 < \infty \tag{2}$$

where the $a_n$'s are the continuous-time Fourier series coefficients for $x(t)$.

From the theory of power series (calculus), it's known that in order for (2) to converge (sum being less than infinity) a necessary condition is

$$ \lim_{n \to \infty} |a_n|^2 = 0 \tag{3}$$

(3) implies therefore that, for finite power periodic signal's CTFS coefficients $a_n$ also goes to zero as $n$ goes to infinity.


Is frequency resolution ("bin size" inverse) always directly proportional to sampling duration (window length)? I am new to DSP and signal processing in general (digital or analogue). In implementing an FFT/Goertzel I have noticed that the resolution ("bin size" inverse) is directly proportional to the sampling duration. In other words; the tighter you want the resolution, the longer you have to sample for (regardless of sample period / how quickly or slowly you sample (above the Nyquist rate obviously)). Is there any way around this (using some technique other than FFT/Goertzel) or is this an inherent scientific principle? If so then does the principle apply to analogue filters as well?
Would a Goertzel algorithm be the best solution for detecting the presence of a particular frequency, with the best resolution (tightest filter)?

For orthogonality of basis vectors, the DFT (or FFT) basis vectors all have a different but exact integer number of full  cycles in the DFT length, e.g. if the FFT result bin for one basis vector has K full cycles, the next bin up will be for K+1 full cycles.

So if you want a DFT to separate 2 frequencies of sinusoids, the closer they are in frequency, the longer a time interval it takes for the two waveforms to become separated by exactly 1 full period or cycle, the longer a DFT (of FFT, or Goertzel, or wavelet, etc.) it takes to fit that time interval for measurement, and have at least a full result bin of separation.

Even if you don’t want to wait until two frequencies of sinusoids differ by one full cycle, say there is so little noise that you can tell apart 2 the sinusoids if they differ by only 0.1 degree (pi/1800.0) of phase change within some window, the rule still applies.  The closer the 2 frequencies are that you want to tell apart or resolve, the longer you will have to wait for them to differ by that 0.1 degree of phase shift or change that you can barely detect (regression fit within some statistical error bound or distribution, etc.)

Otherwise you can’t tell them apart, and thus claim a resolution or level of accuracy that can.

And the more noise, the bigger the difference you will have to wait for (maybe more than 1 degree, maybe even up to more than 2 full periods of difference within the measurement time window).

Note that even with analog filters (or IIRs), generally the narrower the bandwidth, the higher the Q, the longer the impulse response, thus the longer the measurement delay to respond to or detect that narrower frequency range.

So the “tightest” fit requires either more time, and/or low noise (including interference and measurement error), or a very high signal-to-noise ratio (or of the assumption of such, and hope that that assumption is actually correct.)


Upsampling - What purpose does the interpolation filter have? I want to apply some nonlinear processing to a signal, namely: I want to implement a tube emulation which adds warmth/harmonic distortion to a digital audio signal. I am worried about aliasing, so I figured out the following processing sequence:

Upsampling of the sample rate from $f_s$ to $2 \cdot f_s$
Applying an interpolation lowpass filter with $f_g = 0.5 f_s$ and a passband gain of $2$
Applying the nonlinear processing to the upsampled signal
Again applying a filter with with $f_g = 0.5 f_s$
Decimation back to $fs$.

My question arises at point 2, the interpolation filter. If I google this topic, it results in fascinating concepts like Lagrange interpolation, spline interpolation and polynominal interpolation and the math behind those seem far from trivial. Moreover, I don't understand why this is even needed - what is the reason why one would utilize one of those concepts rather than just running the signal through a run-off-the-mill IIR with sufficient stopband attenuation? The only thing I can understand is why the passband gain needs to be non-unity.

Also, although not main point of my question, feel free to point out any logical flaws that might be in my processing sequence. Thank you!

In step 1 you insert a zero-valued sample between each pair of successive original samples. This is a gain 0.5 operation. For example consider 0 Hz. You can calculate the signed 0 Hz amplitude or direct current (DC) offset simply as the mean of the sample values. Inserting the zero samples halves the mean, so the gain is 0.5, which you need to compensate for by a gain of 2.

Your 2x-oversampled samples represent values of the continuous band-limited signal exactly at the times of the original samples or exactly half-way between them. As you suggest, in step 2 you can use a discrete-time filter, which can be a finite-impulse-response (FIR) or an infinite-impulse-response (IIR) filter. If you would have an arbitrary real oversampling ratio such as 3.12897684124, you couldn't use a discrete-time filter, but would have to resort to continuous-time interpolation, which can be done using piece-wise polynomial splines.


About using Walsh matrix as spreading sequence I'm asking about some details about using the Walsh matrix as spreading sequence code. For example, suppose I'm using a $4\times 4$ Walsh matrix given by

$$H = \begin{pmatrix}
1 & 1 & 1 & 1 \\
1 & -1 & 1 & -1 \\
1 & 1 & -1 & -1 \\
1 & -1 & -1 & 1
\end{pmatrix}$$
and I have data $x$ which will be spread on that Walsh matrix.

My question: what is the length of $x$ supposed to be? Must it be less than 4? Or can it be more than 4 too? In other words, suppose I have a signal $x = 10$ , it's supposed that is easy to spread the signal $x$ on the first column of $H$. The question is,  What's about if the signal $x = 101101$ , can we do the spread on the first column too? where the signal $x$ has length larger than the first column.

Yes you can, because you do spreading for each bit of data. So, regardless the length of your data, you are going to spread every bit on your code. In case if you have $x = 10$ you will have data after spreading with length 8, and if $x=101101$ you will have data  to send with length 4 x 6 = 24. Just try to understand the process of spreading (I think it's by XOR), then you do spread for each bit of data, not for all data.

Good luck


How to add more cowbell? Remember that Christopher Walken skit on Saturday Night live where he told the band “it needs more cowbell” https://en.wikipedia.org/wiki/More_cowbell ?  This site automated it http://www.morecowbell.dj/ using echonest http://developer.echonest.com/ .  That’s all fine.  Question is: “What is the actual algorithm behind the beat detection and cowbell insertion?”

The beat and onset detection algorithms used at the Echo Nest are probably variants/improvements of the techniques developed by Tristan Jehan in his Ph.D. This is not the only approach, and I would recommend you to try first:

Getting an onset detection function using spectral flux or complex amplitude.
Using this algorithm to detect beats (you can improve it using the "comb template" method described here to score tempo candidates - I believe this combination of Ellis' dynamic programming approach with a more sophisticated candidate scoring function is what is used in Sonic Visualizer's default beat detection vamp plug-in).

As for adding the cowbell, I think there are a bunch of heuristics that could work. First of all, you can use the breaking down of the song into sections (as provided by the EN analyzer or by an algorithm like this), and process each section individually - to have changes in the cowbell pattern coincide with changes into chorus/verse. For each section, you can align the detected onsets on a grid of 16 or 32 sixteenth note to get a discrete representation of the onsets (a bit similar to what you would see on a drum machine or sequencer); and then average these patterns to get a summary of the rhythmic structure of the section. If the section is too chaotic - use a cowbell hit on each 4th or 8th note... Possible refinements: tune the cowbell so that it matches the key of the song, stop it in the sections that have a very low loudness...


Is This Interpretation of Auxiliary Independent Vector Analysis (AuxIVA) Correct? I am implementing auxiliary function based Independent Vector Analysis (AuxIVA) from  Nobutaka Ono's original paper Stable and fast update rules for independent vector analysis based on auxiliary function technique. AuxIVA is also presented in other freely available papers, for example this one. It is an extension of plain Independent Vector Analysis where an auxiliary function of the original cost function is minimized.

There are four equations that make up the update rules of AuxIVA:

$$ r_k = \sqrt{\sum_{\omega=1}^{N_\omega} {\lvert w_k^h(\omega)x(\omega)\rvert}^2} $$
$$ V_k(\omega) = {1\over{N_t}}\sum_{t=1}^{N_t}{{G'(r_k) \over r_k} x(\omega) x^h(\omega) } $$
$$ w_k(\omega) = (W(\omega)V_k(\omega))^{-1}e_k$$
$$ w_k(\omega) = {w_k(\omega) \over {\sqrt {w_k^h(\omega)V_k(\omega)w_k(\omega)}}} $$

where

$ N_s $ is the number of sources.

$ k = 1,2,...,N_s$ is the source index.

$ e_k = [0,0,...0,1,0...]^T $, is a column vector with its only non-zero entry being a 1 at the k-th index.

$ G'(r_k)$ is the first derivative of the AuxIVA contrast function and since I'm using $G(r) = r$, $ G'(r_k) = 1$.

$ N_w $ is the number of frequency bins and $ N_t $ is the number of time bins in the STFT of the mixture.

The input I am using is structured like this:

$X$ is a tensor $N_s  * N_t * N_\omega $ which is the short-time Fourier transform (STFT) of $N_s$ mixtures.
$W$ is the mixing tensor that will be estimated by repeatedly applying the 4 update rules. Its dimensions are $N_s*N_s*N_\omega$

In Matlab notation, $w^h_k(\omega) = W[k,:,\omega] $ and $x(\omega)=X[:,:,\omega]$.

When I evaluate $r_k$ using the first equation , I get a vector of size $ 1 * N_t $.

To evaluate $V_k(\omega)$ using the second equation, I evaluate $ {G'(r_k) \over r_k}$, which is just $r_k$ because of the assumption $G(r_k) = r_k$. I then do point-wise multiplication between $ {G'(r_k) \over r_k} $ and $x(\omega)$ then post-multiply this product by $x^h(\omega)$. Because $x(\omega)$ has size $N_s*N_t$ and $ {G'(r_k) \over r_k} $, the point-wise multiplication between them a weighted copy of $x(\omega)$. When this weighted copy is post-multiplied by $x^h(\omega)$, a weighted covariance matrix of size $N_s*N_s$ is generated.

In the third equation I post-multiply $W(\omega)$, that has size $N_s*N_s$, by $V_k(\omega)$ from the previous equation, that has size $N_s*N_s$. I then invert this matrix to get $(W(\omega)V_k(\omega))^{-1}$ and post-multiply this by $e_k$ to get $w_k(\omega)$ that is a vector of size $1*N_s$

Finally, in the last equation, I evaluate ${\sqrt {w_k^h(\omega)V_k(\omega)w_k}}$, that is a scalar, and divide the vector $w_k(\omega)$ by this scalar to get a scaled version of $w_k(\omega)$.

I'd like to know if my interpretation of the four equations is correct. I do not have an in-depth understanding of the derivation of these equations so I'm relying on little information present in these papers to interpret and implement the algorithm.

I think your understanding is right, and it has already been implemented in many open sourced codes, like https://github.com/LCAV/pyroomacoustics and https://github.com/n-narisawa/bss

However, the performance of AUXIVA is not as good as original IVA according to my experiements. Do you have some new conclusions about the performance and the real-time implementation of AUXIVA? Thanks.


Representing a positive 0-limiting continuous time domain function, in terms of spectral samples, outside of a window Bit of an odd question

Assume you have a known function $x(t)$ which is continuous, but approaches but is not equal to, 0, on the interval $[-\infty,\infty]$ . You know also that it approaches 0 exponentially fast on some interval $[T_a, T_b]$ e.g. the Normal or Cauchy Distribution, or some pulse with characteristics you already know with the properties as described. You also know this time domain function's closed form Fourier Transform $X(f)$, which is also continuous:

$$ X(f) \triangleq \mathscr{F}\{ x(t) \} = \int\limits_{-\infty}^{\infty} x(t) \, e^{-j 2 \pi f t} \ dt $$

Say now, you wanted to reconstruct an approximate $\hat{x}(t) \approx x(t)$ from only spectral samples $f_k$ coming from the Fourier domain, and a given $t$. This is perhaps in the form of an IFT of an FT. It is acceptable that this reconstruction of $\hat{x}(t)$ be equal to $0$ outside of the given time domain range $[T_a, T_b]$ at a specified error which is determined beforehand.

I have only dabbled a little in signal processing, but my current understanding suggest that one can only represent this reconstruction $\hat{x}(t)$signal by one of two ways

Since you know the closed form of $x(t)$, you can take a Fourier series with coefficients calculated on the given interval. This will enable one to take the Fourier Transform of this FS, and then apply the Inverse Fourier Transform to gain the desired representation. By using the FS however, it means that your reconstructed signal is inherently assumed to be periodic outside of your given interval, such that if one were to provide a $t$ such that $t<T_a$ and $t>T_b$  you would not get 0, but instead a value, in the period you're in, that corresponds to a value inside $[T_a, T_b]$
Since you know the close from of $X(f)$, you can directly sample in the frequency domain for a given spectral sampling rate. This however would still give the same result as 1. because the sampling impulse train is periodic, while $X(f)$ is continuous, meaning one would need an infinite number of samples, which is definitely not desirable!

The possible idea I had for a solution is to say that $x(t)$ is multiplied in the time domain with $w(t)=\operatorname{rect}(t/T)$ for some period $T$. Assume for simplicity now that our $[T_a, T_b]$ is symmetric about 0. This would mean we would require convolution with a $\operatorname{sinc}$ function in the frequency domain. But since the $\operatorname{sinc}$ function is continuous, our convolved $W(f) \star X(f)$ would then be continuous which then means we would need again infinitely many samples!

I have a feeling something is missing in my understanding but I've searched everywhere for something that hints at an answer, but everything talks about band-limited signals, not time-limited signals, and any that regard time-limited signals only applies to a solution on a specified interval.

So specifically I am curious how to use the additional information (that we know either 1. $x(t)$ is 0 outside a given interval, or 2. we have the closed form of $x(t)$). Also note that this assumed a linear shift invariance for $t$ in $x(t)$. Perhaps an additional approximate integration step between each sample in the frequency domain is required? Or some extension with Poisson Summation?

Any advice is greatly appreciated!

edit:
- thanks Robert for the suggested clarifications

Finite $X_k$ samples will inevitably lead to a periodic approximation. The question, as is, is unfortunately and quite interestingly, unanswerable.

Indeed, under the standard formalisms, the reconstructed signal $\hat x(t)$ will be peridic:
$$
F{X(f)|||(f)}=\sum^{\inf}_{i=-\inf} X(i\Delta f)e^{i2\pi i\Delta f t}\Delta f = \hat x(t)
$$

If we truncate the signal $x(t)$ to be time-bounded, and using a similar argument than the Paley Wiener theorem to the direct Fourier Transform, we will obtain a exponentially converging spectrum on $X(f)$ when $f$ tends to infinity.

Based on the last, a logarithmic approximation could be sketched, in which the $f_i$ are being preferred when closest to zero, and sparse for increasing values:
$$
\hat x(t)=\sum^{\inf}_{f_i=-\inf} X(f_i)e^{i2\pi f_i t}\Delta f_i
$$

And under such aproximation, we could let the first terms of the sum to be again, approximated with known spectrum shapes (as $sinc(f)$ as you said). But that is not a train of frequency samples. This is a different work, and perhaps a different question?.


peak limiting/audio compression formula needed I'm looking for a formula to effectively compress an audio waveform to limit peaks.  This isn't an "automatic volume control" application where one would control amplifier gain to maintain a volume level, but rather I want to limit ("soft" truncate) individual peaks.  (I know this introduces harmonics, but I'm trying to analyze the data, not listen to it.)

My (very crude) formula so far is:

factor = (10 * average / level) + exp(-sqrt(0.1 * level / average))

Where level is the instantaneous sound level, average is the historical average sound level, and factor is a multiplier used to produce the "adjusted" level (factor times level).

Further, this multiplier is only applied if it computes to a value less than 1.  Otherwise level is left unadjusted.

The intent is to limit the adjusted level to some multiple (about 15x with this formula) of the historical average.  This formula is sorta what I need, but exhibits a "dip" as the numbers get larger.  That is, the adjusted level (ie, factor times level) increases up to a point with increasing unadjusted level but then, rather than going asymptotic, begins to actually get smaller.  (In fact, the first factor was added primarily to prevent the formula from going to zero with extremely high values.)

(The reason for wanting to limit the values this way is primarily so transient noise doesn't seriously upset the running average of the sound level.  But when you're analyzing snores "transient noise" is quite significant, so I can's simply squelch it.)

So, can anyone suggest something better?  (It seems that asymptotic behavior is easy to produce when you don't want it, but hard when you do.)

Two problems here: how to get a reliable estimate of the level, and how to compress the data.

Use robust statistics on the original (not peak-limited) data like median or quantiles instead of a running average to make your "typical level" detection robust to outliers.
$k \times \tanh(\frac{x}{k})$ works nice as a $C^\infty$ compression formula and it is actually what is happening in some audio circuits (using OTAs). To get adaptive compression that preserves the dynamics of the original signal and just remove transients, make k track the smoothed "average" level.

Blue: original signal
Green: 2 x median of the absolute value over a sliding window as a "typical level" detection
Red: tanh compression (formula given above with k equal to the level plotted in green)


Break a signal into segments, add them together and then perform an FFT on the result If you break a 4 second signal into 0.5 second segments, add them together and then perform an FFT on the resultant 0.5 second data.

Will the FFT / frequency spectrum be the simlar to an FFT on the original 4 second data?

This would mean that you can store long amounts of data in smaller memory if all you cared about was retaining the frequency components... I'm not sure if I am missing something though?

That's not true (unless the signal is periodic, and those pieces are the base periods of it).

If you divide a signal $x_L[n]$ of a length $L$ into smaller pieces $x_i[n]$ of length $N$ each,

$$ x_L[n] = x_1[n] + x_2[n-N] + x_3[n-2N] + x_4[n-3N]$$

and sum them into a shorter signal $x_N[n]$ of length N

$$ x_N[n] = x_1[n] + x_2[n] + x_3[n] + x_4[n]$$

and take $N$-point DFT/FFT of their resultant $x_N[n]$, you won't be able to reconstruct the original long singnal $x_L[n]$ back.

Indeed it's not even about FFT. You would not be able to retrieve the components $x_i[n]$ from their sum $x_N[n]$. So since there is no way to reconstruct back $x_L[n]$ in this manner.


When does the convolution of $2$ signals equal zero? For example, how can I determine if the convolution of $x(t)$ with $y(t)$ is equal to $0$?

Time-domain convolution is frequency-domain multiplication. If at all frequencies at least one of the signals is zero-valued in frequency domain, then the convolution of the two signals will be zero-valued at all frequencies, and at all times. Except for a zero signal, no finite-length signal has a continuous run of frequency domain zeros, so your choices are limited to infinite-length signals such as periodic signals with no coinciding non-zero harmonics and ideal lowpass–highpass filter impulse response pairs.


How to get finer grained spectral decomposition for audio signals? I am not a "proper" DSP guy at all, so apologies in advance if I am making any basic errors in this question. I would really appreciate it if you could point out any errors I make.

I have an audio signal sampled at (say) 44100 Hz. I want to get a spectrogram from this signal and display it. For this question, I will be using the formula as displayed in https://ccrma.stanford.edu/~jos/mdft/DFT_Definition.html .

My approach is as follows:

Take (say) 1024 of the frames at a time (which I will denote by $x$), and run an FFT routine over them, probably after applying a Hann window.
This will give me the frequency domain representation of the input window. I will denote the output of the FFT routine by $X$.

$X$ will effectively contain 512 values, with the value at index $n$ being the Fourier coefficient of the frequency $n * 22050 / 512$ Hz. Thus, say, the 300th sample of $X$ will correspond to the frequency 12919.92 Hz. (Is this understanding correct?)

Now my actual problem: I am not really interested in frequencies above 5000 Hz since my input will be human voice. Instead of $X$ covering a frequency range of 22050 Hz, I would like to "zoom in" so that it covers a range of only 5000 Hz. I.e., I want a finer grained spectral decomposition till a lower frequency of 5000 Hz. My question: how do I do this?

An approach I can think of: currently my sampling rate is 44100 Hz. If I can "reduce" the sampling rate to 11025 Hz, then the normal FFT approach will work and the output of the FFT routine will cover only frequencies till approximately 5000 Hz. I can accomplish this by either,

Averaging every 4 frames of my original input.
Taking only every 4th frame of my original input.
Taking the max of every 4 frames of my original input.

Will any of the above 3 approaches work? Is there any other way of getting a finer grained spectrum?

Edit: it struck me later that I could achieve the result by using a larger "bin" . I.e., 8192 frames instead of 1024 frames as input to my FFT routine. Is this the best approach available?

if you are not experiencing computational limits (like you have all of the CPU speed and memory that your heart desires), what i would suggest is just do a much longer DFT or FFT.  instead of 1024 points, do it with 16384 points or 65536 points or more.  you will have lotsa data up in the top two octaves that you don't need and can decide to ignore.

if you do not want to have all of that data above 5512 Hz, then what i would recommend is formally downsampling by a factor of 4 and running the FFT on the downsampled data.  to downsample legitimately, first you have to apply a good low-pass filter and kill all of the content above 5 kHz.  after that LPF, then simply pick out every 4th sample, discarding the 3 samples in between.  your data will be reduced by a factor of 4, your sample rate will be 11025 Hz and your Nyquist will be 5512.5 Hz and you can send that data into an FFT and look at those results.


Prevent negative image with BPSK modulation in a complex signal In a baseband complex signal, I need to generate a BSPK modulated signal at a given frequency.  If I generate a CW at this frequency, all is good.  However, when I then BPSK modulate the CW, I get a signal at both the positive and negative frequencies.

Is there a way of implementing BPSK modulation without this negative frequency image, or must I simply implement a filter to remove the image?

An actual BPSK signal can be expressed as
$$\begin{align}
s(t) &= \operatorname{Re}\left[\sum_{n=-\infty}^\infty
(-1)^{b_n}p(t-nT)e^{j(2\pi f_c t + \theta)}\right]\\
&= \left[\sum_{n=-\infty}^\infty
(-1)^{b_n}p(t-nT)\right]\cos(2\pi f_c t + \theta)
\end{align}$$
where

$b_n \in {0,1}$ is the $n$-th data bit (note that $(-1)^{b_n} = \pm 1$),
$p(t)$ is the baseband pulse shape. Typically $p(t)$ is nonzero only for
$t \in \left[-\frac{T}{2},\frac{T}{2}\right)$, e.g.
$\operatorname{rect}\left(\frac{t}{T}\right)$.
$f_c$ is the carrier frequency which is sometimes chosen to be a harmonic of
the bit rate $T^{-1}$
$\theta$ is the carrier phase which might be arbitrary (depending
on initial conditions in the local oscillator when it is turned on) or carefully controlled to have
value $0$ or $\pm \frac{\pi}{2}$ etc. by MATLABi people.

Notice that the term in square brackets in the second displayed equation
above is a series of nonoverlapping pulses $\pm p(t-nT)$ and that
the value of $b_n$ affects the BPSK signal phase only during
$\left[nT-\frac{T}{2},nT+\frac{T}{2}\right)$.

With that as a preliminary setting of the stage, note that the negative
frequency image comes about from taking the real part in the first
displayed equation above.  We have
$$\operatorname{Re}\left[e^{j(2\pi f_c t + \theta)}\right]
= \frac{e^{j(2\pi f_c t + \theta)} + e^{-j(2\pi f_c t + \theta)}}{2}
= \frac{e^{j(2\pi f_c t + \theta)} + e^{j(2\pi (-f_c) t - \theta)}}{2}.$$
Furthermore,

  no filtering can remove the negative frequency image.

You might say "Well, let's not take the real part at all
since that's the source of the trouble!" and say
that
$$\sum_{n=-\infty}^\infty
(-1)^{b_n}p(t-nT)e^{j(2\pi f_c t + \theta)}$$ is a BPSK signal
containing positive frequencies only. Unfortunately, that is
a complex-valued signal that cannot be transmitted on a single
wire or radiated from a single antenna.  So you need to live with
the negative frequency image: it is there to make your complex
signal more real.


Absolute convergence of periodic sinc interpolation An $N$-periodic complex discrete-time sequence $[x_0, \dots, x_{N-1}]$ can be resampled to an $M$-periodic sequence $[y_0, \dots, y_{M-1}]$ with $M>N$, using sinc interpolation:

$$\begin{align}y_m &= \sum_{n=-\infty}^\infty \operatorname{sinc}\left(\frac{Nm}{M} - n\right)x_{n\operatorname{mod}N} \\&= \sum_{k=-\infty}^{\infty}\sum_{n=0}^{N-1}\operatorname{sinc}\left(\frac{Nm}{M} - n - Nk\right)x_n\end{align}\tag{1}$$

where $\operatorname{mod}$ denotes the modulo operation and:

$$\operatorname{sinc}(x) = \begin{cases}1&\text{if }x=0,\\\frac{\sin(\pi x)}{\pi x}&\text{otherwise}.\tag{2}\end{cases}$$

Eq. 1 can be seen as resampling an $N$-periodic continuous-time signal from samples $x_n$ at times $n + Nk$ to samples $y_m$ at times $\frac{Nm}{M}$.

For example, a $2$-periodic complex discrete-time sequence $[x_0, x_1]$ can be resampled to a $4$-periodic sequence $[y_0, y_1, y_2, y_3]$:

$$\text{Eq. 1, }N=2,\,M=4$$
$$\Rightarrow\left\{\begin{align}y_0 &= x_0\\
y_1 &= \sum_{k=-\infty}^\infty\Bigg(\operatorname{sinc}\left(2k+\frac{3}{2}\right)x_0 + \operatorname{sinc}\left(2k+\frac{1}{2}\right)x_1\Bigg)\\
y_2 &= x_1\\
y_3 &= \sum_{k=-\infty}^\infty\Bigg(\operatorname{sinc}\left(2k+\frac{1}{2}\right)x_0 + \operatorname{sinc}\left(2k+\frac{3}{2}\right)x_1\Bigg)
\end{align}\right.\tag{3}$$

The two series in Eq. 3 converge conditionally, with for example these possible rearrangements of the first series that give conflicting results if $x_0 \ne x_1$:

$$\begin{gather}\sum_{k=0}^\infty\bigg(f(-k) + f(k+1)\bigg)\\= \frac{x_0 + x_1}{2},\\
\sum_{k=0}^\infty\bigg(f(-k) + f(2k+1) + f(2k+2)\bigg)\\= \frac{(x_1-x_0)\ln(2)}{2\pi} + \frac{x_0 + x_1}{2},\\
\sum_{k=0}^\infty\bigg(f(-k) + f(3k+1) + f(3k+2) + f(3k+3)\bigg)\\= \frac{(x_1-x_0)\ln(3)}{2\pi} + \frac{x_0 + x_1}{2},
\end{gather}\tag{4}$$

with shorthand $f(k) = \operatorname{sinc}\left(2k+\frac{3}{2}\right)x_0 + \operatorname{sinc}\left(2k+\frac{1}{2}\right)x_1$.

Under which condition does the series given by Eq. 1 converge absolutely?

I'm top-editing this since it answers the question directly.
The sinc series is fundamentally a $C/x$, so you can extract as many absolutely convergent series out of it as you want, but what is left over is still only conditionally convergent.  Also, you can rescale $x$ and it is still a $C/x$ series.
Saying you have a summation to or from infinity is an informality.  Formally, you have a finite sum to some value, and take the limit as that value goes to infinity.
Therefore, your first and second series should have been done like this:
$$ \lim_{L \to \infty} \sum_{k=-L}^{L} f(k) =  \lim_{L \to \infty} \left[  \sum_{k=0}^{L} f(-k) + f(k+1) \right]  $$
$$ =  \lim_{L \to \infty} \left[ \sum_{k=0}^{L} \left( f(-k) + f(2k+1) + f(2k+2) \right) + \sum_{k=0}^{L+1} f(-k-L-1) \right] $$
Likewise, your third should have added this:
$$ \sum_{k=0}^{L+1} \left( f(-k-L-1) + f(-k-2L-3) \right) $$
Sometimes it takes a while to get around to where you should have been in the first place.  I'm deleting the rest.  Whoever is curious can find it in the edit history.

Proceeding informally....
First rearrange it:
$$
\begin{aligned}
y_m &= \sum_{n=0}^{N-1} x[n] \sum_{k=-\infty}^{\infty}  \operatorname{sinc} \left( \frac{Nm}{M} - n - Nk \right) \\
&= \sum_{n=0}^{N-1} x[n] W_m[n]
\end{aligned}
$$
One way to look at that is a resampled value is a linear combination (weighted average) of the sample points.
Another way is that you now have $N$ separate infinite series, all of the form:
$$
\begin{aligned}
W_m[n] &=  \sum_{k=-\infty}^{\infty}  \operatorname{sinc} \left( \frac{Nm}{M} - n - Nk \right) \\
 &=  \sum_{k=-\infty}^{\infty}  \frac{ \sin \left( ( Nm/M - n - Nk ) \pi\right) }{ (Nm/M - n - Nk) \pi }  \\
\end{aligned}
$$
Even $N$ Case:
$$ W_m[n] =  \sin \left( ( Nm/M - n ) \pi\right) \sum_{k=-\infty}^{\infty}  \frac{ 1 }{ (Nm/M - n - Nk) \pi }  $$
Odd $N$ Case:
$$ W_m[n] =  \sin \left( ( Nm/M - n ) \pi\right) \sum_{k=-\infty}^{\infty}  \frac{ (-1)^k }{ (Nm/M - n - Nk) \pi }  $$
Clearly, both are cases of $C/x$ series and not absolutely convergent.  If $Nm/M$ is an integer all the terms are zero except perhaps the zeroth terms.
As for the second comment, if I remember correctly (and I've already proven I didn't remember well), doing it formally does away with all the rearrangement tricks.  And yes, If I remember correctly, absolutely convergent series are immune to rearrangement tricks.
This too:
A series converges if and only if the sequence of partial sums converges.
A sequence converges if and only if for any given $\epsilon$ there exist a $\delta$ so for every $k > \delta$ the absolute value of the difference of the limit and the sequence value is less than $\epsilon$.
Stamp it on your forehead for formal occasions.
Disclaimer:  Been a long time ...

As clearly as I think I can say it:
The only conditions for which the series in Olli's Eq (1) will converge absolutely is when all the terms heading towards infinity are zero, since then their absolute values are zero.  This happens when all the $x_n$ are zero (the trivial solution) or $Nm/M$ is an integer.  Both the even and odd cases under any different conditions can be rearranged to be summations of alternating monotonically decreasing sequences, therefore they converge conditionally since they diverge absolutely.

Epilogue:
There is no need to do the infinite summation at all.  Direct closed form expressions exist for the odd and then even case based on the interpolation functions found when considering an inverse DFT as a continuous function.  The derivation of the functions can be found in the epilogue of my answer here:
How to get Fourier coefficients to draw any shape using DFT?
The derivation is based on the definitions of the DFT, the inverse DFT, and a finite geometric summation.
Resampling the continuous function at $M$ evenly spaced (in the cycle domain) points can be done by a simple variable substitution.
$$ t = \frac{m}{M} 2\pi $$
The direct sample set to sample set equations are then as follows.
Odd case:
$$ y_m = \sum_{n=0}^{N-1} x[n]
\left[ \frac{   \sin \left( N \left( \frac{m}{M} - \frac{n}{N} \right) \pi \right) }
            { N \sin \left(   \left( \frac{m}{M} - \frac{n}{N} \right) \pi \right) }
\right] $$
Even case, evenly split Nyquist bin:
$$ y_m = \sum_{n=0}^{N-1} x[n]
\left[ \frac{   \sin \left( N \left( \frac{m}{M} - \frac{n}{N} \right) \pi \right) }
            { N \sin \left(   \left( \frac{m}{M} - \frac{n}{N} \right) \pi \right) }
\right]         \cos \left(   \left( \frac{m}{M} - \frac{n}{N} \right) \pi \right) $$
These are mathematically equivalent to taking the DFT of size $N$, zero padding it at the Nyquist frequency to size $M$ (splitting the Nyquist bin in the even case), then taking the inverse DFT to recover a $M$ point upsampled sequence.  All the upsampled points lie on the underlying continuous interpolation function no matter what the point count.
For the $N=2$, $M=4$ case:
$$
\begin{aligned}
y_0 &= x_0 ( 1 ) + x_1 ( 0 ) = x_0 \\
y_1 &= x_0 \left( \frac{ \sin( \pi / 2 ) }{ 2 \sin( \pi / 4 ) } \cos( \pi / 4 ) \right) + x_1 \left( \frac{ \sin( -\pi / 2 ) }{ 2 \sin( -\pi / 4 ) } \cos( -\pi / 4 ) \right) \\
&= \frac{1}{2} ( x_0 + x_1 ) \\
y_2 &= x_0 ( 0 ) + x_1 ( 1 ) = x_1 \\
y_3 &= x_0 \left( \frac{ \sin( 3 \pi / 2 ) }{ 2 \sin( 3 \pi / 4 ) } \cos( 3 \pi / 4 ) \right) + x_1 \left( \frac{ \sin( \pi / 2 ) }{ 2 \sin( \pi / 4 ) } \cos( \pi / 4 ) \right) \\
&= \frac{1}{2} ( x_0 + x_1 )
\end{aligned}
$$
Which should be the results you are expecting.
An infinite number of sinc functions can now take the day off.

Suppose that instead of doing halfsies on the Nyquist bin we apportioned them as $(1/2+g)$ and $(1/2-g)$, this would alter the continuous interpolation function as follows.
$$
\begin{aligned}
D(t_n) &= \left( \frac{1}{2} + g \right) e^{i(N/2) t_n } + \left( \frac{1}{2} - g \right) e^{i(-N/2) t_n } + \sum_{l=0}^{N-2} e^{i ( l - N/2 + 1 ) t_n } \\
&= \cos \left( \frac{N}{2} t_n \right) + i 2 g\sin \left( \frac{N}{2} t_n \right)  + \frac{  \sin( t_n  N /2 ) } { \sin( t_n / 2 ) } \cos( t_n / 2 ) - \cos( t_n  N /2 )  \\
&= \frac{  \sin( N t_n/2 ) }{ \sin( t_n / 2 ) } \cos( t_n / 2 ) + i 2g\sin \left( \frac{N}{2} t_n \right)
\end{aligned}
$$
The extra term introduced is purely imaginary.  That can be folded in, but I prefer to leave it separate when put back into the function definition.
$$
\begin{aligned}
z(t) &= \sum_{n=0}^{N-1} x[n] \left[ \frac{ \sin( N (t - \frac{n}{N}2\pi) / 2 ) } { N \tan( (t - \frac{n}{N}2\pi) / 2 ) } + i \frac{2g}{N}\sin \left( N (t - \frac{n}{N}2\pi) / 2  \right)  \right]
\end{aligned}
$$
If is obvious that any non-zero value of $g$ will add "energy" to the signal, thus the $g=0$ solution, corresponding to halfsies on the Nyquist bin, is the most natural solution, or lowest energy, out of a whole family of solutions of periodic bandlimited at $N/2$ functions.
The more significant convincer for me is that it also introduces imaginary values into what is other wise fully real set of weighting values.
Whether R B-J's series converges uniquely to this "natural" solution, or the "natural solution" is a unique solution (it is not) are two totally separate issues.

Olli, I hope this makes you smile.
Start with the discrete resampling formula for the odd $N$ case.
$$ y_m = \sum_{n=0}^{N-1} x[n]
\left[ \frac{   \sin \left( N \left( \frac{m}{M} - \frac{n}{N} \right) \pi \right) }
            { N \sin \left(   \left( \frac{m}{M} - \frac{n}{N} \right) \pi \right) }
\right] $$
Since the sequence of $N$ points is periodic $( x[n] = x[n+N] )$ and all the points are covered, we can shift the summation range to be zero centered.
$$ L = (N-1) / 2 $$
Also, the $m$th point can be located on the $n$ scale.
$$ w = m \frac{N}{M} = \frac{m}{M} N $$
Since the $M$ resampled points are evenly spaced along the cycle, they too can be arbitrarily shifted to be zero centered, though strictly that is not necessary.
Since "$t$" has already been used above, the scale of the domain of the continuous interpolation function, both will get new names.  "$z(t)$" and "$Y(\omega)$" describe the same function.  Plug all the defined values in.
$$
\begin{aligned}
y_m = Y(w) &= \sum_{n=-L}^{L} x[n]
\left[ \frac{   \sin \left(             \left( w - n \right) \pi \right) }
            { N \sin \left( \frac{1}{N} \left( w - n \right) \pi \right) }
\right] \\
&= \sum_{n=-L}^{L} x[n]
\left[
\frac{\frac{\sin \left( \left( w - n \right) \pi \right)}{ \left( \omega - n \right) \pi }}
{\frac{\sin \left( \frac{1}{N} \left( w - n \right) \pi \right)}{\frac{1}{N} \left( w - n \right) \pi }}
\right] \\
&= \sum_{n=-L}^{L} x[n]
\left[
\frac{\operatorname{sinc} \left( w - n \right) }
     {\operatorname{sinc} \left( \frac{1}{N} \left( w - n \right) \right)}
\right] \\
\end{aligned}
$$
Now it's time to take the big step, that is, big stroll out to infinity.  The cycle of $N$ points grows until one cycle spans negative to positive infinity.  As it gets bigger, the circular nature gets more remote.
$$
\begin{aligned}
\lim_{N \to \infty} y_m &= \lim_{N \to \infty}  Y(w) \\
&= \lim_{N \to \infty} \sum_{n=-L}^{L} x[n]
\left[
\frac{\operatorname{sinc} \left( w - n \right) }
     {\operatorname{sinc} \left( \frac{1}{N} \left( w - n \right) \right)}
\right] \\
&=  \sum_{n=-\infty}^{\infty} x[n]
\left[
\frac{\operatorname{sinc} \left( w - n \right) }
     {1}
\right] \\
&=  \sum_{n=-\infty}^{\infty} x[n] \operatorname{sinc} \left( w - n \right) \\
&=  \sum_{n=-\infty}^{\infty} x[n] \operatorname{sinc} \left( \frac{Nm}{M} - n \right)
\end{aligned}
$$
Now look at that.  The Whittaker–Shannon interpolation formula has been derived from scratch and we are right at your starting point.
The even case can be done similarly and ends up with the same formula.

Definition of DFT of $N$ samples
Inverse DFT used as Fourier Series Coefficients for interpolation function
Dirichlet Kernel form of interpolation function
Interpolation function used for $M$ samples
Even and Odd Discrete Weighted Average Resampling Formulas
N goes to infinity
Whittaker–Shannon emerges
Whittaker–Shannon applied to a repeating sequence of $N$
Convergence questioned

I hope realizing using step 7 to achieve what step 2 has already answered will put a smile on R B-J as well.  Your proof lies there.
For $ N = 2 $
$$
\begin{aligned}
y_m &= \sum_{n=0}^{1} x[n]
\left[ \frac{   \sin \left( 2 \left( \frac{m}{M} - \frac{n}{2} \right) \pi \right) }
            { 2 \sin \left(   \left( \frac{m}{M} - \frac{n}{2} \right) \pi \right) }
\right]         \cos \left(   \left( \frac{m}{M} - \frac{n}{2} \right) \pi \right) \\
&= \sum_{n=0}^{1} x[n] \cos^2 \left( \left( \frac{m}{M} - \frac{n}{2} \right) \pi \right) \\
&= x_0 \cos^2 \left( \frac{m}{M} \pi \right) + x_1 \sin^2 \left( \frac{m}{M} \pi \right)
\end{aligned}
$$
For $ x_0 = 1 $ and $ x_1 = -1 $
$$
\begin{aligned}
y_m &=  \cos^2 \left( \frac{m}{M} \pi \right) - \sin^2 \left( \frac{m}{M} \pi \right) \\
&=  \cos \left( \frac{m}{M} 2 \pi \right)
\end{aligned}
$$
I'm going to have to be done with this for a while.  Neat stuff.

Olli, thanks for the bounty points.
This little exercise has deepened my understanding of W-S considerably.  I hope that is true for you and Robert (and others) too.
It is still a precarious foundation though.  I wanted to convince myself that it would work for a sinusoid of any frequency.  To wit:
$$ x[n] = M \cos( \alpha n + \phi ) $$
$$
\begin{aligned}
x(t) &= \sum_{n=-\infty}^{\infty} x[n] \operatorname{sinc}(t-n) \\
&= \sum_{n=-\infty}^{\infty} M \cos( \alpha n + \phi ) \operatorname{sinc}(t-n) \\
&= \sum_{n=-\infty}^{\infty} M \cos( \alpha t + \phi - \alpha( t - n ) ) \operatorname{sinc}(t-n) \\
&= \sum_{n=-\infty}^{\infty} M \left[ \cos( \alpha t + \phi ) \cos( \alpha( t - n ) ) + \sin( \alpha t + \phi ) \sin( \alpha( t - n ) ) \right] \operatorname{sinc}(t-n) \\
&= M \cos( \alpha t + \phi ) \sum_{n=-\infty}^{\infty} \cos( \alpha( t - n ) ) \operatorname{sinc}(t-n) \\
& \qquad \qquad + M \sin( \alpha t + \phi ) \sum_{n=-\infty}^{\infty} \sin( \alpha( t - n ) ) \operatorname{sinc}(t-n) \\
&= M \cos( \alpha t + \phi ) \cos( \alpha( t - t ) ) + M \sin( \alpha t + \phi ) \sin( \alpha( t - t ) ) \\
&= M \cos( \alpha t + \phi ) \cdot 1  + M \sin( \alpha t + \phi )  \cdot 0  \\
&= M \cos( \alpha t + \phi )
\end{aligned}
$$
I seem to have accomplished my goal.  However, there is nothing in this proof that prohibits $\alpha \ge \pi$, though that is a condition for the validity of the theorem.  So, knowing that, you are okay.  If you didn't know that, the formula itself does not reveal it.  To me, that's troubling.

Reply to R B-J:
First off, no where is it stipulated that $x[n]$ must be real.  Even for a real valued function, you don't have to split the Nyquist bin halfsies to get a real interpolation function.  Just pick $g$ to be a multiple of $i$ above.
Suppose you have the function:
$$ z(\tau) = \sum_{k=-L}^{L} c_k e^{ik\tau} $$
Its band limit is $L$ or less.  Every k term, except 0, can be paired up with it's conjugate bin and the sum can be decomposed into a cosine and sine term.
let $ A = \frac{c_k + c_{-k}}{2} $ and $ B = \frac{c_k - c_{-k}}{2} $
$$
\begin{aligned}
c_k e^{ik\tau} + c_{-k} e^{-ik\tau} &= (A+B) e^{ik\tau} + (A-B) e^{-ik\tau} \\
&= 2A \cos(\tau) - i 2B \sin(\tau)
\end{aligned}
$$
For a regular bin, we can only say $X[k] = c_k$ if $k+N>L$, otherwise I have more than one k in the bin and cannot separate them.  At the Nyquist bin $X[k] = c_k + c_{-k}$
Think in terms of degrees of freedom.  For a complex signal, $c_k + c_{-k}$ has four and the Nyquist bin two.  Therefore there are two free.  Just enough to put a complex parameter on the Sine function at the Nyquist frequency.  With a real signal, $c_k + c_{-k}$ has two degrees of freedom and the Nyquist bin value restricts one of those leaving one left over.  Just enough for a real valued parameter times the Sine function to remain a real valued signal.
I showed earlier the translation between not doing halfsies and the consequence on the interpolation function.  Nothing prohibits that and it doesn't increase the bandwidth of the solution one iota.

R B-J asks:
//  //"But we do know A will be zero in the halfsies and W-S reconstructions."// how do you know that?  //
The halfsies is easy.  Without loss of generality, consider the $N=2$ case.
$$ x[n] = [1,-1] $$
$$ \frac{1}{N} X[k] = [0,1] $$
Halfsies on the Nyquist of 1.  Doing an unfurled inverse DFT with split Nyquist:
$$ x[n] = \frac{1}{2} e^{i\pi n} + \frac{1}{2} e^{-i\pi n} = \cos(\pi n) $$
Now allow $n$ to be real, call it $t$ to indicate the change.  This defines an interpolation function (still called $x$).
$$ x(t) = \cos(\pi t) $$
For every other even N, the unnormalized DFT will be (0,0,0,....,N), so the result remains the same.
For the W-S summation, look at the section where "omega" temporarily lived, the "Sinc is the limit of the Dirichlet Kernel" section.  The left side $y_m=Y(w)$ is known to be $ \cos( \pi w ) $.  I even did the specific $N=2$ case after the dependency list.  Just set "M=2" which makes $w = m$.  The limit reached at the end of the second chunk gives your summation.  Just reverse the order of the equation and you get:
$$ \sum_{n=-\infty}^{\infty} (-1)^n \operatorname{sinc}( w - n ) = \cos( \pi w ) $$
The fact that your summation is the limit of something is why proving it differently has been hard.
I think your time reversal argument is good, too.  The sampled points are time reversible on the discrete $n$ scale, but that does not mean the source x(t) is, but it does mean Y(w) is.
P.S.  From now on, when a fresh context can be established, I'm going to use $\tau$ for a $ 0 \to 2\pi $ cycle scale, $t$ to be on the sampling scale ($=n$).


Confusion in CT Fourier Transform Proof I am confused trying to understand the Proof of Fourier Transform from Oppenheim book Signals and Systems. I am pasting the equations directly from the book:

  $$\widetilde{x}(t)=\sum_{k=-\infty}^{+\infty}a_ke^{jk\omega_0t} \tag{4.3} $$
  $$a_k=\frac{1}{T}\int_{-T/2}^{T/2}\widetilde{x}(t)e^{-jk\omega_0t}dt \tag{4.4} $$

  where $\omega_0=2\pi/T$. Since $\widetilde{x}(t)=x(t)$ for $|t| < T/2$, and also, since $x(t)=0$ outside this interval, eq. $(4.4)$ can be
  rewritten as

  $$a_k=\frac{1}{T}\int_{-T/2}^{T/2}x(t)e^{-jk\omega_0t}dt=\frac{1}{T}\int_{-\infty}^{+\infty}x(t)e^{-jk\omega_0t}dt$$

My question is In the last equation on the L.H.S, there is no
$\widetilde{x}(t)$ but it was mentioned in the eq. $(4.4)$
And then we have: $$Ta_k= \left.\frac{2 \sin(\omega T_1)}{\omega}\right\rvert_{\omega =
k\omega_0}$$ That is, with $\omega$ thought of as a continuous variable, the function ($2\sin\omega T_1/ \omega$) represents envelope of $Ta_k$.
Here I want to know what is meant by envelope?I want to add more details from the book exactly " Specifically we have plot $T_0a_k$ rather than $a_k$ and we have also modified the horizontal spacing in each plot. here is the plot "

c) I want to know that how $T_0a_k$ will be drawn in continuous?

@robertbristow-johnson can you please explain exactly the answer?

oh dear.  it's a lot of work and it should be in some textbook.

okay, using your O&S notation:

$$ \tilde{x}(t) = \sum\limits_{k=-\infty}^{+\infty} a_k \ e^{j k \omega_0 t} $$

where

$$ \tilde{x}(t) = \tilde{x}(t+T) \quad \forall t $$

$$ \omega_0 \triangleq \frac{2 \pi}{T} $$

and

$$ a_k = \frac{1}{T} \int\limits_{t_0}^{t_0 + T} \tilde{x}(t) \ e^{-j k \omega_0 t} \ dt $$

and where $t_0$ can be any real value:  $-\infty < t_0 < +\infty$

we choose $t_0$ to be $t_0 = -\frac{T}{2}$ so that

$$ a_k = \frac{1}{T} \int\limits_{-\frac{T}{2}}^{+\frac{T}{2}} \tilde{x}(t) \ e^{-j k \omega_0 t} \ dt $$

That's Fourier Series.

A Riemann sum is a sum of areas of rectangles that, for decently well-behaved functions, converges to the area of an integral in the limit that the rectangles get thinner and thinner:

$$ \int\limits_{t_1}^{t_2} f(t) \ dt  =  \lim_{N \to +\infty} \sum\limits_{n=0}^{N-1} f\left(t_1 + n\frac{t_2 - t_1}{N} \right) \ \frac{t_2 - t_1}{N} $$

for $t_1 < t_2$.

this (Riemann integration) works for most ordinary functions that are decently well-behaved.  for nastier functions, another method of defining integration may have to be used, but we're not dealing with that here.

so do you understand Fourier series and Riemann integration?

now, define

$$ \tilde{x}(t) =
\begin{cases}
x(t), & \text{if } |t| < \frac{T}{2}  \\
x(t - mT), & \text{if } |t| > \frac{T}{2} \text{ and where } m \triangleq \left\lfloor \frac{t}{T} + \frac{1}{2} \right\rfloor
\end{cases} $$

$\lfloor u \rfloor$ is the $\operatorname{floor}(u)$ function that is equal to the largest integer no larger the argument $u$.  (it's always rounding down.)  $\left\lfloor u + \frac{1}{2} \right\rfloor$ is always rounding to nearest.

you can see that $\tilde{x}(t)$ satisfies periodicity:  $ \tilde{x}(t) = \tilde{x}(t+T) \quad \forall t $.  and you can also see that

$$ x(t) = \lim_{T\to +\infty} \tilde{x}(t) \quad \forall t $$

so we're taking a finite section of $x(t)$ and periodically extending that to get $\tilde{x}(t)$.  but as the period gets longer and longer and goes to $\infty$ the periodic extension becomes the same as $x(t)$.

now if you plug in $a_k$ into the Fourier sum:

$$ \tilde{x}(t) = \sum\limits_{k=-\infty}^{+\infty} \frac{1}{T} \int\limits_{-\frac{T}{2}}^{+\frac{T}{2}} \tilde{x}(u) \ e^{-j k \omega_0 u} \ du \ e^{j k \omega_0 t} $$

and let $T \to +\infty$ ...

okay, i got back to this.  (sorry for the delay.)

so we've been saying

$$ \begin{align}
x(t) & = \lim_{T \to \infty} \tilde{x}(t) \\
& = \lim_{T \to \infty} \sum\limits_{k=-\infty}^{+\infty} \frac{1}{T} \int\limits_{-\frac{T}{2}}^{+\frac{T}{2}} \tilde{x}(u) \ e^{-j k \omega_0 u} \ du \ e^{j k \omega_0 t}  \\
& = \lim_{T \to \infty} \sum\limits_{k=-\infty}^{+\infty} \frac{1}{T} \int\limits_{-\frac{T}{2}}^{+\frac{T}{2}} x(u) \ e^{-j k \omega_0 u} \ du \ e^{j k \omega_0 t}  \\
\end{align} $$

(note the subtle difference between the last two lines.  why can we do that?)

and we're saying

$$ \begin{align}
\int\limits_{-\frac{T}{2}}^{+\frac{T}{2}} x(u) \ e^{-j k \omega_0 u} \ du & = \lim_{N \to +\infty} \sum\limits_{n=0}^{N-1} x\left(-\frac{T}{2} + n\frac{+\frac{T}{2} - (-\frac{T}{2})}{N} \right) \ e^{-j k \omega_0 \left(-\frac{T}{2} + n\frac{+\frac{T}{2} - (-\frac{T}{2})}{N} \right)} \   \frac{+\frac{T}{2} - (-\frac{T}{2})}{N} \\
& = \lim_{N \to +\infty} \sum\limits_{n=0}^{N-1} x\left(-\frac{T}{2} + n\frac{T}{N} \right) \ e^{-j k \omega_0 \left(-\frac{T}{2} + n\frac{T}{N} \right)} \   \frac{T}{N} \\
\end{align} $$

putting the two limits together:

$$ \begin{align}
x(t) & = \lim_{T \to \infty} \sum\limits_{k=-\infty}^{+\infty} \frac{1}{T} \int\limits_{-\frac{T}{2}}^{+\frac{T}{2}} x(u) \ e^{-j k \omega_0 u} \ du \ e^{j k \omega_0 t}  \\
& = \lim_{T \to \infty} \sum\limits_{k=-\infty}^{+\infty} \frac{1}{T} \lim_{N \to +\infty} \sum\limits_{n=0}^{N-1} x\left(-\frac{T}{2} + n\frac{T}{N} \right) \ e^{-j k \omega_0 \left(-\frac{T}{2} + n\frac{T}{N} \right)} \   \frac{T}{N} \ e^{j k \omega_0 t}  \\
& = \lim_{T \to \infty} \lim_{N \to +\infty} \sum\limits_{k=-\infty}^{+\infty}   \sum\limits_{n=0}^{N-1} x\left(-\frac{T}{2} + n\frac{T}{N} \right) \ e^{-j k \omega_0 \left(-\frac{T}{2} + n\frac{T}{N} \right)} \   \frac{T}{N} \ e^{j k \omega_0 t} \frac{1}{T}  \\
& = \lim_{N \to +\infty} \sum\limits_{n=0}^{N-1} \lim_{T \to \infty}  \sum\limits_{k=-\infty}^{+\infty}    x\left(-\frac{T}{2} + n\frac{T}{N} \right) \ e^{-j k \omega_0 \left(-\frac{T}{2} + n\frac{T}{N} \right)} \   \frac{1}{N} \ e^{j k \omega_0 t}  \\
\end{align} $$

(answer not done yet, still being developed.)


Find gaussian in data I have some time series data (blue line lower plot) with a gaussian added into it (red line upper plot). I would like to estimate where the gaussian is in the data, and the height of the gaussian (these parameters can be bounded).
Is there an algorithm that can estimate and find where the gaussian is in the data (blue line, upper plot)?
Or a better question might be, is there a way I could come up with a probability curve for finding a gaussian in the data if I don't know what height or it's location (mean).
The model would be:
$$y(x)=A c_0 e^{-A c_0 \frac{(x-x_0)^2}{x_0}}$$
or
$$y(x)=\frac{A^2}{{(1+e^{-A c_0 \frac{(x-x_0)^2}{x_0}})}^2} c_0 e^{-A c_0 \frac{(x-x_0)^2}{x_0}}$$

Note: this is just a simulated gaussian, normally I have just the blue line in the upper plot.
Explanation of upper plot:
Red line, original data
Blue line, original data with gaussian added (the gaussian is in the lower plot has been added to the blue line in the upper plot)
I only have the blue line in the upper plot, I am trying to detect the location in the x direction of the gaussian (where it peaks parameter x_0) and maybe the height A
The gaussian is my signal, the blue line in the upper plot is my noise
the redline is an example of noise + gaussian and a typical time series that I work with.

Still not making any sense to me.  This is what I have parsed from what you are saying:
Red Line = Noise + Gaussian
Blue Line = Red Line + Lower Gaussian
Blue Line = Noise + Gaussian + Lower Gaussian
Blue Line = Noise
Red Line = Original
Blue Line = What you got
Clear as mud.  What you got isn't what is original.  You arbitrarily chose a spot to add a contrived Gaussian to the red line to get the blue line, yet that is what you started with.
You call the red line a gaussian + noise, but it doesn't include the gaussian you added.
Let's start with this.  Show the data set you are starting with.  The original data, the time series as you will measure it.
I'll take another look in the morning, and probably delete this answer, but it was way too long for a comment.
In the meanwhile, maybe a look at this:  Estimate peak width from a vector that is a superposition of unknown number of identical Gaussian peaks with different heights?


Identifying steps in walking accelerometer I have some recorded walking data from an android phone and want to determine when someones foot hits the ground using the acceleration data. I'm making an assumption here that the person will "bob" vertically as they walk (i.e. when their foot pushes off their body will rise, and as a foot is approaching the ground their body will fall)

Given that, I have taken the $X$-axis acceleration data (which corresponds to the persons vertical axis), removed gravity, and filtered out the noise giving me something like this:

My question is: where in this signal does a foot hit the ground?

I have 3 options:

Would it be all of the valleys in the signal?

This would be correct if I were measuring position information but acceleration is a bit different. I think the marked points are actually the moments of maximum deceleration (or, acceleration towards the ground), so at these times the body is still in motion and not stationary

Would it be the points where acceleration is 0?

0 acceleration would mean the person/device is moving at a constant speed, which I believe would be the case upon foot landing because I'm assuming the speed is 0 and therefore not accelerating. But this would also be true on the "rising" part of the walk cycle - as a person pushes off the ground they will accelerate in +X, but once they reach the peak of the rise they would no longer be accelerating and you would probably see a 0 again

Would it be the 0 acceleration points, but only after a valley?

A valley would, I think, represent point of max acceleration towards the ground (i.e. theyre on the way down), so a 0 acceleration value after a valley would presumably be once theyve stopped accelerating downwards and the foot has landed?

I'm thinking the answer is #3 but wanted to see if my logic is correct, or whether my interpretation of acceleration graphs is completely wrong

EDIT: this data is recorded using a chest harness carrying the phone, so the phone is not on any of the limbs. The main assumption being made is that the torso will move up and down with each step. Also, in the phones coordinate system positive X is pointing towards the persons head

EDIT2: Using a low pass butterworth filter, instead of a bandpass, to eliminate any high frequency noise. The blue vertical lines are where I visually observe the foot hitting the ground so I can use it as a guide (even though it will be off by a few milliseconds):

And here is the filtered signal (red) over the raw signal (black):

It would be simple to make an acquisition of a single step and identify the phases of the step. For that, you shall consider the geometry of your setup, and therefore answer these:

The accelerometer has (or not) a relative speed to the ground equal to ZERO?
The accelerometer measurement axis is aligned to gravity?

My guess is that, when one foot is on the ground, the accelerometer shall measure only gravity for a short period. Therefore, I would look for segments where the jerk is close to ZERO in the vertical axis, AND the magnitude of it is close to -1g.

EDIT:
All this new information you have given us, like sensor close to chest, gravity, alignment, etc, demands new dynamic models to represent the measurements you have there. But any model, with respective solution, can have counter-examples to demonstrate that the model is inadequate at all. For example, if the user carries a heavy rucksack, there is a new model to measure the acceleration close to the neck. This is different from the acceleration tied to the feet. That said, the body shall damp the actual acceleration at the feet to a filtered and delayed measurement at the neck. Nevertheless, the readings shall have the same low frequency components as at the feet. This is more less the same I see at your graph, the one with blue lines. The blue lines are as periodic as the low freq. components of the red line, but for a phase. You should consider to review your necessity of precise instant of touch on the floor in accordance with your application.


How does a FIR equiripple filter behave close to DC? I am multiplying two sine waves with the same frequency (f), but might have a phase difference 0 <-> 90 degrees. The product is a two frequency sinusoidal wave with f1 = f-f = DC and f2 = f+f = 2f. I now want to filter out the 2f component. I am currently sampling at Fs = 32 x f. I get a 40th order equiripple FIR filter in Matlab which is ok, but if I later am to increase the sampling rate to 64xf, 128xf etc the filter order gets very high.

I assume however that f1 will not be 0 Hz exactly but somewhere close to 0 Hz. Over this region the ripple in the passband might be unwanted. Can the ripple change that much over such a narrow range? (within whatever limit one has set as the maximum ripple in the passband).

I am using the Filter Design & Analysis tool in Matlab and saw that it does not seem like an equiripple filter has any ripple close to DC (0<->1Hz)? Is this correct?

I need a filter that can attenuate the 2f component with 60dB and need as little computing power as possible. (working on a MCU without FPU using fixed-point).

Is there any other filter that can do a better job than equiripple for what I described? I am open for other types of filters as IIR if someone can convince me of their advantages.

If you want to know the exact DC gain of your filter take the sum of the coefficients. This is true for an equiripple filter as well as any other FIR filter.  There is no such thing as ripple at DC, because "ripple" is a term used to describe a magnitude variation over some range of frequencies.

If what you are looking for is a flat passband then equiripple can fit the bill.  When using the fdatool in Matlab, notice that the equiripple FIR filter allows you to specify a parameter Apass.  This paramter lets you adjust the maximum allowable amplitude of the ripples in the passband as shown below.  By default it is set to 1, which means that you will have a maximum ripple of 1 dB in the passband.

If you are not constrained to FIR filters, then you might want to consider an IIR filter, depending on your needs.  You can usually achieve the desired result with a lower filter order when using IIR.  This means less computation and less delay.


Why Butterworth filter always starts my signal from from zero mark (amplitude)? I am trying to filter my force signal using Butterworth low pass filter.

so this is the signal, I  have to filter (noise frequency is around 33 Hz). So when I use Butterworth filter with following command in MATLAB (with sixth order):

[B,A]=butter(6,5/100,'low');
forceS=filter(B,A,forceS);

the result I get is

So the problem is the filtered signal start from zero mark instead of around 1.5. How can i start the filtered signal from the same value as unfiltered signal has?

Thanks

Filters do have a delay (a lag) since they do not act immediately on your signal. Also all samples before the time 0 are zeros, thus in general you will start from the "zero mark", as you said (just imagine your filter equation with all zeros).

There are ways to make a filter to have a zero lag. It is done by so called zero-phase filtering, also known as forward-backward filtering. The way you do it is by filtering your signal twice - in forward and in backward direction. Obviously this can work only for offline applications. In MATLAB you can do it very easily using the filtfilt function.

Please find below a code corresponding to your case.

fs = 100;
T = 2;
t = 0:1/fs:T;
f = 2;
s = 1.5*cos(2*pi*f*t) + 0.7*sin(2*pi*33*t);

[B,A]=butter(6, 5/fs,'low');
S=filter(B, A, s);
Sff=filtfilt(B, A, s);

plot(t, s)
hold on
plot(t, S)
plot(t, Sff)
grid on
legend({'Original signal', 'Signal filtered using filter', 'Signal filtered using filtfilt'})

That yields what I suspect you want:


Question on sketching the amplitude spectrum of signals I can't seem to understand the sketches of the amplitude spectrums provided in my tutorial solutions.

For this part , I understand that X(f) is basically sampling (T/Tp)sinc(Tf). But from what I understand , the amplitude spectrum is a plot of |X(f)| against f. With that being said, shouldn't the sketch in the solution of the sinc function have no portions that go below the x axis? Since |sinc| means that you flip the negative portions up.

Same problem as question 1, why are there negative magnitudes for the frequency components 12pi and -12pi? Given x(t) , X(f) would just be a sum of dirac delta functions , and any negative signs can be changed to a phase. As such , there shouldn't be any negative magnitudes in the amplitude spectrum.

I'm suspecting that maybe there's a difference between a continuous-frequency spectrum and an amplitude spectrum, but question 2 tells me otherwise.

Thanks for the help!

It's really very simple, but should have been clarified in your tutorial. It may turn out that a signal's Fourier Transform is real (I'll let you figure out when that happens). In that case, the spectrum $X(f)$ can be drawn in a single plot using positive and negative numbers, without having to split it up into magnitude and phase spectrums.

You can, of course, represent a real number in polar form. If the number is positive, its phase is zero, and if it is negative, its phase is $\pi$. What this means is that, whether you plot the real Fourier transform using negative and positive numbers, or you split it into magnitude and phase plots, you're conveying the exact same information.